Mon. Not. R. Astron. Soc. 000, 000–000 (0000) Printed 18 August 2009 (MN LATEX style file v2.

2)

Chaotic Transport and Chronology of Complex Asteroid Families

B. Novaković1⋆ , K. Tsiganis2⋆ and Z. Knežević1⋆

1 Astronomical Observatory, Volgina 7, 110 60 Belgrade 38, Serbia
2 Section of Astrophysics, Astronomy & Mechanics, Department of Physics, Aristotle University of Thessaloniki, GR 54 124 Thessaloniki, Greece

18 August 2009

ABSTRACT
We present a transport model that describes the orbital diffusion of asteroids in chaotic regions
of the 3-D space of proper elements. Our goal is to use a simple random-walk model to study
the evolution and derive accurate age estimates for dynamically complex asteroid families. To
this purpose, we first compute local diffusion coefficients, which characterize chaotic diffu-
sion in proper eccentricity (ep ) and inclination (Ip ), in a selected phase-space region. Then, a
Monte-Carlo-type code is constructed and used to track the evolution of random walkers (i.e.
asteroids), by coupling diffusion in (ep ,Ip ) with a drift in proper semi-major axis (ap ) induced
by the Yarkovsky/YORP thermal effects. We validate our model by applying it to the family
of (490) Veritas, for which we recover previous estimates of its age (∼ 8.3 Myr). Moreover,
we show that the spreading of chaotic family members in proper elements space is well repro-
duced in our randomk-walk simulations. Finally, we apply our model to the family of (3556)
Lixiaohua, which is much older than Veritas and thus much more affected by thermal forces.
We find the age of the Lixiaohua family to be 155 ± 36 Myr.
Key words: celestial mechanics, minor planets, asteroids, methods: numerical

1 INTRODUCTION young families (i.e. age < 10 Myr), is to integrate the orbits of the
family members backwards in time, until the orbital orientation an-
As first noted by Hirayama (1918), asteroids can form prominent
gles cluster around some value. As such a conjunction of the orbital
groupings in the space of orbital elements. These groups, nowadays
elements can occur only immediately after the disruption of the par-
well-known as asteroid families, are believed to have resulted from
ent body, the time of conjunction indicates the formation time. The
catastrophic collisions among asteroids, which lead to the ejection
method was successfully applied by Nesvorný et al. (2002, 2003)
of fragments into nearby heliocentric orbits, with relative velocities
to estimate the ages of the Karin cluster (5.8±0.2 Myr) and of the
much lower than their orbital speeds. To date, several tens of fami-
Veritas family (8.3±0.5 Myr). This method is however limited to
lies have been discovered across the whole asteroid Main Belt (e.g.
groups of objects residing on regular orbits.
Bendjoya & Zappalà 2002; Nesvorný et al. 2006). Also, families
For older families (i.e. age > 100 Myr), one can make use
have been identified among the Trojans (Milani 1993; Beaugé &
of the fact that asteroids slowly spread in semi-major axis due to
Roig 2001), and most recently, proposed to exist in the Transnep-
the action of Yarkovsky thermal forces (Farinella & Vokrouhlický
tunian region (Brown et al. 2007). Studies of asteroid families are
1999). As small bodies drift faster than large bodies, the distribu-
very important for planetary science. Families can be used, e.g. to
tion of family members in the (ap , H) plane – where ap is the
understand the collisional history of the asteroid Main Belt (Bottke
proper semi-major axis and H is the absolute magnitude – can be
et al. 2005), the outcomes of disruption events over a size range in-
used as a clock. An improved version of this method, which ac-
accessible to laboratory experiments (e.g. Michel et al. 2003; Durda
counts for the initial ejection velocity field and the action of YORP
et al. 2007), to understand the mineralogical structure of their par-
thermal torques, has been successfully applied to several families
ent bodies (e.g. Cellino et al. 2002) and the effects of related dust
by Vokrouhlický et al. (2006a) and Vokrouhlický et al. (2006b).
“showers” on the Earth (Farley et al. 2006). Obtaining the relevant
Again, it is not straightforward to apply this method to families lo-
information is, however, not easy. One of the main complications
cated in the chaotic regions of the asteroid belt.
arises from the fact that the age of a family is, in general, unknown.
Thus, accurate dating of asteroid families is an important issue in Milani & Farinella (1994) suggested that asteroid families,
the asteroid science. which reside in chaotic zones, can be approximately dated by
A number of age-determination methods have been proposed chaotic chronology. This method is based on the fact that the age
so far. Probably the most accurate procedure, particularly suited for of the family cannot be greater than the time needed for its most
chaotic members to escape from the family region. In its original
form, this method provides only an upper bound for the age. Re-
⋆ E-mail: [email protected] (BN); [email protected] (KT); cently, Tsiganis et al. (2007) introduced an improved version of this
[email protected] (ZK) method, based on a statistical description of transport in the phase

c 0000 RAS
2 B. Novaković et al.
space. Tsiganis et al. (2007) successfully applied it to the family latter can be obtained by adding in the fictitious bodies, selected in
of (490) Veritas, finding an age of 8.7±1.7 Myr, which is statisti- such a way that they occupy the same region of the proper elements
cally the same as that of 8.3±0.5 Myr, obtained by Nesvorný et al. space as the real family members. Since we wish to study just these
(2003)1 . Despite these improvements, the chaotic chronology still local diffusion properties and the effect of the use of variable coef-
suffered from two important limitations. It did not account for the ficients in our chronology method, we are going to follow here this
variations in diffusion in different parts of a chaotic zone, which can strategy.
significantly alter the distribution of family members (i.e. the shape The 3-D space of proper elements, occupied by the selected
of the family). Moreover, it did not account for Yarkovsky/YORP family, is divided into a number of cells. Then, in each cell, the dif-
effects, thus being inadequate for the study of older families. fusion coefficientsqare calculated for bothqrelevant action variables,
a a
In this paper we extend the chaotic chronology method, by namely J1 ∼ 21 aJp e2p and J2 ∼ 12 aJp sin2 Ip (aJ denotes
constructing a more advanced transport model, which alleviates Jupiter’s semi-major axis, ep the proper eccentricity and Ip the
the above limitations. We first use the Veritas family as bench- proper inclination of the asteroid). This is done by calculating the
mark, since its age can be considered well-defined. Local diffu- time evolution of the mean squared displacement h(∆Ji )2 i (i=1,2)
sion coefficients are numerically computed, throughout the region in each action, the average taken over the set of bodies (real or fic-
of proper elements occupied by the family. These local coefficients titious) that reside in this cell. The diffusion coefficient is then de-
characterize the efficiency of chaotic transport at different loca- fined as the least-squares-fit slope of the h(∆Ji )2 i(t) curve, while
tions within the considered zone. A Monte-Carlo-type model is the formal error is computed as in Tsiganis et al. (2007).
then constructed, in analogy to the one used by Tsiganis et al. The simulation of the spreading of family members in the
(2007). The novelty of the present model is that it assumes vari- space of proper actions and the determination of the age of the fam-
able transport coefficients, as well as a drift in semi-major axis due ily is done using a Markov Chain Monte Carlo (MCMC) technique
to Yarkovsky/YORP effects, although the latter is ignored when (e.g. Gentle 2003; Berg 2004). At each step in the simulation the
studying the Veritas family. Applying our model to Veritas, we find random walkers can change their position in all three directions, i.e.
that both (a) the shape of its chaotic component and (b) its age the proper semi-major axis ap and the two actions J1 and J2 . Al-
are correctly recovered. We then apply our model to the family of though no macroscopic diffusion occurs in proper semi-major axis,
(3556) Lixiaohua, another outer-belt family but much older than the random walker can change its ap value due to the Yarkovsky
Veritas and hence much more affected by the Yarkovsky/YORP effect, while the changes in J1 and J2 are controlled by the local
thermal effects. We find the age of the Lixiaohua family to be values of the diffusion coefficients. In the case of normal diffu-
∼ 155 Myr. sion the transport properties in action space are determined by the
We note that, depending on the variability of diffusion coeffi- solution of the Fokker-Planck equation (see Lichtenberg & Lieber-
cients in the considered region of proper elements, this new trans- man (1983)). The MCMC method is in fact equivalent to solving a
port model can be computationally much more expensive than the discretized 2-D Fokker-Planck equation with variable coefficients,
one applied in Tsiganis et al. (2007). This is because, if the val- combined here with a 1-D equation for the Yarkovsky-induced dis-
ues of the diffusion coefficients vary a lot accross the considered placement in ap . The latter acts as a drift term, contributing to the
region, one would have to calculate them in many different points. variability of diffusion in J1 and J2 .
However, even so, this computation needs to be performed only The rate of change of ap due to the Yarkovsky thermal force, is
once. Then, the random-walk model can be used to perform multi- given by the following equation (e.g. Vokrouhlický 1999; Farinella
ple runs at very low cost, e.g. to test different hypotheses about the & Vokrouhlický 1999):
original ejection velocities field or about the physical properties of
the asteroids. On the other hand, for “smooth” diffusion regions in da
= k1 cos γ + k2 sin2 γ (1)
which the coefficients only change by a factor of 2-3 across the con- dt
sidered domain, the model can be simplified. In such regions, the where the coefficients k1 and k2 depend on parameters that de-
age of a family can be accurately determined even by assuming an scribe physical and thermal characteristics of the asteroid and γ de-
average (i.e. constant over the entire region) diffusion coefficient, notes the obliquity of the body’s spin axis. For km-sized asteroids,
as we show in Section 3. the drift rate is inversely proportional to their radius. This simpli-
fied Yarkovsky model assumes that the asteroid follows a circular
orbit, and thus linear analysis can be used to describe heat diffusion
across the asteroid’s surface. The obliquity of the spin axis and the
2 THE MODEL
angular velocity of rotation (ω) of the asteroid are subject to ther-
Our study begins by selecting the target phase-space region. This is mal torques (YORP) that change their values with time, according
done by identifying the members of an asteroid family crossed by to the following equations:
resonances, from a catalog of numbered asteroids. Apart from the
dω dγ g(γ)
largest Hirayama families, for the other smaller and more compact = f (γ) , = (2)
dt dt ω
ones, in the current catalog one typically finds up to several hun-
dred members. Thus the chaotic component of the family consists (e.g. Vokrouhlický & Čapek 2002; Čapek & Vokrouhlický 2004),
of a few tens to a few hundreds of asteroids. Although this may be where the functions f and g describe the mean strength of the
adequate to compute the average values of the diffusion coefficients YORP torque and depend on the asteroid’s surface thermal con-
(as in Tsiganis et al. 2007), a detailed investigation of the local dif- ductivity (Čapek & Vokrouhlický 2004).
fusion characteristics requires a much larger sample of bodies. The The length of the jump in ap that a random-walker under-
takes at each time-step dt in the MCMC simulation, is determined
by equations (1)-(2), in their discretized form. Of course, a set
1 Of course, Tsiganis et al. (2007) and Nevorný et al. (2003) used a chaotic of values of the physical parameters must be assigned to each
and a regular subsets of the family, respectively. body. As the majority of the Veritas family members are of C-

c 0000 RAS, MNRAS 000, 000–000
 Chaotic Transport and Chronology of Complex Asteroid Families 3
type (Di Martino et al. 1997; Mothé-Diniz et al. 2005), while the 3 THE RESULTS
Lixiaohua family members seem to be C/X-type (Lazzaro et al.
In this section we use our model to study the evolution of the
2004; Nesvorný et al. 2005), the following values for these pa-
chaotic component of two outer-belt asteroid families: (490) Ver-
rameters (Brož et al. 2005; Brož 2006) are adopted: thermal con-
itas and (3556) Lixiaohua. In both cases, a number of mean motion
ductivity K = 0.01 − 0.5 [W (m K)−1 ], specific heat capacity
resonances (MMR) cut-through the family, such that a significant
C = 1000 [J (K kg)−1 ], and the same value for surface and bulk
fraction of members follow chaotic trajectories. On the other hand,
density ̺ = 1500 [kg m−3 ]. In Tedesco et al. (2002), the geomet-
their ages differ significantly, according to previous estimates. In
ric albedos (pv ) of several Veritas and Lixiaohua family members
this respect, the Yarkovsky effect can be neglected in the study of
are listed, yielding a mean pv = 0.068 ± 0.018 for Veritas and
Veritas, but not in the study of Lixiohua.
pv = 0.049 ± 0.027 for Lixiaohua. The rotation period, P , is cho-
We begin by performing an extensive study of the local dif-
sen randomly from a Gaussian distribution peaked at P = 8 h,
fusion properties in the chaotic region of the Veritas family. Then,
while the distribution of initial obliquities, γ, is assumed to be uni-
the MCMC model is used to simulate the evolution of the chaotic
form2 . To assign the appropriate values of absolute magnitude H
members and to derive an estimate of the age of the family. The
to each body, we need to have an estimate of the cumulative distri-
results are compared to the ones given by the model of Tsiganis et
bution N (< H) of family members. A power-law approximation
al. (2007). Finally, we apply the MCMC model to the Lixiaohua
is used (e.g. Vokrouhlický et al. 2006b)
family and derive estimates of its age, for different values of the
N (< H) ∝ 10βH (3) Yarkovsky-related physical parameters.

where β depends on the considered interval for H; e.g. for the Ver-
3.1 The Veritas family
itas family, we find β = 0.74 ± 0.03 for H ∈[11.5,13.5] and
β = 0.23 ± 0.03 for H ∈[13.5,15.5]. Having the values of H The Veritas family is a comparatively small and compact outer-belt
and pv , the radius R of a body can be estimated, using the relation family, spectroscopically different from the background population
(e.g. Carruba et al. 2003) of asteroids. In terms of dynamics, it occupies a very interesting
−H
and complex region, crossed by several mean motion resonances.
10 5 Application of the Hierarchical Clustering Method (HCM) Zappalà
R (km) = 1329 √ (4)
2 pv et al. (1995) to the AstDys catalog of synthetic proper elements
(numbered asteroids http://hamilton.unipi.it/astdys as of XXXX),
At each time-step in the MCMC simulation, a random-walker yields 409 family members, for a velocity cut-off of vc = 40 m s−1
p a jump in J1 and J2 , whose length is given by ∆Ji
suffers as in Tsiganis et al. (2007).
=µ D(Ji )dt/2 (i = 1, 2), where µ is a random number from a Although the family appears now to extend beyond ap =
Gaussian distribution (Tsiganis et al. 2007). Since the values of the 3.18 AU (see Fig. 1), the main dynamical groups remain practically
diffusion coefficients D(Ji ) vary in space, the maximum allowable the same (see Tsiganis et al. 2007, for a detailed description of the
jump, for a given dt, changes from cell to cell. In our simulations, groups). Since the scope of this paper is to present a refined trans-
the values of D(J1 ) and D(J2 ) used for each body, are given by: port model, we will briefly describe here only the main relevant
features, referring to a forthcoming paper for a renewed analysis of
dR dL
D(Ji ) = DLi + DRi , (5) the Veritas family itself.
dL + dR dL + dR
The main chaotic zone, where appreciable diffusion in proper
where dL and dR denote the distances of the random-walker from elements is observed, is located around 3.174 AU (Fig. 1) and is
the two nearest nodes (left and right) and DL and DR denote the associated with the action of the (5,-2,-2) three-body mean motion
corresponding values of the diffusion coefficients at these nodes.3 resonance (MMR); see Fig. 2 for the typical short-term evolution
For a correct determination of the age of the family, the ran- of such a resonant asteroid. The family members that reside in this
dom walkers have to be placed initially in a region, whose size is resonance can disperse over the observed range in ep and sin Ip on
as close as possible to the size that the real family members oc- a ∼ 10 Myr time-scale. This is exactly the group of bodies (group
cupied, immediately after the family-forming event. This is in fact A) that was used by Tsiganis et al. (2007), to compute the age of
a source of uncertainty for our model. In our calculations we as- Veritas.
sumed the initial spread of the family in (ap , ep ) and (ap , sin Ip )
to be accurately represented by a Gaussian equivelocity ellipse (see
3.1.1 The local diffusion coefficients
Morbidelli et al. 1995), computed such that (i) the spread in ap
of the whole family and (ii) the spread in ep and sin Ip of family As the number of bodies in group A is small, we need to generate a
members that follow regular orbits is well reproduced. uniform distribution of fictitious bodies, in order to compute local
diffusion coefficients across the observed range in (ep , sin Ip ). For
this reason we start by selecting ∼ 25, 000 initial conditions (fic-
titious bodies), covering the same region as the real Veritas family
2
members, in the space of osculating elements. We note that the ac-
According to Paolicchi et al. (1996) a size-rotation relation suggests that tual number of bodies used in the calculations of the coefficients
smaller fragments are rotating somewhat faster than the larger ones (see
is much smaller than that (see below). The orbits of the fictitious
also Donnison 2003). Also, Kryszczyńska et al. (2007) claim that the poles
are not isotropically distributed, as general theoretical considerations may
bodies are integrated for a time-span of 10 Myr, using the ORBIT9
predict. These facts are not considered here, but could become important integrator (version 9e), in a model that includes the four major plan-
when studying very old families. ets (Jupiter to Neptune) as perturbing bodies. The indirect effect of
√
3 A geometric mean D(J )= D D could be used instead of an arith-
i L R the inner planets is accounted for by applying a barycentric correc-
metic one; we actually found negligible differences. tion to the initial conditions. This model is adequate for studying

c 0000 RAS, MNRAS 000, 000–000
4 B. Novaković et al.

0.075 3 3 -2 5 -2 -2 7 -7 -2

3.177
0.07
3.176
0.065 3.175

ap [AU]
ep

3.174
0.06
3.173

0.055 3.172

3.171
50 100 150 200 250 300
0.05
t [kyr]
3.15 3.155 3.16 3.165 3.17 3.175 3.18 3.185 3.19
ap [AU]

0.168 3 3 -2 5 -2 -2 7 -7 -2
0.166
350
0.164 300
0.162 250
σ(5,-2,-2) [deg]
sin Ip

0.16 200

0.158 150
100
0.156
50
0.154
0
0.152 50 100 150 200 250 300
3.15 3.155 3.16 3.165 3.17 3.175 3.18 3.185 3.19 t [kyr]
ap [AU]

Figure 1. Distribution of the real Veritas family members in the (ap ,ep ) Figure 2. The time evolution of the mean semi-major axis a (top) and the
and (ap ,sinIp ) planes. The superimposed ellipses represent equivelocity resonant angle (bottom), for a fictitious body inside the (5, -2, -2) three-body
curves, computed according to the equations of Gauss (e.g. Morbidelli et MMR. Note the correlation between the two quantities. Temporary capture
al. 1995), for a velocity of 40m s−1 , true anomaly f = 30◦ and argument at the resonance border occurs when a is at maximum (resp. minimum) and
of pericentre ω = 150◦ . The vertical dashed lines represent approximate σ(5,−2,−2) circulates in a positive (resp. negative) sense.
borders of the main three-body MMRs, as indicated by the corresponding
labels.

number of cells in ap and a small number of cells in ep and sin Ip ,

outer-belt asteroids. Note that the integration time used here is in except in the (5,-2,-2) region. The efficiency of the computation is
fact longer than the known age of the Veritas family. This is done in improved if we use a moving-average technique (i.e. overlapping
order to study the convergence of the computation of the diffusion cells), instead of a large number of static cells, because in the latter
coefficients, with respect to the integration time-span. case we would need a significantly larger number of fictitious bod-
For each body mean elements are computed on-line, by ap- ies. We selected the size of a cell in each dimension as well as an ap-
plying digital filtering, and proper elements are subsequently com- propriate step-size, by which we shift the cell through the family, as
puted according to the analytical theory of Milani & Knežević follows: for ap , the cell-size was ∆ap = 5×10−4 AU and the step-
(1990, 1994). Synthetic proper elements (Knežević & Milani 2000) size 10−4 AU; for J1 the cell-size was ∆J1 = 10−4 and the step-
are also calculated, for comparison and control. Since the mapping size 4 × 10−5 ; finally, for J2 , the cell-size was ∆J2 = 2 × 10−4
from osculating to proper elements is not linear, the distribution of and the step-size 10−4 . Thus, the total number of (overlapping)
the fictitious bodies in the space of proper elements is not uniform, cells used in our computations was 2, 196.
which can complicate the statistics. A smaller sample of ∼ 10, 000 The time evolution of the mean squared displacement in J1
bodies, with practically uniform distribution in proper elements, is and J2 is shown in Fig. 3, for a representative cell in the (5,-2,-2)
therefore chosen. Thus, our statistical sample, on which all com- resonance. The evolution is basically linear in time, as it should be
putations are based, is in fact ∼ 25 times larger than the actual for normal diffusion. The slope of the fitted line defines the value
population of the family. of the local diffusion coefficient. When performing such computa-
As explained in the previous section, we computed local dif- tions, one needs to know (i) what is the shortest possible integration
fusion coefficients, by dividing the space occupied by the Veritas time-span, and (ii) what is the smallest possible number of fictitious
family in a number of cells. Our preliminary experiments suggested bodies per cell, for which reliable values of the coefficients can be
that, while a large number of cells is needed to accurately represent obtained. For several different groups of fictitious bodies (i.e. dif-
the dependence of the coefficients on ap , the same is not true for ferent cells), we calculated the diffusion coefficients using differ-
ep and sin IP , except for the wide chaotic zone of the (5, −2, −2) ent values of the integration time-span, between 1 and 10 Myr. Our
MMR. Thus, we decided to follow the strategy of using a large results suggest, that an integration time of ∼ 4 Myr is sufficient

c 0000 RAS, MNRAS 000, 000–000
 Chaotic Transport and Chronology of Complex Asteroid Families 5

1.6e-007

1.4e-007
5 -2 -2 ap = 3.174 AU
2.5e-014
1.2e-007
2e-014
1e-007
1.5e-014
<(∆J1)2>

D(J1)
8e-008
1e-014
6e-008
5e-015
4e-008
0
0 1 2 3 4 5 6 7 8 9 10 11
2e-008
tint [Myr]

0
0 2e+006 4e+006 6e+006 8e+006 1e+007
t [yr]
1.6e-007

1.4e-007
5 -2 -2 ap = 3.174 AU
2.5e-014
1.2e-007
2e-014
1e-007
1.5e-014
<(∆J2)2>

D(J2)

8e-008
1e-014
6e-008
5e-015
4e-008
0
0 1 2 3 4 5 6 7 8 9 10 11
2e-008
tint [Myr]

0
0 2e+006 4e+006 6e+006 8e+006 1e+007
t [yr]

Figure 3. The mean squared displacement in J1 and J2 as a function of Figure 4. The values of the diffusion coefficients D(J1 ) and D(J2 ) as a
time, for a cell inside the (5,-2,-2) MMR. The evolution is basically linear function of the integration time span. Each curve on these graphs corre-
in time, as can be seen from the respective fit given on each graph, with a sponds to a different cell.
superimposed small amplitude oscillation of a period of sim5 Myr.

local minima are D(J1 ) = (0.73±0.01)×10−14 yr−1 , and D(J2 )

to obtain reliable values, as shown in Fig. 4. This saturation time = (1.16±0.02)×10−14 yr−1 . This form of dependence in ap is in
is about half the known age of the Veritas family. Hence, for this agreement with the results of Nesvorný & Morbidelli (1999), where
case, computing diffusion coefficients is practically as expensive it was shown that the dynamics in this resonance are similar to those
as studying the evolution of the family by long-term integrations. of a modulated pendulum (see also Morbidelli (2002)), with an is-
However, the saturation time is related to the resonance in ques- land of regular motion persisting at the center of the resonance.
tion and not to the age of the family, which could be much longer. However, the size of the island decreases with decreasing ep , so
Thus, as a matter of principle, the computational gain can become that the diffusion coefficients may decrease as the resonance center
important when dealing with much older families. Resonances of is approached, but do not go to zero.
similar order are characterized by similar Lyapunov and diffusion As also seen in Fig. 5, the (5,-2,-2) is by far the most im-
times (see (Murray & Holman 1997)), and thus similar computation portant resonance in the Veritas region, associated with the widest
time-spans (i.e. a few Myr) should be used, for various resonances chaotic zone. Bodies inside this resonance exhibit a complex be-
throughout the belt. havior, as already noted by Knežević & Pavlović (2002). Most res-
The dependence of the diffusion coefficients on the number onant bodies show oscillations in proper semi-major axis around ap
of bodies considered in each cell was also tested. For a number of = 3.174 AU, but some are temporarily trapped near the resonance’s
different cells, we calculated the coefficients, using from 10 to 100 borders (see also Fig. 2), at ap = 3.172 AU or ap = 3.1755 AU. This
bodies for the computation of the corresponding avergages. Our “stickiness” can be important, as it can affect the diffusion rate. In
results suggest that at least 50 bodies per cell are needed, for an fact, we find that the ap values of these bodies shift towards the
accurate computation. resonance borders, where slower diffusion rates are also measured.
The values of the diffusion coefficients in J1 and J2 , along The diffusion properties are quite different at the (3,3,-2) (lo-
with their formal errors, are given as functions of ap in Fig. 5. The cated at ap ≈ 3.168 AU) and (7,-7,-2) (at ap ≈ 3.18 AU) MMRs,
largest diffusion rate is measured in the (5,-2,-2) MMR, which cuts which are of higher order in eccentricity with respect to the (5,-2,-
through the family at ap ∼ 3.174 AU. Both coefficients increase as 2) MMR (see Nesvorný & Morbidelli 1999). The values of D(J1 )
the center of the resonance is approached, but show local minima for the (3,3,-2) resonance (see Fig. 5a) are also increasing as the
at ap ≈3.174 AU, which is approximately the location of the cen- center of the resonance is approached, but no local minimum is
ter. The maximum values are D(J1 ) = (1.28±0.02)×10−14 yr−1 , seen near the center of the resonance, at least in this resolution.
and D(J2 ) = (1.36±0.02)×10−14 yr−1 , while the values at the The maximum values are only D(J1 ) = (8.30±0.09)×10−16 yr−1

c 0000 RAS, MNRAS 000, 000–000
6 B. Novaković et al.
1.4e-014

3 3 -2 7 -7 -2
1.2e-014
1e-015 1e-015
1.8e-014
1e-014 5e-016 5e-016 0.165
1.6e-014
D(J1) [yr-1]

8e-015 0 0
3.167 3.168 3.169 3.17 3.1795 3.18 3.1805

6e-015 0.16 1.4e-014

sin Ip
4e-015
1.2e-014
2e-015 0.155
1e-014
0
3.165 3.17 3.175 3.18 3.185
ap [AU]
0.15 8e-015
1.4e-014
6e-015
3 3 -2 7 -7 -2
1.2e-014 0.055 0.06 0.065 0.07
4e-017 4e-017
ep
1e-014
2e-017 2e-017
D(J2) [yr-1]

8e-015 0 0
3.167 3.168 3.169 3.17 3.1795 3.18 3.1805

6e-015

4e-015 2e-014
0.165
2e-015
1.8e-014
0
3.165 3.17 3.175 3.18 3.185 0.16 1.6e-014
ap [AU]
sin Ip

1.4e-014
Figure 5. The local diffusion coefficients D(J1 ) and D(J2 ) (with their cor- 0.155
responding error-bars) in the Veritas family region, shown here as functions 1.2e-014
of the proper semi-major axis ap . The embedded rectangles are magnifica-
tions of the graph, in the vicinity of the (3, 3, -2) and (7, -7, 2) MMRs. Note 1e-014
0.15
that D(J2 ) is practically zero everywhere, except in the (5,-2,-2) region.
8e-015
0.055 0.06 0.065 0.07
and D(J2 ) = (0.22±0.03)×10 −16 −1
yr , clearly much smaller than
in the (5,-2,-2) MMR. Note that D(J2 ) has practically zero value. ep
The region of the (7,-7,-2) resonance is even less exciting. D(J1 )
is almost constant across the resonance, with a very small value
(4.1±0.06 × 10−16 yr−1 ), and D(J2 ) is practically zero. Figure 6. The local diffusion coefficients D(J1 ) (top) and D(J2 ) (bottom)
The above results suggest that the (5,-2,-2) MMR is essentially inside the (5,-2,-2) region, shown here in grey-scale as functions of ep and
the only resonance in the Veritas region characterized by apprecia- sin Ip . The values are averaged over the resonant range in semi-major axis.
ble macroscopic diffusion. We now focus on the variation of the dif-
fusion coefficients with respect to (ep , sin Ip ) along this resonance.
A set of six values (ap ,J1 ,J2 ,P ,γ,H) is assigned to each random
As shown in Fig. 6, the dependence of the diffusion rate on the ini-
walker in the simulation. All bodies are initially distributed uni-
tial values of the actions4 is very complex. The values of D(J1 )
formly inside a region of predefined size in δJi (0) and semi-major
vary from (0.60±0.01)×10−14 yr−1 to (1.66±0.02)×10−14 yr−1
axes in the range [3.172, 3.176]. The age of the family, (τ ), is de-
while, for D(J2 ), they vary from (0.63±0.02)×10−14 yr−1 to
fined as the time needed for 0.3% of the random walkers to leave
(2.31±0.03)×10−14 yr−1 . Consequently, chaotic diffusion along
an ellipse in the (J1 , J2 ) plane, corresponding to a 3-σ confidence
this resonance can in principle produce asymmetric “tails” in the
interval of a 2-D Gaussian distribution. We note that, for Veritas,
distribution of group-A members. Note, however, that the coeffi-
the mobility in semi-major axis due to Yarkovsky is very small and
cients only vary by a factor of 2 − 3 and that their average values
practically insignificant for what concerns the estimation of its age,
are essentially the same as in Tsiganis et al. (2007).
since the family is young and distant from the Sun. For older fami-
lies one should also define appropriate borders in ap .
3.1.2 The MCMC Simulations - Chronology of Veritas Using the values of the coefficients obtained above, and the
values of σ(J1 )=(2.31±0.22)×10−4 , σ(J2 )=(3.97±0.30)×10−4 ,
We can now use the MCMC method to simulate the evolution and
calculated from the distribution of the real group A members, we
determine the age of the Veritas family, assuming that all the dy-
simulate the spreading of group A and estimate its age. Of course,
namically distinct groups originated from a single brake-up event.
the model depends on some free parameters: the initial spread of
the group in (δJ1 (0), δJ2 (0)), the time-step, dt, and the number
4 The values of the diffusion coefficients are calculated for a regular grid of random-walkers, n. Therefore, the dependence of the age, τ , on
in J1 and J2 and then translated into proper elements space. these parameters was checked. Uncertainties in the values of the

c 0000 RAS, MNRAS 000, 000–000
 Chaotic Transport and Chronology of Complex Asteroid Families 7
current borders of the group (i.e. the confidence ellipse) and the
values of the diffusion coefficients were taken into account, when -4 -4
n=2000, δJ1(0)=1.25 x 10 , δJ2(0)=5.6 x 10
calculating the formal error in τ .
11
Different sets of simulations were performed, the results of 10
which are given in Fig. 7(a)-(d). Each “simulation” (i.e. each point 9

τ [Myr]
8
in a plot) actually consists of 100 different realizations (runs) of
7
the MCMC code. In each run, the values of D(Ji ) and σ(Ji ) were 6
varied, according to the previously computed distributions of their 5
0 1000 2000 3000 4000 5000 6000
values. The values of the free parameters were the same for all runs
dt [yr]
in a given simulation.
The first set of simulations was performed in order to check
how the results depend on the time step, dt. Five simulations were
made, with dt ranging from 1000 to 5000 yr (Fig. 7a). The stan- -4 -4
dt=2000, δJ1(0)=1.25 x 10 , δJ2(0)=5.6 x 10
dard deviation of τ is relatively small, suggesting that τ is roughly
11
independent of dt. According to this set of simulations, the age of 10
the family is τ = hτ i ± ∆τ = 8.5 ± 1.3 Myr, where hτ i is the 9

τ [Myr]
8
mean value and ∆τ the standard error of the mean. This estimate
7
is in excellent agreement with those of Nesvorný et al. (2003) and 6
Tsiganis et al. (2007). 5
0 1000 2000 3000 4000 5000 6000
The second set of simulations was performed in order to check n
how the results depend on the number of random walkers, n. As
shown in Fig. 7b), τ is weakly dependent of n, the variation of the
mean value of τ is slightly larger than in the previous case. The age
of the family, according to this set, is τ = 8.6±1.3 Myr. n=2000, dt=2000 yr, δJ1(0)=2.30 x 10
-4

The third group of simulations was performed in order to 11

check how the results depend on the assumed initial spread of group 10
9
A in δJ2 (0), a parameter which is poorly constrained from the re-
τ [Myr]

8
spective equivelocity ellipse in Fig. 1. We fixed the value of δJ1 (0) 7
to 2.3×10−4 , which can be considered an upper limit, according 6

to Fig. (1)5 . Six simulations were performed, with δJ2 (0) ranging
5
2 4 6 8 10 12
from 3.5 ×10−4 to 11.0×10−4 , and the results are shown in Fig. δJ2(0) x 10
4

7c. The values of τ tend to decrease as δJ2 (0) increases. This is to

be expected, since increasing the initial spread of the family, while
targeting for the same final spread, should take a shorter time for a
given diffusion rate. The results yield τ = 8.8±1.1 Myr. -4
n=2000, δJ1(0)=2.3 x 10 , δJ2(0)=11.3 x 10
-4

As a final check, we performed a set of simulations with

11
δJ1 (0) = 2.3×10−4 and δJ2 (0) = 11.0×10−4 . These values cor- 10
respond to equivelocity ellipses that contain almost all regular and 9
τ [Myr]

8
(3,3,-2)-resonant family members, except for the very low inclina- 7
tion bodies (sin Ip < 0.156). Five sets of runs, for five different 6

values of dt, were performed (see Fig. 7d). As in our first set of 5
0 1000 2000 3000 4000 5000 6000
simulations, τ is practically independent of dt. On the other hand, dt [yr]
τ turns out to be smaller than in the previous simulations, since the
assumed values for δJ1 (0) and δJ2 (0) are quite large. Even so, we
find τ = 7.6±1.1 Myr, which is still an acceptable value.
Figure 7. Dependence of τ and σ(τ ) on three free parameters: (a) time step
Combining the results of the first three sets of simulations and
dt; (b) number of random walkers n; (c) initial size of the family δJ2 (0).
taking into account all uncertainties, we find an age estimate of τ =
The bottom panel, 7(d), shows the dependence of τ on dt for a run in which
8.7±1.2 Myr for the Veritas family. This result is very close to the the maximum values of δJ1 (0) and δJ2 (0) were assumed. The values of
one found by Tsiganis et al. (2007), the error though being smaller the free parameters that are constant in each group of simulations are indi-
by ∼ 30%. cated on top of each plot. The horizontal lines and dashed areas denote the
In addition to the determination of the family’s age, we would age of the Veritas family (and the corresponding errorbar), as obtained by
like to know how well the MCMC model reproduces the evolu- Nesvorný et al. (2003).
tion of the spread of group-A bodies, in the (J1 ,J2 ) space. For this
purpose we compared the evolution of group A for 10 Myr in the
future, as given (i) by direct numerical integration of the orbits, and Fig. 8(a), the random-walkers of the MCMC simulation (triangles)
(ii) by an MCMC simulation with variable diffusion coefficients. practically cover the same region in (J1 ,J2 ) as the real group-A
Figures 8(a)-(b) show the outcome of this comparison. As shown in members (circles). Moreover, the time evolution of the ratio of the
standard deviations σ(J2 )/σ(J1 ), which characterizes the shape of
the distribution, is reproduced quite well, as shown in Fig. 8(b).
5 Assuming that an equivelocity ellipse in (ap ,ep ), large enough to encom- An additional MCMC simulation with constant (i.e. average
pass both the regular part of the family and the (3,3,-2) bodies, is a better with respect to ep and sin Ip ) coefficients of diffusion was also per-
constraint. formed. As shown in Fig. 8(b), the value of σ(J2 )/σ(J1 ) in this

c 0000 RAS, MNRAS 000, 000–000
8 B. Novaković et al.

0.0115 0.215
+10 Myr
0.011 0.21

0.205
0.0105
0.2

ep
J2

0.01
0.195
0.0095
0.19

0.009 (5900) Jensen

0.185 (3330) Gantrisch

0.0085 0.18
0 0.0005 0.001 0.0015 0.002 0.0025 0.003 3.13 3.135 3.14 3.145 3.15 3.155 3.16 3.165 3.17
J1 ap [AU]

0.184
(3330) Gantrisch
0.182

0.18

0.178
sin Ip

0.176

0.174

0.172 (5900) Jensen

0.17

3.13 3.135 3.14 3.145 3.15 3.155 3.16 3.165 3.17

ap [AU]

16

15

14
H [mag]

Figure 8. Comparison of the time evolution of group A Veritas members 13

in (i) a numerical integration of their orbits for 10 Myr, and (ii) an MCMC
realization, corresponding to the same time interval. Top: the final distribu- 12
tion, as taken from the integration (circles) and the MCMC run (triangles). (5900) Jensen
Bottom: time evolution of the ratio σ(J2 )/σ(J1 ). The thick solid (resp. 11
(3330) Gantrisch
dashed) curve corresponds to the MCMC run with variable (resp. constant)
diffusion coefficients, while the thin, “noisy” one corresponds to the numer- 3.13 3.135 3.14 3.145 3.15 3.155 3.16 3.165 3.17
ical integration of the equations of motion. ap [AU]

Figure 9. Top and middle panels: The same as Fig. 1, but for the Lixi-
simulation appears to slowly deviate from the one measured in the
aohua family. The different symbols correspond to different values of the
previous MCMC simulation, as time progresses. However, this de-
Lyapunov time, TL . Triangles correspond to TL > 105 y, empty circles to
viation is not very large, also compared to the result of the numeri- 2 × 104 ≤ TL ≤ 105 y and solid circles to TL < 2 × 104 y. The equive-
cal integration. Thus, we conclude that an MCMC model with con- locity curves were obtained for vc = 40 m s−1 , f = 80◦ and ω = 300◦ .
stant coefficients is adequate for deriving a reasonably accurate es- Bottom: The distribution of Lixiaohua members in the (ap , H) plane. The
timate of the age of a family, provided that the variations of the local grey-shaded region denotes the extent of the main chaotic zone (MCZ) in
diffusion coefficients are not very large. Given this result, we de- ap . The two largest members of the family are also denoted.
cided to use average coefficients for the Lixiaohua case. Note also
that, in the Veritas case, the observed deviation in σ(J2 )/σ(J1 )
between the two MCMC models is reflected in the error of τ (i.e. ily members that follow chaotic trajectories. At the same time, a
1.7 Myr vs. 1.2 Myr). clear ‘V’-shaped distribution is observed in the (ap , H) plane (see
Fig. 9), suggesting that the family is old-enough for Yarkovsky to
have significantly altered its size in ap . Nesvorný et al. (2005) in
3.2 The age of the Lixiaohua family
this way estimated the age of this family to τ = 300 ± 200 Myr.
The Lixiaohua family is another typical outer-belt family, crossed Here, we use our MCMC method to derive a more accurate estimate
by several MMRs. This results into a significant component of fam- of its age, taking into account also the Yarkovsky/YORP effects.

c 0000 RAS, MNRAS 000, 000–000
 Chaotic Transport and Chronology of Complex Asteroid Families 9
The distribution of the family members in (ap , ep ) and 9e-15
(ap , sin Ip ) is shown in Fig. 9. Using vc = 50 m s−1 we find
8e-15
263 bodies linked to the family. The shape of this family is, as
in the Veritas case, intriguing. For a > 3.15 AU, the family ap- 7e-15
pears to better fit inside the equivelocity ellipse shown in the fig- 6e-15

D(J1) [yr-1]
ure, with only a few bodies showing a significant excursion in ep 5e-15
and sin Ip . On the other hand, for a < 3.145 AU, the family mem-
4e-15
bers occupy a wider area in ep and sin Ip . Throughout the family
region we find thin, “vertical”, strips of chaotic bodies, with Lya- 3e-15
punov times6 TL < 2 × 104 y. These strips are associated to differ- 2e-15
ent mean motion resonances. The most important chaotic domain 1e-15
(hereafter Main Chaotic Zone, MCZ) is the one centered around
0
ap ≈ 3.146 AU; indicated by the grey-shaded area in both plots of 3.13 3.135 3.14 3.145 3.15 3.155 3.16 3.165 3.17
Fig. 9. A number of two- and three-body MMRs can be associated ap [AU]
to the formation of the MCZ, such as the 17:8 MMR with Jupiter
and the (7, 9, -5) three-body MMR. 4e-15
Note that the two largest members of this family, (3330)
3.5e-15
Gantrisch and (5900) Jensen, are located just outside the MCZ, as
indicated by their larger values of TL (> 2 × 104 y). In fact, a 3e-15
significant group of bodies just outside the MCZ (see Fig. 9a) has
2.5e-15
D(J2) [yr-1]
higher values of TL but similar spread in ep as the MCZ bodies.
This suggests that bodies around the chaotic zone could have once 2e-15
resided therein, evolving towards high/low values of ep by chaotic 1.5e-15
diffusion. Numerical integrations of the orbits of selected Lixiao-
hua members for 100 Myr indeed confirmed that bodies could enter 1e-15
(or leave) the MCZ. We believe that the distribution of family mem- 5e-16
bers on either side of the MCZ is strongly indicative of an interplay
0
between Yarkovsky dirft in semi-major axis and chaotic diffusion 3.13 3.135 3.14 3.145 3.15 3.155 3.16 3.165 3.17
in ep and sin Ip , induced by the overlapping resonances; bodies ap [AU]
can be forced to cross the MCZ, thus receiving a “kick” in ep and
sin Ip , before exiting on the other side of the zone.
A population of ≈ 5, 000 fictitious bodies was selected and Figure 10. The same as Fig. 5, but for the Lixiaohua family. Note that while
several diffusive peaks are visible in D(J1 ), only the region of the MCZ
used for calculation of the local diffusion coefficients. Here, we re-
(ap ∼ 3.146 AU) shows significant diffusion in J2 .
strict ourselves in calculating coefficients as functions of ap only
(i.e. averaged in ep and sin Ip ). As shown in Fig. 10, there are sev-
eral diffusive zones, corresponding to the low-TL strips of Fig. 9. the Yarkovsky effect. Given the uncertainty in determining the ini-
However, as in the Veritas case, only one zone appears diffusive tial size of the family, we repeated the computations for two more
in both actions; D(J2 ) is practically zero everywhere outside the sets of δJ1 (0), δJ2 (0). The results of these computations are pre-
MCZ (a ∼ 3.146 AU), and significant dispersion in proper ele- sented in Fig. 11. We find an upper limit of ∼ 230 Myr for the age
ments is observed only in this zone. of the family and a lower limit of ∼ 100 Myr. Taking into account
Given the above results, we conclude that the MCZ family only the runs performed for our “nominal” initial size and includ-
members can be used to estimate the age of the family, much like ing the Yarkovsky effect, we find the age of Lixiaohua family to be
the Veritas (5,-2,-2) resonant bodies. Given the fact that random- 155 ± 36 Myr. This value lies towards the lower end but comfort-
walkers can drift in ap , the way of computing the age is accordingly ably within the range (300 ± 200 Myr) given by Nesvorný et al.
modified. A large number of random walkers, uniformly distributed (2005).
across the whole family region (i.e. the equivelocity ellipses) is We note that, when the Yarkovsky effect is taken into account,
used. The simulation again stops when 0.3% of MCZ-bodies are the age of the family turns out to be longer by ∼ 30 Myr (see also
found to be outside the observed (J1 , J2 ) borders of the family. Fig. 11). This is a purely dynamical effect, related to the fact that
However, the number of MCZ bodies is not constant during the bodies can drift towards the MCZ from the adjacent non-diffusive
simulation, because bodies initially outside (resp. inside) the MCZ regions. However, as more bodies enter the MCZ near its center
can enter (resp. leave) that region. Thus, the aforementioned per- in ep and sin Ip , it takes longer for 0.3% of random walkers to
centage is calculated with respect to the corresponding number of diffuse outside the confidence ellipse in (ep , sin Ip ). At the same
MCZ bodies at each time-step. time, bodies that are initially inside the MCZ can also drift outside,
The size of the MCZ in the space of proper actions is given by to lower/higher values of ap , thus slowing down or even stop diffus-
σ(J1 ) = (8.49±0.51)×10−4 and σ(J2 ) = (3.17±0.28)×10−4 . ing in ep and Ip . This also explains the large spread in (ep , sin Ip )
In order to compute the age of the family we performed 2400 observed for family members located just outside the MCZ.
MCMC runs. This was repeated three times, for three different val-
ues of thermal conductivity (see above) and once more, neglecting
4 CONCLUSIONS
6 For details on the computation of Lyapunov times see e.g. Tsiganis et al. We have presented here a refined statistical model for asteroid
(2003) transport, which accounts for the local structure of the phase-space,

c 0000 RAS, MNRAS 000, 000–000
10 B. Novaković et al.
that its age is known by independent means (e.g. by applying the
method of Nesvorný et al. (2003) to the regular members of the
220
family). A large number of MCMC runs can be performed at low
200 computational cost, thus allowing a thorough analysis of the phys-
ical parameters of family members or the properties of the original
180 ejection velocities field that better reproduce the currently observed
τ [Myr]

shape of the family.

160

140
ACKNOWLEDGEMENTS
120
K.T. would like to thank Z.K. and B.N. for their warm hospitality
100 during his stay in Belgrade. The work of B.N. and Z.K. has been
3e-08 5e-08 7e-08 9e-08 1.1e-07 supported by the Ministry of Science and Technological Develop-
Size of the initial ellipse ment of the Republic of Serbia (Project No 146004 ”Dynamics of
Celestial Bodies, Systems and Populations”).
Figure 11. The age of Lixiaohua family, as measured by our MCMC runs.
The mean value and standard deviation (error-bar) of the age are given, as
functions of the assumed size of the initial equivelocity ellipse. Four sets of
points are shown, corresponding to our estimate (i) including Yarkovsky and REFERENCES
assuming different values of the thermal conductivity (K = 0.01 = dotted
line, K = 0.1 = dashed line and K = 0.5 = solid thin line), and (ii) Beaugé, C., Roig, F., 2001, Icarus, 153, 391
neglecting Yarkovsky (thick solid line). Bendjoya, Ph., Zappalà, V., 2002, in: W. F. Bottke Jr., A. Cellino,
P. Paolicchi, and R. P. Binzel (eds.), Asteroids III, University of
Arizona Press, Tucson, 613
by using variable diffusion coefficients. Also, the model takes into Berg, A.B., 2004, Markov Chain Monte Carlo Simulations and
account the long-term drift in semi-major axis of asteroids, induced Their Statistical Analysis, World Scientific, Singapore
by the Yarkovsky/YORP effects. This model can be applied to sim- Bottke, W.F., Durda, D.D., Nesvorný, D., Jedicke, R., Morbidelli,
ulate the evolution of asteroid families, also giving rise to an ad- A., Vokrouhlickỳ, D., Levison, H.F., 2005, Icarus 179, 63
vanced version of the ”chaotic chronology” method for the deter- Brož, M., Vokrouhlický, D., Roig, F., Nesvorný, D., Bottke, W.F.,
mination of the age of asteroid families. Morbidelli, A., 2005, MNRAS, 359, 1437
We applied our model to the Veritas family, whose age is well Brož, M., 2006, Yarkovsky effect and the dynamics of Solar sys-
constrained from previous works. This allowed us to assess the tem. Ph.D. Thesis, Faculty of Mathematics and Physics, Charles
quality and to calibrate our model. We first analyzed the local dif- University, Prague.
fusion characteristics in the region of Veritas. Our results showed Brown, M.E., Barkume, K.M., Ragozzine, D., Schaller, E.L.,
that local diffusion coefficients vary by about a factor of 2 − 3 2007, Nature, 446, 294
across the (ep , sin Ip ) region covered by the (5,-2,-2) MMR. Thus, Čapek, D., Vokrouhlický, D., 2004, Icarus, 172, 526
although local coefficients are needed to accurately model (by the Carruba, V., Burns, J.A., Bottke, W.F., Nesvorný, D., 2003, Icarus
MCMC method) the evolution of the distribution of group-A mem- 162, 308
bers, average coefficients are enough for a reasonably accurate es- Cellino, A., Bus, S.J., Doressoundiram, A., Lazzaro, D., 2002, in:
timation of the family’s age. We note though that the variable co- W. F. Bottke Jr., A. Cellino, P. Paolicchi, and R. P. Binzel (eds.),
efficients model reduces the error in τ by ∼ 30%, but requires Asteroids III, University of Arizona Press, Tucson, 633
a computationally expensive procedure. Using the variable coeffi- Di Martino, M., Migliorini, F., Zappalà, V., Manara, A., Barbieri,
cients MCMC model, we found the age of the Veritas family to be C., 1997, Icarus, 127, 112
τ = (8.7±1.2) Myr; a result in very good agrement with that of Donnison, J.R., 2003, MNRAS, 338, 452
Nesvorný et al. (2003) and Tsiganis et al. (2007). Durda, D.D., Bottke, W.F., Nesvorný, D., Enke, B.L., Merline,
We used our model to estimate also the age of the Lixiao- W.J., Asphaug, E., Richardson, D.C., 2007, Icarus, 186, 498
hua family. This family is similar to the Veritas family in many re- Farinella, P., Vokrouhlickỳ, D., 1999, Science, 283, 1507
spects; it is a typical outer-belt family of C-type asteroids, crossed Farley, K.A., Vokrouhlický, D., Bottke, W.F., Nesvorný, D., 2006,
by several MMRs. Like the Veritas family, only the main chaotic Nature, 439, 295
zone (MCZ) shows appreciable diffusion in both eccentricity and Gentle, E.J., 2003, Random Number Generation and Monte Carlo
inclination. On the other hand, this is a much older family and the Methods. Springer, New York.
Yarkovsky effect can no longer be ignored. This is evident from the Hirayama, K., 1918, Astron. J., 31, 185
distribution of family members, adjacent to the MCZ. Our model Knežević, Z., Milani, A., 2000, Cel. Mech. Dyn. Astron., 78, 17
suggests that the age of this family is between 100 and 230 Myr, Knežević, Z., Pavlović, R., 2002, Earth Moon Planets 88, 155
the best estimate being 155 ± 36 Myr. Note that the relative error Kryszczyńska, A., La Spina, A., Paolicchi, P., Harris, A. W., Bre-
is ≈ 23%, i.e. close to the ≈ 20% that Tsiganis et al. (2007) found iter, S., Pravec, P., 2007, Icarus, 192, 223
for the Veritas case, using a constant coefficients MCMC model. Lazzaro, D., Angeli, C.A., Carvano, J.M., Mothé-Diniz, T., Duf-
The transport model that we presented here can be used to fard, R., Florczak, M., 2004, Icarus, 172, 179
study other families that reside in complex phase-space regions, Lichtenberg, A.J., Lieberman, M.A. 1983. Regular and Stochastic
such as those of (20) Massalia and (778) Theobalda. An important Motion. Springer- Verlag New York, N.Y., U.S.A.
advantage of the model is that it can be used to estimate the phys- Michel, P., Benz, W., Richardson, D.C., 2003, Nature, 421, 608
ical properties of a dynamically complex asteroid family, provided Milani, A., 1993, Cel. Mech. Dyn. Astron., 57, 59

c 0000 RAS, MNRAS 000, 000–000
 Chaotic Transport and Chronology of Complex Asteroid Families 11
Milani, A., Farinella, P., 1994, Nature 370, 40
Milani, A., Knežević, Z., 1990, Cel. Mech. Dyn. Astron., 98, 211
Milani, A., Knežević, Z., 1994, Icarus, 107, 219
Mothé-Diniz, T., Roig, F., Carvano, J.M., 2005, Icarus, 174, 54
Morbidelli, A., Zappalà, V., Moons,M., Cellino, A., Gonczi, R.,
1995, Icarus 118, 132
Morbidelli, A., 2002, Modern Celestial Mechanics: Aspects of
Solar System Dynamics. Taylor & Francis, London.
Murray, N., Holman, M., 1997, Astron. J., 114, 1246
Nesvorný, D., Morbidelli, A., 1999, Celest. Mech. Dynam. As-
tron., 71, 243
Nesvorný, D., Bottke, W.F., Dones, L., Levison, H.F., 2002, Na-
ture, 417, 720
Nesvorný, D., Bottke, W.F., Levison, H.F., Dones, L., 2003, As-
trophys. J., 591, 486
Nesvorný, D., Jedicke, R., Whiteley, R.J., Ivezić, Ž., 2005, Icarus
173, 132
Nesvorný, D., Bottke, W.F., Vokrouhlický, D., Morbidelli, A.,
Jedicke, R., 2006, in: Daniela, L., Sylvio Ferraz, M., and An-
gel, F. Julio (eds.) Proceedings of the 229th Symposium of the
IAU, Asteroids, Comets, Meteors, Cambridge University Press,
Cambridge, 289
Paolicchi, P., Verlicchi, A., Cellino, A., 1996, Icarus, 121, 126
Tedesco, E.F., Noah, P.V., Noah, M., Price, S.D., 2002, Astron. J.,
123, 1056
Tsiganis, K., Varvoglis, H., Morbidelli, A., 2003, Icarus 166, 131
Tsiganis, K., Knežević, Z., Varvoglis, H., 2007, Icarus, 186, 484
Vokrouhlický, D., 1999, Astron. Astrophys., 344, 362
Vokrouhlický, D., Čapek, D., 2002, Icarus, 159, 449
Vokrouhlický, D., Brož, M., Morbidelli, A., Bottke, W.F.,
Nesvorný, D., Lazzaro, D., Rivkin, A.S., 2006, Icarus, 182, 92
Vokrouhlický, D., Brož, M., Bottke, W.F., Nesvorný, D., Mor-
bidelli, A., 2006, Icarus, 182, 118
Zappalà, V., Bendjoya, Ph., Cellino, A., Farinella, P., Froeschlé,
C., 1995, Icarus, 116, 291

c 0000 RAS, MNRAS 000, 000–000