# Did the Hilda collisional family form during the late heavy bombardment?

M. Bro, D. Vokrouhlický, A. Morbidelli, D. Nesvorný, W. F. Bottke

► To cite this version:

M. Bro, D. Vokrouhlický, A. Morbidelli, D. Nesvorný, W. F. Bottke. Did the Hilda collisional family form during the late heavy bombardment?. Monthly Notices of the Royal Astronomical Society, 2011, 414, pp.596. 10.1111/j.1365-2966.2011.18587.x. hal-00632637

HAL Id: hal-00632637

https://hal.science/hal-00632637v1

Submitted on 23 Sep 2021

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons CC BY 4.0 - Attribution - International License
# Did the Hilda collisional family form during the late heavy bombardment?

M. Brož,<fnref n="1" />* D. Vokrouhlický,<fnref n="1" /> A. Morbidelli,<fnref n="2" /> D. Nesvorný<fnref n="3" /> and W. F. Bottke<fnref n="3" />

Accepted 2011 February 22. Received 2011 February 3; in original form 2011 January 5

## ABSTRACT

We model the long-term evolution of the Hilda collisional family located in the 3/2 mean-motion resonance with Jupiter. Its eccentricity distribution evolves mostly due to the Yarkovsky/YORP effect and assuming that (i) impact disruption was isotropic and (ii) albedo distribution of small asteroids is the same as for large ones, we can estimate the age of the Hilda family to be $4_{-1}^{+0}$ Gyr. We also calculate collisional activity in the J3/2 region. Our results indicate that current collisional rates are very low for a 200-km parent body such that the number of expected events over gigayears is much smaller than 1.

The large age and the low probability of the collisional disruption lead us to the conclusion that the Hilda family might have been created during the late heavy bombardment (LHB) when the collisions were much more frequent. The Hilda family may thus serve as a test of orbital behaviour of planets during the LHB. We have tested the influence of the giant-planet migration on the distribution of the family members. The scenarios that are consistent with the observed Hilda family are those with fast migration time-scales $\sim$0.3–3 Myr, because longer time-scales produce a family that is depleted and too much spread in eccentricity. Moreover, there is an indication that Jupiter and Saturn were no longer in a compact configuration (with period ratio $P_S/P_J > 2.09$) at the time when the Hilda family was created.

**Key words:** methods: numerical – celestial mechanics – minor planets, asteroids: general.

## 1 INTRODUCTION

There are many independent lines of evidence that the orbits of planets of the Solar system were not the same all the time, but that they have changed substantially over billions of years. The major arguments are based on the observed orbital distribution of Kuiper belt objects (Malhotra 1995; Levison et al. 2008) or small but non-negligible eccentricities and inclinations of the giant planets (Tsiganis et al. 2005). Observations of Jupiter's Trojans (Morbidelli et al. 2005), main-belt asteroids (Minton & Malhotra 2009; Morbidelli et al. 2010), the amplitudes of secular oscillations of the planetary orbits (Brasser et al. 2009; Morbidelli et al. 2009), or the existence of irregular moons (Nesvorný, Vokrouhlický & Morbidelli 2007) provide important constraints for planetary migration scenarios.

Asteroids are a fundamental source of information about the evolution of the planetary system. Some of the resonant groups, i.e. those which are located in the major mean-motion resonances with Jupiter, might also have been influenced by planetary migration, because their current distribution does not match the map of the currently stable regions. For instance, there are two stable islands

denoted by A and B in the J2/1 resonance and only the B island is populated (Nesvorný & Ferraz-Mello 1997).

In this work we focus on the Hilda asteroid family in the 3/2 resonance with Jupiter. We exploit our ability to model long-term evolution of asteroid families, which is usually dominated by the Yarkovsky effect on the orbital elements (Bottke et al. 2001), often coupled to the YORP effect on the spin rate and obliquity (Vokrouhlický et al. 2006b). Chaotic diffusion in eccentricity and sometimes interactions with weak mean-motion or secular resonances (Vokrouhlický et al. 2006a) also play important roles. In case of asteroids inside strong mean-motion resonances, one has to account for the 'resonant' *Yarkovsky effect*, which causes a systematic drift in eccentricity (Brož & Vokrouhlický 2008). This is different from usual non-resonant orbits where the Yarkovsky effect causes a drift in semimajor axis.

The Hilda collisional family – a part of the so-called Hilda group in the 3/2 mean motion resonance with Jupiter – was already briefly discussed by Brož & Vokrouhlický (2008). However, the modelling presented in that paper was not very successful, since the resulting age of the family seemed to be too large (exceeding 4 Gyr). This was an important motivation for our current work. We think that we missed an important mechanism in our previous model, namely perturbations arising from the migration of the giant planets and also an appropriate treatment of the YORP effect. Indeed, the age

<footnote n="1">¹Institute of Astronomy, Charles University, Prague, V Holešovičkách 2, 18000 Prague 8, Czech Republic</footnote>

<footnote n="2">²Observatoire de la Côte d'Azur, BP 4229, 06304 Nice Cedex 4, France</footnote>

<footnote n="3">³Department of Space Studies, Southwest Research Institute, 1050 Walnut St., Suite 300, Boulder, CO 80302, USA</footnote>

<footnote>*E-mail: mira@sirrah.troja.mff.cuni.cz</footnote>
≥4 Gyr suggests that the planetary migration might have played a direct role during the early evolution of the Hilda family. In this paper we thoroughly test this hypothesis.

The paper is organized as follows. First we study the observed properties of the J3/2 resonance population in Section 2. Our dynamical model of the Hilda family (without migration first) is described in Section 3. Then we estimate the collisional activity in the J3/2 region in Section 4. The results of our simulations of the giant-planet migration are presented in Section 5. Finally, Section 6 is devoted to conclusions.

## 2 CURRENT ASTEROID POPULATION IN THE J3/2 RESONANCE

Asteroids located in the 3/2 mean motion resonance with Jupiter have osculating semimajor axes around (3.96 ± 0.04) au, i.e. beyond the main asteroid belt. Contrary to the Kirkwood gaps (associated with J3/1, J7/3 or J2/1 resonances), this resonance is populated by asteroids while its neighbourhood is almost empty. The Hilda collisional family we are going to discuss in detail is a small part of the whole J3/2 resonant population.

Our identification procedure of the J3/2 resonant population was described in the previous paper, Brož & Vokrouhlický (2008). Using the AstOrb catalogue of orbits (version JD = 245 5500.5, 2010 October 31) we identified 1787 numbered and multi-opposition bodies with the librating critical argument

$$ \sigma = \frac{p+q}{q} \lambda' - \frac{p}{q} \lambda - \omega, \quad (1) $$

where $p = 2$, $q = 1$, $\lambda'$ is the mean longitude of Jupiter, $\lambda$ the mean longitude of the asteroid and $\omega$ the longitude of perihelion of the asteroid.

In order to study the detailed distribution of the bodies librating inside the resonance, we have to use pseudo-proper resonant elements defined as approximate surfaces of sections (Roig, Nesvorný & Ferraz-Mello 2002), i.e. the intersection of the trajectory with a plane defined by

$$ |\sigma| < 5^{\circ}, \quad \frac{\Delta\sigma}{\Delta t} > 0, \quad |\omega - \omega'| < 5^{\circ}. \quad (2) $$

These conditions correspond to the maximum of the semimajor axis $a$ over several oscillations and the minimum of the eccentricity $e$ or the inclination $I$. We need to apply a digital filter to $\sigma(t)$ prior to using equation (2), namely filter A from Quinn, Tremaine & Duncan (1991), by sampling 1 yr and with a decimation factor of 10, to suppress fast ≈80 yr oscillations, which would otherwise disturb slower ≈280 yr oscillations associated with resonant libra- tions. Finally, we apply an averaging of the sections $a, e, I$ over 1-Myr running window and these averages are the pseudo-proper elements $a_p, e_p, I_p$. The accuracy of the pseudo-proper elements is of the order of $10^{-4}$ au for $a_p$ and $10^{-4}$ for $e_p$ or $\sin I_p$, which is much smaller than those of the structures we are interested in.

The overall dynamical structure of the J3/2 resonance is deter- mined by secular resonances $v_5, v_6$ at high eccentricities $e_p \gtrsim 0.3$ and secondary resonances at lower values of $e_p \lesssim 0.13$ (according to Morbidelli & Moons 1993; Nesvorný & Ferraz-Mello 1997; Ferraz- Mello et al. 1998; Roig & Ferraz-Mello 1999). They destabilize the orbits at the borders of a stable island. The orbits inside the island exhibit very low chaotic diffusion rates, so bodies can remain there for 4 Gyr (without non-gravitational perturbation).

Next we apply a hierarchical clustering method (Zappalà et al. 1994) to detect significant clusters. We use a standard metric in the

<figure><img src="image_2.png" /><figcaption>Figure 1. The number $N$ of the Hilda family members versus the selected cut-off velocity $v_{\text{cut-off}}$.</figcaption></figure>

pseudo-proper element space ($a_p, e_p, \sin I_p$):

$$ \delta v = na\sqrt{\frac{5}{4}\left(\frac{\delta a_p}{a_p}\right)^2 + 2(\delta e_p^2) + 2(\delta \sin I_p)^2}. \quad (3) $$

In the following, we do not discuss the known Schubart family, which was sufficiently analysed elsewhere (Brož & Vokrouhlický 2008), but we focus on the family associated with (153) Hilda. A suitable cut-off velocity for the Hilda family seems to be $v_{\text{cut-off}} = 140 \text{ m s}^{-1}$, because the number of members does not change substantially around this value (see Fig. 1). The number of members at this cut-off is 400.

The resulting plots ($a_p, H$), ($e_p, H$) and ($I_p, H$) of the Hilda family show very interesting features (see Fig. 2). The distribution of semimajor axis and inclination seems rather uniform and almost independent of absolute magnitude $H$, but the eccentricities of small asteroids (i.e. with high $H$) are clearly concentrated at the outskirts of the family and depleted in the centre.

In order to explain the distribution of asteroids in the ($e_p, H$) plane we have to recall that asteroids orbiting about the Sun are affected by non-gravitational forces, mostly by the Yarkovsky/YORP effect, i.e. the recoil force/torque due to anisotropic emission of thermal radiation. We consider the concentrations in the ($e_p, H$) plane to be a strong indication of the ongoing Yarkovsky/YORP evolution, because they are very similar to those observed among the several main-belt families in the ($a_p, H$) plane and successfully modelled by Vokrouhlický et al. (2006b). The difference between these two cases stems from the fact that the main-belt families are non-resonant and the Yarkovsky/YORP effect thus increases or decreases the semimajor axis (depending on the actual obliquity of the spin axis), while in our resonant case, the same perturbation results instead in a systematic increase or decrease of eccentricity. A detailed modelling of the $e$-distribution is postponed to Section 3.5.

The central part of the ($e_p, H$) distribution, from $e = 0.17–0.23$, seems rather extended. The large asteroids ($H < 12.5$ mag) are spread over this interval of eccentricities even though their Yarkovsky drift rates must have been small. Only 2–4 of them are likely to be interlopers, because there is a very small number of background asteroids in the surroundings of the family (see Fig. 3). We think this shape might actually be the result of the initial-size- independent perturbation that the family distribution received by the migration of the giant planets (which we discuss in Section 5.1).

Regarding the ($a_p, H$) distribution, the largest asteroid (153) Hilda is offset with respect to the centre, but this is a natural outcome of the definition of the pseudo-proper elements – fragments that fall to the left of the libration centre are mapped to the right, which creates the offset.

The geometric albedos for Hilda family objects are poorly known. There are only six measured values for the family members: 0.064,
<figure><img src="image_3.png" /><figcaption>**Figure 2.** The Hilda family displayed in resonant semimajor axis $a_p$ (left), eccentricity $e_p$ (middle) and inclination $\sin I_p$ (right) versus absolute magnitude $H$. The libration centre is located at $a \simeq 3.96$ au and all the bodies are displayed to the right of it. The ‘ears’ in $(e_p, H)$, i.e. both the concentration of small asteroids on the outskirts of the family and their depletion at the centre are very prominent here. The thin vertical lines denote the central part of the $(e_p, H)$ distribution discussed in the text. The family has 400 members at $v_{\text{cut-off}} = 140$ m s$^{-1}$.</figcaption></figure>

<figure><img src="image_4.png" /><figcaption>**Figure 3.** The J3/2 region displayed in $(a_p, I_p)$ plot. A very prominent Schubart cluster (studied by Brož & Vokrouhlický 2008) is visible around $\sin I_p \simeq 0.05$. The close surroundings of the Hilda family, where only a low number of bodies is present, are highlighted by grey rectangles.</figcaption></figure>

0.046, 0.038, 0.089, 0.044, 0.051 (Davis & Neese 2002). Given the small number of values and the possibility of selection effects, we prefer to assume that the family members have a mean value $p_V = 0.044$, which corresponds to the whole J3/2 population. The size of the parent body can be then estimated to be $D_{\text{PB}} = (200 \pm 20)$ km. We employ two independent methods to determine the diameter $D_{\text{PB}}$: (i) we sum the volumes of the observed bodies larger than an assumed completeness limit $D_{\text{complete}} = 10$ km and then we prolong the slope of the size-frequency distribution down to $D = 0$ to account for unobservable bodies (see Brož & Vokrouhlický 2008), which results in $D_{\text{PB}} \simeq 185$ km; (ii) we also use a geometric method developed by Tanga et al. (1998) which gives $D_{\text{PB}} \simeq 210$ km. A test with different albedo values will be described in Section 3.6.

The size-frequency distribution $N(>D)$ versus $D$ of the Hilda family is steeper than that of background J3/2 population, but shallower than for usual main-belt families (Fig. 4). Interestingly, the slope $\gamma = -2.4 \pm 0.1$ of the distribution $N(>D) = CD^\gamma$ is close to a collisional equilibrium calculated by Dohnanyi (1969).

Colour data extracted from the Sloan Digital Sky Survey Moving Object Catalogue version 4 (Parker et al. 2008) confirm that the Hilda family belongs to the taxonomic type C, because most of the spectral slopes are small. Recall that the whole J3/2 population exhibits a bimodal distribution of slopes, i.e. it contains a mixture of C- and D-type asteroids.

<figure><img src="image_5.png" /><figcaption>**Figure 4.** Cumulative size distributions of the J3/2 population and the Hilda family. The polynomial fits of the form $N(>D) = CD^\gamma$ are plotted as thin lines, together with the respective values of the $\gamma$ exponent. Several main-belt families are plotted for comparison: Eos (with slope $\gamma = -2.8$), Euno-mia ($-5.0$), Hygiea ($-3.8$), Koronis ($-2.8$), Themis ($-2.9$), Tirela ($-3.3$), Veritas ($-3.4$) and Vesta ($-5.4$).</figcaption></figure>

## 3 THE HILDA FAMILY MODEL WITH RADIATION FORCES

To understand the long-term evolution of the Hilda family, we construct a detailed numerical model, extending efforts in Brož & Vokrouhlický (2008), which includes the following processes: (i) impact disruption, (ii) the Yarkovsky effect, (iii) the YORP effect, (iv) collisions and spin-axis reorientations. We describe the individual parts of the model in the following subsections.

### 3.1 Impact disruption

To obtain initial conditions for the family just after the breakup event, we need a model for the ejection velocities of the fragments. We use a very simple model of an isotropic ejection from the work of Farinella, Froeschlé & Gonczi (1994). The distribution of velocities ‘at infinity’ follows the function

$$dN(v)dv = C v(v^2 + v_{\text{esc}}^2)^{-(\alpha+1)/2}dv, \quad (4)$$

with the exponent $\alpha$ being a free parameter, C a normalization constant and $v_{\text{esc}}$ the escape velocity from the parent body, which
is determined by its size $D_{PB}$ and mean density $\rho_{PB}$ as $v_{esc} = \sqrt{(2/3)\pi G \rho_{PB}} D_{PB}$. The distribution is usually cut at a selected maximum allowed velocity $v_{max}$ to prevent outliers. The actual values of all these parameters are given in Section 3.5. Typically, the overall distribution of velocities has a peak close to the escape velocity, which is approximately $100 \text{ m s}^{-1}$ for a 200-km parent body. The initial velocities $|v|$ of individual bodies are generated by a straightforward Monte Carlo code and the orientations of the velocity vectors $v$ in space are assigned randomly.

Here, we assume the velocity of fragments is independent of their size, which seems reasonable with respect to the observed uniform distribution of the Hilda family in the ($a_p, H$) and ($I_p, H$) planes (Fig. 2). We also perform tests with non-isotropic distributions in Section 3.7.

We must also select initial osculating eccentricity $e_i$ of the parent body, initial inclination $i_i$, as well as true anomaly $f_{imp}$ and argument of perihelion $\omega_{imp}$ at the time of impact disruption. All of these parameters determine the initial shape of the synthetic ‘Hilda’ family just after the disruption of the parent body. Initial semimajor axis $a_i$ is not totally free, instead it is calculated from the initial semimajor axis of Jupiter $a_{Ji}$ and the Kepler’s law, since the parent body has to be confined to the J3/2 resonance.

## 3.2 Yarkovsky effect in a resonance

The long-term evolution of asteroid orbits is mainly driven by the Yarkovsky thermal effect. The implementation of the Yarkovsky effect in the SWIFT integrator was described in detail in Brož (2006). Only minor modifications of the code were necessary to incorporate spin rate evolution, which is driven by the YORP effect (see Section 3.3).

The thermal parameter we use are reasonable estimates for C/X-type bodies: $\rho_{surf} = \rho_{bulk} = 1300 \text{ kg m}^{-3}$ for the surface and bulk densities, $K = 0.01 \text{ W m}^{-1} \text{ K}^{-1}$ for the surface thermal conductivity, $C = 680 \text{ J kg}^{-1}$ for the heat capacity, $A = 0.02$ for the Bond albedo and $\epsilon_{IR} = 0.95$ for the thermal emissivity parameter.

We can use a standard algorithm for the calculation of the Yarkovsky acceleration which results in a semimajor-axis drift in case of non-resonant bodies. The drift in eccentricity in case of resonant bodies arises ‘automatically’ due to the gravitational part of the integrator. In Fig. 5 we can see a comparison between the expected drift $\Delta a$ in semimajor axis and the resulting drift $\Delta e$ in eccentricity, computed for the Hilda family (see the explanation in appendix A of Brož & Vokrouhlický 2008). The data can be approx-

<figure><img src="image_6.png" /><figcaption>Figure 5. An almost linear relation between the expected drift $\Delta a$ in semimajor axis and the simulated drift $\Delta e$ in eccentricity, computed for 360 members of the Hilda family located inside the J3/2 resonance.</figcaption></figure>

imated by a linear relationship, where the departures from linearity are caused mainly by interactions of drifting orbits with embedded weak secular or secondary resonances.

Note that according to a standard solar model the young Sun was faint (Güdel 2007), i.e. its luminosity 4 Gyr ago was 75 per cent of the current $L_{\odot}$. We can then expect a lower insolation and consequently weaker thermal effects acting on asteroids. Since we assume a constant value of $L_{\odot}$ in our code, the age estimated for the Hilda family (in Section 3.5) can be 12.5 per cent larger.

## 3.3 YORP effect

The magnitude of the Yarkovsky drift sensitively depends on the orientation of the spin axis with respect to the orbital plane and, to a lesser extent, on the angular velocity too. We thus have to account for the long-term evolution of spins of asteroids which is controlled by torques arising from the emission of thermal radiation, i.e. the YORP effect. The implementation of the YORP effect follows Vokrouhlický et al. (2006b). We assume the following relations for the rate of angular velocity and obliquity:

$$ \frac{d\omega}{dt} = f_i(\epsilon), \quad i = 1, \dots, 200, \qquad (5) $$

$$ \frac{d\epsilon}{dt} = \frac{g_i(\epsilon)}{\omega}, \qquad (6) $$

where $f$- and $g$-functions are given by Čapek & Vokrouhlický (2004) for a set of 200 shapes with mean radius $R_0 = 1 \text{ km}$, bulk density $\rho_0 = 2500 \text{ kg m}^{-3}$, located on a circular orbit with semimajor axis $a_0 = 2.5 \text{ au}$. The shapes of the Hilda family members are not known, so we assign one of the artificial shapes (denoted by the index $i$) randomly to each individual asteroid. We only have to scale the $f$- and $g$-functions by the factor

$$ c = c_{\text{YORP}} \left(\frac{a}{a_0}\right)^{-2} \left(\frac{R}{R_0}\right)^{-2} \left(\frac{\rho_{\text{bulk}}}{\rho_0}\right)^{-1}, \qquad (7) $$

where $a$, $R$, $\rho_{bulk}$ are semimajor axis, radius and density of the simulated body, and $c_{YORP}$ is a free scaling parameter, which can account for an additional uncertainty of the YORP model. Because the values of $f$ and $g$ were computed for only a limited set of obliquities (with a step $\Delta\epsilon = 30^{\circ}$) we use interpolation by Hermite polynomials (Hill 1982) of the data in Čapek & Vokrouhlický (2004) to obtain a smooth analytical function each for $f_i(\epsilon)$ and $g_i(\epsilon)$.

If the angular velocity approaches a critical value,

$$ \omega_{\text{crit}} = \sqrt{\frac{8}{3} \pi G \rho_{\text{bulk}}}, \qquad (8) $$

we assume a mass-shedding event, so we keep the orientation of the spin axis and the sense of rotation, but we reset the orbital period $P = 2\pi/\omega$ to a random value from the interval (2.5, 9) h. We also change the assigned shape to a different one, since any change of shape may result in a different YORP effect.

The differential equations (5) and (6) are integrated numerically by a simple Euler integrator. The usual time-step is $\Delta t = 1000 \text{ yr}$. An example of the results computed by the spin integrator for the Hilda family is displayed in Fig. 6. The typical time-scale of the spin-axis evolution is $\tau_{YORP} \simeq 500 \text{ Myr}$. After $\simeq 3$ times $\tau_{YORP}$ most bodies have spin axes perpendicular to their orbits, what maximizes the Yarkovsky drift rate of eccentricity.
<figure><img src="image_7.png" /><figcaption>**Figure 6.** An example of the YORP-driven evolution of obliquities (namely a z-component of the spin-axis unit vector, top panel) and angular velocities $\omega$ (bottom panel) for the members of the synthetic ‘Hilda’ family. At the beginning, all values of $\omega$ were selected positive and spin axes were distributed isotropically. The evolution may force $\omega$ to become negative, which simply corresponds to an opposite orientation of the spin axis. The scaling parameter selected for this run was $c_{YORP} = 0.33$.</figcaption></figure>

## 3.4 Collisions and spin-axis reorientations

In principle, collisions may directly affect the size distribution of the synthetic Hilda family, but we neglect this effect because most of the asteroids are large enough to remain intact.

However, we include spin-axis reorientations caused by collisions. We use an estimate of the time-scale by Farinella, Vokrouhlický & Hartmann (1998):

$$ \tau_{\text{reor}} = B \left( \frac{\omega}{\omega_0} \right)^{\beta_1} \left( \frac{D}{D_0} \right)^{\beta_2}, \quad (9) $$

where $B = 84.5$ kyr, $\beta_1 = 5/6$, $\beta_2 = 4/3$, $D_0 = 2$ m and $\omega_0$ corresponds to period $P = 5$ h. These values are characteristic of the main belt and we use them as an upper limit of $\tau_{\text{reor}}$ for the J3/2 region. Even so, the time-scale is $\tau_{\text{reor}} \simeq 3$ Gyr for the smallest observable ($D \simeq 5$ km) bodies, and reorientations are thus only of minor importance. Note that the probability of the reorientation is enhanced when the YORP effect drives the angular velocity $\omega$ close to zero.

## 3.5 Results for the Yarkovsky/YORP evolution

We start a simulation with an impact disruption of the parent body and create 360 fragments. Subsequent evolution of the synthetic Hilda family due to the Yarkovsky/YORP effect is computed up to 6 Gyr in order to estimate the time-span needed to match the observed family even though the family cannot be older than $\simeq 4$ Gyr, of course. Planets are started on their current orbits. A typical outcome of the simulation is displayed in Fig. 7.

<figure><img src="image_8.png" /><figcaption>**Figure 7.** Eccentricity versus absolute magnitude plot for the synthetic ‘Hilda’ family just after the impact disruption (time $t = 0$, top panel) and after 4 Gyr of evolution due to the Yarkovsky/YORP effect (bottom panel). There is a comparison with the observed Hilda family (grey dots).</figcaption></figure>

Due to the long integration time-span and large number of bodies, we were able to compute only four simulations with the following values of true anomaly and YORP efficiency:

(i) $f_{\text{imp}} = 0^{\circ}, c_{\text{YORP}} = 0$;

(ii) $f_{\text{imp}} = 180^{\circ}, c_{\text{YORP}} = 0$;

(iii) $f_{\text{imp}} = 0^{\circ}, c_{\text{YORP}} = 1$;

(iv) $f_{\text{imp}} = 0^{\circ}, c_{\text{YORP}} = 0.33$.

The remaining parameters were fixed: $e_i = 0.14$, $i_i = 7.8$, $\omega_{\text{imp}} = 30^{\circ}$, $\alpha = 3.25$, $v_{\text{max}} = 300 \text{ m s}^{-1}$, $R_{\text{PB}} = 93.5$ km, $\rho_{\text{PB}} = 1300 \text{ kg m}^{-3}$, $p_V = 0.044$.

We are mainly concerned with the distribution of eccentricities $e_p$, because the observed family has a large spread of $e_p$ values, while the initial synthetic family is very compact. For this purpose we constructed a Kolmogorov–Smirnov test (Press et al. 1999) of the normalized cumulative distributions $N(<e)$:

$$ D_{\text{KS}} = \max_{0 \le e < 1} |N(<e)_{\text{syn}} - N(<e)_{\text{obs}}|, \quad (10) $$

which provides a measure of the difference between the synthetic Hilda family, at a given time, and the observed Hilda family (see Fig. 8 for an example). The results of the KS tests are summarized in Fig. 9 (the first four panels).

There is an easy possibility to assess the sensitivity of results with respect to the $v_{\text{max}}$ parameter too, without the need to compute the simulation again. We simply select bodies fulfilling the condition $v < v'_{\text{max}}$, with $v'_{\text{max}} = 200, 100$ or $50 \text{ m s}^{-1}$, and recompute only the KS statistics for this subset. The results are plotted in Fig. 9 as thin lines. We can state that values lower than $v_{\text{max}} \simeq 100 \text{ m s}^{-1}$ are surely excluded.

As a preliminary conclusion we may say that all simulations point to a large age of the Hilda family. The $e_p$-distributions are most compatible with the observed family for ages $t = (4.0 \pm 1.0)$ Gyr. This suggests that the Hilda family might have experienced the giant-planet migration period which is dated by the late heavy bombardment to $t_{\text{LHB}} \simeq 3.85$ Gyr (Gomes et al. 2005). The large
<figure><img src="image_9.png" /><figcaption>**Figure 8.** Normalized cumulative distributions N(<e) of eccentricities for (i) the observed Hilda family, (ii) the synthetic ‘Hilda’ family at time t = 0 (just after the impact disruption), (iii) evolved due to the Yarkovsky/YORP effect (at time t = 3845 Myr). In this figure we show the best fit for the simulation with parameters f<sub>imp</sub> = 0°, c<sub>YORP</sub> = 0.33. Note that the ‘bent’ shape of the observed distribution corresponds to the ‘ears’ on the (e<sub>p</sub>, H) plot (Fig. 2). There is no perturbation by planetary migration in this particular case.</figcaption></figure>

uncertainty of the age stems from the fact that the runs including the YORP effect ($c_{YORP} \ge 0.33$) tend to produce ages at a lower limit of the interval while the YORP-less runs (with $c_{YORP} = 0$) tend to the upper limit.

## 3.6 Alternative hypothesis: high albedos of small asteroids

We now discuss two scenarios that further reduce the minimal age of the family: (i) high albedos of small asteroids (i.e. larger Yarkovsky/YORP drift); (ii) strongly asymmetric velocity field after impact (like that of the Veritas family).

Albedo is the most important unknown parameter, which can affect results on the Yarkovsky/YORP evolution. Fernández, Jewitt & Ziffer (2009) measured albedos of small Trojan asteroids and found systematically larger values than those for large Trojans. If we assume that the J3/2 asteroids behave similarly to Trojans, we may try a simulation with a rather high value of geometric albedo $p_v = 0.089$ (compared to previous $p_v = 0.044$). Moreover, we decrease density $\rho_{bulk} = 1200 \text{ kg m}^{-3}$, increase maximum velocity of fragments $v_{max} = 500 \text{ m s}^{-1}$ (though the velocity distribution is still determined by equation 4) and select true anomaly $f_{imp} = 90^\circ$ to maximize the spread of $e_p$ values.

The KS test is included in Fig. 9, panel (e). The most probable age is (2.3 ± 0.5) Gyr in this case. However, we do not think that the size-dependent albedo is very plausible because both large and small family members should originate from the same parent body and their albedos, at least just after the disruption, should be similar. Nevertheless, the albedos may change to a certain degree due to space weathering processes (Nesvorný et al. 2005). Unfortunately, we do not have enough data for small asteroids to assess a possible albedo difference between large and small family members.

## 3.7 Alternative hypothesis: strongly asymmetric velocity field

Another possibility to reduce the estimate of the family age is that the original velocity was highly anisotropic. A well-known example from the main belt is the Veritas family. Let us assume that the

anisotropy is of the order of Veritas, i.e. approximately four times larger in one direction. Note that Veritas is a young family and can be modelled precisely enough to compensate for chaotic diffusion in resonances (Nesvorný et al. 2003; Tsiganis, Knežević & Varvoglis 2007). This family is characterized by a large spread of inclinations, which corresponds to large out-of-plane components of velocities. In case of the Hilda family we multiply by 4 the radial components of initial velocities to maximize the dispersion of eccentricities, assuming the most favourable geometry of disruption ($f_{imp} \approx 90^\circ$).

The fit in Fig. 9, panel (f), is seemingly better at the beginning of the simulation, but bodies on unstable orbits are quickly eliminated and the fit gets much worse at $t \approx 500$ Myr. We can see that the synthetic Hilda family is similar to the observed Hilda family quite early (at $t \approx 2.5$ Gyr); however, the best fit is at later times ($t \approx 3.5$ Gyr), so there is no significant benefit compared to isotropic velocity-distribution cases.

# 4 DISRUPTION RATES IN THE J3/2 POPULATION

## 4.1 Present collisional activity

The results presented above show that the Hilda family is old. However, the uncertainty of the age is too large to conclude whether or not the family formed during the late heavy bombardment (LHB) period. An alternative constraint is the collisional lifetime of the parent body. If the probability that the parent body broke in the last 4 Gyr in the current collisional environment is negligible, it would argue that the family broke during the LHB when the collisional bombardment was much more severe. Thus, here we estimate the collisional lifetime of the parent body.

In our case, the target (parent body) has a diameter $D_{target} = 200$ km, a mean impact velocity $V_i = 4.8 \text{ km s}^{-1}$ (Dahlgren 1998), and a probable strength $Q_D^* = 4 \times 10^5 \text{ J kg}^{-1}$ (Benz & Asphaug 1999), and thus the necessary impactor size (Bottke et al. 2005) is

$$d_{disrupt} = (2Q_D^*/V_i^2)^{1/3} D_{target} \approx 65 \text{ km.} \quad (11)$$

The population of $\ge 65$ km projectiles is dominated by main-belt bodies: $n_{project} \approx 160$, according to Bottke et al. (2006), and we have only one 200-km target in the J3/2 region, so $n_{target} = 1$. The intrinsic collisional probability for Hilda versus main belt collisions is $P_i = 6.2 \times 10^{-19} \text{ km}^{-2} \text{ yr}^{-1}$ (Dahlgren 1998) and the corresponding frequency of disruptions is

$$f_{disrupt} = P_i \frac{D_{target}^2}{4} n_{project} n_{target} \approx 10^{-12} \text{ yr}^{-1}. \quad (12)$$

Thus, over the age of the Solar system $T_{SS} \approx 4$ Gyr (after LHB), we expect a very small number of such events $n_{events} = T_{SS}f_{disrupt} \approx 0.004$.

The value of strength $Q_D^*$ used above corresponds to strong targets. Though there is a theoretical possibility that the Hilda parent body was weaker, it does not seem to us likely, because the Hilda family is of the C taxonomic type. Thus, it is rather similar to (presumably stronger) main belt asteroids, than to (likely weaker) D-type objects. Anyway, even if we use an order of magnitude lower strength inferred for weak ice, $Q_D^* \approx 4 \times 10^4 \text{ J kg}^{-1}$ (see Leinhardt & Stewart 2009; Bottke et al. 2010), we obtain $d_{disrupt} \approx 30$ km, $n_{project} \approx 360$ and $n_{events} \approx 0.009$, so the conclusion about the low number of expected families remains essentially the same.
<figure><img src="image_10.png" /><figcaption>**Figure 9.** Kolmogorov–Smirnov tests of the synthetic ‘Hilda’ family: (a) no migration, only initial disruption (at anomaly $f_{\text{imp}} = 0^{\circ}$, $\omega_i = 30^{\circ}$) and subsequent Yarkovsky evolution; (b) the case with $f_{\text{imp}} = 180^{\circ}$; (c) including the YORP effect; (d) YORP with efficiency factor $c_{\text{YORP}} = 0.33$; (e) high albedo values (i.e. small bodies); (f) strongly asymmetric velocity field. The horizontal line denotes the distance $D_{\text{KS}} = 0.165$ for which the probability $p(>D_{\text{KS}})$ that the two eccentricity distributions differ by this amount equals 0.01.</figcaption></figure>

## 4.2 The late heavy bombardment

We now compute the probability that the parent body broke during the LHB. We can think of two projectile populations: (i) transient decaying cometary disc; (ii) D-type asteroids captured in the J3/2. Models like that of Levison et al. (2009) suggest that the decay time-scale of the cometary bombardment is of the order of 10–100 Myr and the flux of impactors integrated over this time-span might have been 100 times larger than it is today. Higher mean collisional velocities, due to projectiles on high-*e* and high-*i* orbits, are also favourable.

In order to estimate collisional activity we use a self-consistent model of the cometary disc from Vokrouhlický, Nesvorný & Levison (2008). Their *N*-body simulations included four giant planets and 27 000 massive particles with a total mass of $M_{\text{disc}} = 35 \, \text{M}_{\oplus}$.

The orbital evolution was propagated by the SyMBA integrator for 100 Myr. Using the output of these simulations, we calculate the mean intrinsic collisional probabilities $P_i(t)$ between the cometary-disc population (at given time $t$) and the current J3/2 population. We use an algorithm described in Bottke et al. (1994) for this purpose. Typically, the $P_i$ reaches $2\text{to } 3 \times 10^{-21} \, \text{km}^{-2} \, \text{yr}^{-1}$ and the corresponding mean impact velocities are $V_{\text{imp}} = 7–10 \, \text{km s}^{-1}$ (see Fig. 10).

The necessary impactor size is slightly smaller than before, $d_{\text{disrupt}} = 40–50 \, \text{km}$ due to larger $V_{\text{imp}}$. To estimate the number of such projectiles we assume that the cometary disc had a size distribution described by a broken power law with differential slopes $q_1 = 5.0$ for $D > D_0$, $q_2 = 2.5 \pm 0.5$ for $D < D_0$, where the diameter corresponding to the change of slopes is $D_0 = 50–70 \, \text{km}$. We then use the following expressions to calculate the number of bodies
<figure><img src="image_11.png" /><figcaption>**Figure 10.** Mean intrinsic collisional probability $P_i$ and mean impact velocity $V_{imp}$ versus time for one of the disc simulations from Vokrouhlický et al. (2008).</figcaption></figure>

larger than the given threshold (Vokrouhlický et al. 2008):

$$D_1 = D_0 \left[ \frac{(q_1 - 4)(4 - q_2)}{(q_1 - 1)(q_1 - q_2)} \frac{M_{\text{disc}}}{M_0} \right]^{\frac{1}{q_1 - 1}}, \quad (13)$$

$$N(>D) = \frac{q_1 - 1}{q_2 - 1} \left(\frac{D_1}{D_0}\right)^{q_1 - 1} \left(\frac{D_0}{D}\right)^{q_2 - 1} - \frac{q_1 - q_2}{q_2 - 1} \left(\frac{D_1}{D_0}\right)^{q_1 - 1}, \quad (14)$$

where $M_0 = \frac{\pi}{6} \rho D_0^3$ and $\rho = 1300 \text{ kg m}^{-3}$. The result of this calculation is $N(>d_{\text{disrupt}}) \approx 0.3$ to $1.7 \times 10^9$. The actual number of bodies in the simulation (27 000) changes in course of time and it was scaled such that initially it was equal to $N(>d_{\text{disrupt}})$. The resulting number of events is

$$n_{\text{events}} = \frac{D_{\text{target}}^2}{4} n_{\text{target}} \int P_i(t) n_{\text{project}}(t) dt \\ \approx 0.05 \text{ to } 0.2, \quad (15)$$

which is 10–50 times larger than the number found in Section 4.1.

Regarding the captured D-type asteroids, they were probably not so numerous and their impact velocities were lower but their collisional probabilities were larger and the population might have had a substantially longer time-scale of decay (Levison et al. 2009). Using the reasonable values $V_i = 4.0 \text{ km s}^{-1}$, $d_{\text{disrupt}} = 70 \text{ km}$, $n_{\text{project}} = 5000$, $P_i = 2.3 \times 10^{-18} \text{ km}^{-2} \text{yr}^{-1}$, $T_{\text{LHB}} \approx 1 \text{ Gyr}$, we obtain the number of events $\approx 0.1$ which is again 25 times larger than the number presented in Section 4.1.

We conclude that the Hilda family was likely created during the LHB when the collisions were much more frequent than in the current collisional environment. We must now test whether the structure of the family is consistent with the giant-planet migration, since it is connected with the LHB.

## 5 PLANETARY MIGRATION

At the LHB-time the planetary migration was most probably caused by the presence of a massive cometary disc. Instead of a full N-body model we use a simpler analytic migration, with an artificial dissipation applied to the planets. This is the only realistic possibility in our case, because we need to test not only a large number of various migration scenarios but also various initial configurations of the synthetic Hilda family.

For this purpose we use a modified version of the symplectic SWIFT-RMVS3 integrator (Levison & Duncan 1994). We account for four giant planets and include the following dissipation term

<table><caption>**Table 1.** Free parameters of our Hilda family model.</caption><thead><tr><th>No.</th><th>Parameter</th><th>Description</th></tr></thead><tbody><tr><td>1</td><td>a<sub>Ji</sub></td><td>Initial semimajor axis of Jupiter</td></tr><tr><td>2</td><td>a<sub>Si</sub></td><td>Saturn</td></tr><tr><td>3</td><td>e<sub>Ji</sub></td><td>Initial eccentricity of Jupiter</td></tr><tr><td>4</td><td>e<sub>Si</sub></td><td>Saturn</td></tr><tr><td>5</td><td>τ<sub>mig</sub></td><td>Migration time-scale</td></tr><tr><td>6</td><td>e<sub>dampJ</sub></td><td>Eccentricity damping for Jupiter</td></tr><tr><td>7</td><td>e<sub>dampS</sub></td><td>Saturn</td></tr><tr><td>8</td><td>e<sub>i</sub></td><td>Initial eccentricity of the parent body</td></tr><tr><td>9</td><td>i<sub>i</sub></td><td>Initial inclination</td></tr><tr><td>10</td><td>f<sub>imp</sub></td><td>True anomaly at the impact disruption</td></tr><tr><td>11</td><td>ω<sub>imp</sub></td><td>Argument of perihelion</td></tr><tr><td>12</td><td>α</td><td>Slope of the velocity distribution</td></tr><tr><td>13</td><td>v<sub>max</sub></td><td>Maximum velocity of fragments</td></tr><tr><td>14</td><td>R<sub>PB</sub></td><td>Radius of parent body</td></tr><tr><td>15</td><td>ρ<sub>PB</sub></td><td>Bulk density</td></tr><tr><td>16</td><td>p<sub>V</sub></td><td>Geometric albedo of fragments</td></tr><tr><td>17</td><td>c<sub>YORP</sub></td><td>Efficiency of the YORP effect</td></tr></tbody></table>

<table><caption>**Table 2.** Fixed (assumed) parameters of the Hilda family model. There are also a number of less important parameters, such as the thermal ones ($\rho_{\text{surf}}$, K, C, A, $\epsilon_{\text{IR}}$) or the collisional ones ($B$, $\beta_1$, $\beta_2$).</caption><thead><tr><th>No.</th><th>Parameter</th><th>Description</th></tr></thead><tbody><tr><td>18</td><td>a<sub>Jf</sub></td><td>Final semimajor axis of Jupiter</td></tr><tr><td>19</td><td>a<sub>Sf</sub></td><td>Saturn</td></tr><tr><td>20</td><td>N(&lt;H)</td><td>(observed) absolute magnitude distribution</td></tr></tbody></table>

applied to the planets in every time-step:

$$v = v \left[ 1 + \frac{\Delta v}{v} \frac{\Delta t}{\tau_{\text{mig}}} \exp \left( -\frac{t - t_0}{\tau_{\text{mig}}} \right) \right], \quad (16)$$

where $v$ denotes a velocity vector of a given planet, $v$ is the absolute value of velocity, $\Delta t$ is the time-step, $\tau_{\text{mig}}$ is the selected migration time-scale, $\Delta v = \sqrt{GM/a_i} - \sqrt{GM/a_f}$ the required total change of velocity (i.e. the difference of mean velocities between the initial and the final orbit), $t$ is the time and $t_0$ is some reference time. If there are no perturbations other than (16), the semimajor axis of the planet changes smoothly (exponentially) from the initial value $a_i$ to the final $a_f$. We use time-step $\Delta t = 36.525 \text{ d}$ and the total time-span of the integration is usually equal to $3\tau_{\text{mig}}$ when planetary orbits practically stop to migrate.

We would like to resemble evolution of planetary orbits similar to the Nice model so it is necessary to use an eccentricity-damping formula, which simulates the effects of dynamical friction (Morbidelli et al. 2010). This enables us to model a decrease of eccentricities of the giant planets to relatively low final values. The amount of eccentricity damping is characterized by a parameter $e_{\text{damp}}$.

Because inclinations of the planets are not very important for what concerns the perturbation of minor bodies (the structure of resonances is mainly determined by planetary eccentricities), we usually start with the current values of inclination of the planets.

We admit that the analytic migration is only a crude approximation of the real evolution, but we can use it as a first check to see which kinds of migration scenarios are allowed and which are not with respect to the existence and structure of the Hilda family.

As a summary we present a list of free and fixed (assumed) parameters of our model in Tables 1 and 2. According to our numerical tests the initial configuration of Uranus and Neptune is not very important, as these planets do not produce significant direct
perturbations on asteroids located in the J3/2 resonance. We thus
do not list the initial semimajor axes and eccentricities of Uranus
and Neptune among our free parameters though we include these
planets in our simulations.

The problem is that we cannot tune all the 17 parameters together,
since the 17-dimensional space is enormous. We thus first select a
reasonable set of impact parameters for the family (No. 8–17 in
Table 1), keep them fixed, and experiment with various values of
migration parameters (No. 1–7). We test roughly 10³ migration
scenarios. Then, in the second step, we vary impact parameters for
a single (successful) migration scenario and check the sensitivity of
results.

## 5.1 Results for planetary migration

In the first test we compute an evolution of the synthetic Hilda
family during planetary migration phase for the following parameter
space (these are not intervals but lists of values): $a_{Ji} = (5.2806$ and
$5.2027)$ au, $a_{Si} = (8.6250, 8.8250, 9.3000)$ au, $e_{Ji} = (0.065, 0.045)$,
$e_{Si} = (0.08, 0.05)$, $\tau_{mig} = (0.3, 3, 30, 300)$ Myr, $e_{dampJ} = 10^{-11}$,
$e_{damps} = 10^{-11}$.$^1$ The values of $a_{Ji}$ and $a_{Si}$ correspond to period
ratios $P_S/P_J$ from 2.09 to 2.39 (the current value is 2.49), i.e. the
giant planets are placed already *beyond* the 2:1 resonance, since
the 2:1 resonance crossing would destroy the Hilda family (Brož
& Vokrouhlický 2008). Impact parameters were fixed except $f_{imp}$:
$e_i = 0.14, i_i = 7.8, f_{imp} = (0^{\circ}, 180^{\circ})$, $\omega_{imp} = 30^{\circ}$, $\alpha = 3.25$, $v_{max} = 300$ m s$^{-1}$, $R_{PB} = 93.5$ km, $\rho_{PB} = 1300$ kg m$^{-3}$.

The synthetic Hilda family has 360 bodies in case of short simu-
lations ($\tau_{mig} = 0.3\text{or}3$ Myr). In case of longer simulations we create
60 bodies only. Their absolute magnitudes (sizes) were thus selected
randomly from 360 observed values. This is a minimum number of
bodies necessary to compare the distributions of eccentricities. We
performed tests with larger numbers of bodies and the differences
do not seem significant.

A comparison of the final orbits of the planets with the current planetary orbits shows we have to exclude some migration simulations (mostly those with Uranus and Neptune on compact orbits). One of the reasons for the unsuccessful scenarios is that a compact configuration of planets is inherently unstable. If the migration time-scale is too large or the eccentricity damping too low, it may result in a violent instability, close encounters between planets and eventually an unrealistic final configuration.

The change in the structure of the synthetic Hilda family due to
migration can be seen in Fig. 11. The family is shifted in semimajor
axis, because it moves together with the resonance with migrating
Jupiter. Moreover, the eccentricities are dispersed while the incli-
nations are barely affected.

We have found that the eccentricity distribution is modified when
secondary resonances occur between the libration frequency $f_{J3/2}$ of an asteroid in the J3/2 resonance and the frequency $f_{1J-2S}$ of the critical argument of Jupiter–Saturn 1:2 resonance (see Kortenkamp, Malhotra & Michtchenko 2004; Morbidelli et al. 2005 for the case of Trojans):

$$nf_{J3/2} = f_{1J-2S}, \quad (17)$$

<figure><img src="image_12.png" /><figcaption>**Figure 11.** A usual evolution of the synthetic ‘Hilda’ family in the pseudo-proper semimajor axis versus eccentricity plot. The initial ($t = 0$ Myr) and final stages ($t = 100$ Myr) are plotted. The migration time-scale was $\tau_{mig} = 30$ Myr in this particular example. We selected this longer time-scale because secular frequencies can then be computed more precisely (see Fig. 12). The arrow indicates a total change of the position of the J3/2 resonance due to migration of Jupiter.</figcaption></figure>

where *n* is a small integer number, *n* = 2, 3 or 4 in our case.<fnref n="2" /> We can
see the evolution of resonant semimajor axes and the correspond-
ing dominant frequencies, computed by means of periodogram, in
Fig. 12.

Because the resonances are localized – they act only at particular
values of semimajor axes of planets – it is *not* necessary to have a
dense grid in $a_{Ji}$, $a_{Si}$ parameters to study the dependence of the syn-
thetic Hilda family shape on $a_{Ji}$, $a_{Si}$. Essentially, there are only three
situations when the Hilda family is strongly perturbed, otherwise
the spread in *e* does not change much in course of time.

A very simple test, which allows us to quickly select allowed migration scenarios, is the number of remaining Hilda family members. We may assume that the depletion by dynamical effects was probably low (say 50 per cent at most), otherwise we would obtain a much larger parent body than $D \simeq 200$ km, which has a much lower probability of collisional disruption. The fractions of the remaining bodies $N_{left}/N_{initial}$ versus initial conditions for planets are displayed in Fig. 13.

The small number of remaining bodies $N_{left}$ indicates that perturbations acting on the synthetic family were too strong. It means either the family had to be formed later (when fewer and weaker secondary resonances are encountered) to match the observed family or this migration scenario is not allowed. The same applies to the dispersion of *e*-distribution (see below): if it is too large compared to the observed Hilda family, the synthetic Hilda had to be formed later or the scenario is not allowed. Our results indicate that

(i) a faster migration time-scale $\tau_{mig} \simeq 0.3-30$ Myr is preferred
over slower time-scales;

(ii) Jupiter and Saturn were not in the most compact configuration ($a_{Ji} = 5.2806$ au, $a_{Si} = 8.6250$ au) at the time when the Hilda family was created.

## 5.2 A sensitivity to the impact-related parameters

Another important test was devoted to the impact parameters, which
were varied in relatively larger steps: $e_i = (0.12, 0.15)$, $i_i = (6^{\circ}, 8$,
$8^{\circ}, 8)$, $f_{imp} = (45^{\circ}, 90^{\circ}, 135^{\circ})$, $\omega_{imp} = (60^{\circ}, 90^{\circ})$, $\alpha = (2.25, 4.25)$,
$v_{max} = (200, 400)$ m s$^{-1}$, $R_{PB} = (83.5, 103.5)$ km, $\rho_{PB} = (1000$,
$2000)$ kg m$^{-3}$. Note that the selected impact parameters are rather
extreme, reason that we do not expect them to ever be out of these

<footnote n="1">¹ In order to increase the statistics we ran simulations multiple times with different initial conditions for Uranus and Neptune: $a_{Ui} = (18.4479, 12.3170)$ au, $a_{Ni} = (28.0691, 17.9882)$ au, $e_{Ui} = (0.06, 0.04)$, $e_{Ni} = (0.02, 0.01)$.</footnote>

<footnote n="2">² We also looked for secondary resonances connected with the 4:9, 3:7 and 2:5 Jupiter–Saturn resonances, but we found no significant effects.</footnote>
<figure><img src="image_13.png" /><figcaption>**Figure 12.** Top panel: the frequency $f_{1J-2S}$ of the Jupiter–Saturn 1:2 mean motion critical argument (thick grey curve) versus time $t$. The frequency changes due to the migration of planets with the time-scale $\tau_{\text{mig}} = 30$ Myr. We also computed dominant frequencies $f_{J3/2}$ of librations in the J3/2 resonance for three selected members of the synthetic Hilda family (black curves). We do not plot the frequency itself but a selected multiple of it $nf_{J3/2}$. Captures in the secondary resonances of type $nf_{J3/2} = f_{1J-2S}$ are then clearly visible when the frequencies are equal. For the test particle number 1 it occurs between 4 and 10 Myr, particle 2 was captured from 21 to 32 Myr and particle 3 from 54 Myr till the end of the simulation. Bottom panel: the corresponding changes of the pseudo-proper semimajor axes $a_p$ versus time $t$ due to the secondary resonances. The three test particles from the top panel are shown (black curves) together with the remaining members of synthetic ‘Hilda’ family (grey curves). Note that some particles may be pushed to the border of the stable libration zone and then escape from the J3/2 resonance.</figcaption></figure>

bounds. The total number of simulations is 384. The migration parameters were fixed (they correspond to one successful migration scenario): $a_{Ji} = 5.2806$ au, $a_{Si} = 8.8250$ au, $e_{Ji} = 0.065$, $e_{Si} = 0.08$, $\tau_{\text{mig}} = 3$ Myr, $e_{\text{dampJ}} = 10^{-11}$, $e_{\text{dampS}} = 10^{-11}$.

This time, we decided to use a simple quantity to discuss the results, namely the eccentricity dispersion $\sigma_e$ of the synthetic family at the end of the giant-planet migration. The most frequent values of the dispersion are $\sigma_e = 0.015–0.04$ (see the histograms in Fig. 14). Further evolution by the Yarkovsky/YORP effect would increase the dispersions up to $\sigma_e = 0.045–0.06$, while the observed dispersion of the Hilda family is $\sigma_e = 0.046$.

It is notable that the histograms look similar for all the impact parameters, there is even no apparent correlation between them. The explanation for this ‘lack of dependence’ is that the eccentricity distribution is mainly determined by the perturbations of the giant planets. A given planetary evolution therefore gives a characteristic value of $\sigma_e$ whatever the impact parameters are. The dispersion in $\sigma_e$ values is due to the fact that the planetary evolutions that we have computed change widely from one simulation to another. Though planet migration was prescribed analytically, there are mutual interactions between planets and random captures in resonances (or jumps across resonances) which may affect the eccentricity distribution of the synthetic Hilda family. An extreme case is shown in Fig. 15. In this particular simulation, Jupiter and Saturn were captured in the mutual 3:7 resonance for 0.5 Myr which resulted in a large eccentricity dispersion $\sigma_e = 0.044$ of the synthetic family. Our conclusion is that the impact parameters are less important than the parameters related to migration.

## 5.3 Matching results together

Even though we do not perform a joint integration which includes both the planetary migration and Yarkovsky/YORP effect, we try to match the previous results from Sections 5.1 and 3.5. We do it by using a straightforward Monte Carlo approach: (i) we take the pseudo-proper eccentricities $e_{\text{mig}}$ of bodies at the end of planetary migration from Section 5.1; (ii) we compute the total Yarkovsky/YORP drifts $\Delta e_{\text{YE}}$ in eccentricity from Section 3.5; (iii) we assign every body a drift randomly ($e_{\text{final}} = e_{\text{mig}} + \Delta e_{\text{YE}}$) and this way we construct an

<figure><img src="image_14.png" /><figcaption>**Figure 13.** The number of simulations *N* versus the fraction of remaining bodies $N_{\text{left}}/N_{\text{initial}}$ from the synthetic ‘Hilda’ family. The histograms are plotted for four different time-scales of migration $\tau_{\text{mig}}$ and six different initial configurations of Jupiter and Saturn ($a_{Ji}$, $a_{Si}$; we indicate period ratios $P_{Si}/P_{Ji}$ instead of semimajor axes here). The ranges of the remaining free parameters are mentioned in the text. We only plot the successful migration scenarios with $\Delta v_{\text{planets}} \le 2000 \text{ m s}^{-1}$, where $\Delta v_{\text{planets}} = \sum_{i=1}^{4} \delta v_i$ is a sum of the velocity differences $\delta v$ (defined similarly as in the HCM metric, equation 3) between the final simulated orbit of the *i*th planet and the currently observed one. This way we join differences in orbital elements *a*, *e*, *I* into a single quantity which has the dimension of velocity.</figcaption></figure>
<figure><img src="image_15.png" /><figcaption>**Figure 14.** Eccentricity dispersions of the synthetic 'Hilda' families at the end of the giant-planet migration for the various initial conditions of the impact disruption: initial eccentricity $e_i$, inclination $i_i$, true anomaly $f_{imp}$, argument of perihelion $\omega_{imp}$, exponent $\alpha$, maximum velocity $v_{max}$, radius of the parent body $R_{PB}$ and its bulk density $\rho_{PB}$. The values of the remaining parameters related to migration are mentioned in the text. Note that there is no evolution by the Yarkovsky/YORP effect in this simulation. The dotted vertical line denotes the value $\sigma_e = 0.046$ of the observed Hilda family.</figcaption></figure>

<figure><img src="image_16.png" /><figcaption>**Figure 15.** An example of the orbital evolution of Jupiter and Saturn with a rare temporary capture in the mutual 3:7 resonance (bottom panel). This sort of evolution leads to a large spread of pseudo-proper eccentricities of the synthetic 'Hilda' family by the end of the migration (top panel).</figcaption></figure>

evolved synthetic family.<fnref n="3" /> Finally, we compare the synthetic family to the observed Hilda family by computing a Kolmogorov–Smirnov test for $N(<e_{final})$ and $N(<e_{obs})$ distributions.

To avoid problems with the small number of bodies (60 in case of planetary migration), we perform the above procedure 100 times, always with a different random seed for the assignment of the $\Delta e_{YE}$. We then take a median of the 100 KS statistics as a result for one particular run. The resulting histograms of the median $D_{KS}$ for various initial conditions are shown in Fig. 16.

We confirm the conclusions from Section 5.1 – those migration scenarios that preserve the largest number of family members (i.e.

high $N_{left}$) are the same, for which we can find a good fit of eccentricity distribution (low $D_{KS}$). Moreover, it seems we can exclude also the time-scale of migration $\tau_{mig} = 30$ Myr since the total number of successful simulations is significantly smaller in this case.

## 6 CONCLUSIONS

Results of this paper can be summarized as follows.

(i) The Hilda family evolves mainly due to the Yarkovsky/YORP effect and the observed large spread of eccentricities indicates the age $4^{+0}_{-0}$ Gyr.

(ii) The collisional disruption of a $D \simeq 200$ km parent body is unlikely in the current environment. Instead, it rather occurred during the LHB when collisions with comets dominated and were up to 50 times more frequent. Another possible source of projectiles is the population of D-type asteroids captured in the J3/2 resonance (Levison et al. 2009).

(iii) In case the Hilda family was created during giant-planet migration, which seems to us likely, the major perturbations of the family were due to secondary resonances between libration frequency and the frequency of Jupiter–Saturn 1:2 critical argument.

(iv) On the basis of our simulations, we argue that the migration was relatively fast (with time-scale $\tau_{mig} \simeq 0.3–3$ Myr) and Jupiter and Saturn were relatively closer to the current configuration (with period ratio $P_S/P_J \ge 2.13$ or more) at the moment when the Hilda family was created, otherwise the family would be ‘destroyed’ by migration. Slower migration time scales are only allowed for larger values of $P_S/P_J$ ratios.

The Hilda family thus proved to be one of the oldest families in the main asteroid belt.

There are emerging indications that orbital evolution of planets was rather violent and close encounters between planets were present (Nesvorný et al. 2007; Brasser et al. 2009). This might be still consistent with our model of the Hilda family, but of course we have to assume that the family formed after severe perturbations

<footnote n="3">³ Note that gravitational perturbations, caused by planetary migration, are independent of size (mass), so a large body may be easily found at the outskirts of the family. This is another reason for the random assignment of Yarkovsky/YORP drifts.</footnote>
<figure><img src="image_17.png" /></figure>

**Figure 16.** The number of simulations *N* versus the Kolmogorov–Smirnov distance *D*<sub>KS</sub> between the synthetic and the observed Hilda family. The simulations differ by the time-scale of migration τ<sub>mig</sub> and the initial conditions for Jupiter and Saturn (*a*<sub>Ji</sub>, *a*<sub>Si</sub>). We only plot the successful migration scenarios with Δ*v*<sub>planets</sub> ≤ 2000 m s<sup>-1</sup> and the number of bodies left *N*<sub>left</sub> > *N*<sub>initial</sub>/2. The dotted vertical line denotes the distance *D*<sub>KS</sub> for which the probability *p*(*v* > *D*<sub>KS</sub>) that the two eccentricity distributions differ by this amount equals 0.01.

in the J3/2 region ended. A more complicated migration scenario like that of 'jumping Jupiter' (Morbidelli et al. 2010) even seems favourable in our case because Jupiter and Saturn very quickly reach a high period ratio (P<sub>S</sub>/P<sub>J</sub> ≳ 2.3, i.e. the planets are quite close to their current orbits). Then, the perturbations acting on the J3/2 region are already small and the flux of impactors becomes high just after the jump. The Hilda family thus might have formed exactly during this brief period of time.

Regarding future improvements of our model, knowledge of ge-
ometric albedos for a large number of small asteroids may signifi-
cantly help and decrease uncertainties. The WISE infrared mission
seems to be capable of obtaining these data in near future.

ACKNOWLEDGMENTS

We thank Hal Levison for his code on eccentricity damping, David Čapek for sending us the YORP effect data in an electronic form and an anonymous referee for constructive comments.

The work of MB and DV has been supported by the Grant Agency of the Czech Republic (grants 205/08/P196 and 205/08/0064) and the Research Program MSM0021620860 of the Czech Ministry of Education. We also acknowledge the provision of computers to us at the Observatory and Planetarium in Hradec Králové.

REFERENCES

Benz W., Asphaug E., 1999, Icarus, 142, 5

Bottke W. F., Nolan M. C., Greenberg R., Kolvoord R. A., 1994, Icarus,
107, 255

Bottke W. F., Vokrouhlický D., Brož M., Nesvorný D., Morbidelli A., 2001, Sci, 294, 1693

Bottke W. F., Durda D. D., Nesvorný D., Jedicke R., Morbidelli A., Vokrouhlický D., Levison H. F., 2005, Icarus, 175, 111

Bottke W. F., Nesvorný D., Grimm R. E., Morbidelli A., O'Brien D. P., 2006, Nat, 439, 821

Bottke W. F., Nesvorný D., Vokrouhlický D., Morbidelli A., 2010, AJ, 139, 994

Brasser R., Morbidelli A., Gomes R., Tsiganis K., Levison H. F., 2009, A&A, 507, 1053

Brož M., 2006, PhD thesis, Charles Univ.

Brož M., Vokrouhlický D., 2008, MNRAS, 390, 715

Čapek D., Vokrouhlický D., 2004, Icarus, 172, 526

Dahlgren M., 1998, A&A, 336, 1056

Davis D. R., Neese C., eds, 2002, Asteroid Albedos. EAR-A-5-DDR-ALBEDOS-V1.1. NASA Planetary Data System

Dohnanyi J. W., 1969, J. Geophys. Res., 74, 2531

Farinella P., Froeschlé C., Gonczi R., 1994, in Milani A., Di Martino M., Cellino A., eds, Asteroids, Comets, Meteors 1993. Kluwer, Dordrecht, p. 205

Farinella P., Vokrouhlický D., Hartmann W. K., 1998, Icarus, 132, 378

Fernández Y. R., Jewitt D., Ziffer J. E., 2009, AJ, 138, 240

Ferraz-Mello S., Michtchenko T. A., Nesvorný D., Roig F., Simula A., 1998, P&SS, 46, 1425

Gomes R., Levison H. F., Tsiganis K., Morbidelli A., 2005, Nat, 435, 466

Güdel M., 2007, Living Rev. Solar Phys., 4, 3

Hill G., 1982, Publ. Dom. Astrophys. Obser. Victoria BC, 16, 67

Kortenkamp S. J., Malhotra R., Michtchenko T., 2004, Icarus, 167, 347

Leinhardt Z. M., Stewart S. T., 2009, Icarus, 199, 542

Levison H. F., Duncan M., 1994, Icarus, 108, 18

Levison H. F., Morbidelli A., Vanlaerhoven Ch., Gomes R., Tsiganis K., 2008, Icarus, 196, 258

Levison H. F., Bottke W. F., Gounelle M., Morbidelli A., Nesvorný D., Tsiganis K., 2009, Nat, 460, 364

Malhotra R., 1995, AJ, 110, 420

Minton D. A., Malhotra R., 2009, Nat, 457, 1109

Morbidelli A., Moons M., 1993, Icarus, 102, 316

Morbidelli A., Levinson H. F., Tsiganis K., Gomes R., 2005, Nat, 435, 459

Morbidelli A., Brasser R., Tsiganis K., Gomes R., Levison H. F., 2009, A&A, 507, 1041

Morbidelli A., Brasser R., Gomes R., Levison H. F., Tsiganis K., 2010, AJ, 140, 1391

Nesvorný D., Ferraz-Mello S., 1997, Icarus, 130, 247

Nesvorný D., Bottke W. F., Levison H. F., Dones L., 2003, ApJ, 591, 486

Nesvorný D., Jedicke R., Whiteley R. J., Ivezić Ž., 2005, Icarus, 173, 132

Nesvorný D., Vokrouhlický D., Morbidelli A., 2007, AJ, 133, 1962

Parker A., Ivezić Ž., Jurić M., Lupton R., Sekora M. D., Kowalski A., 2008, Icarus, 198, 138

Press W. H., Teukolsky S. A., Vetterlink W. T., Flannery B. P., 1999, Numerical Recipes in Fortran 77. Cambridge Univ. Press, Cambridge

Quinn T. R., Tremaine S., Duncan M., 1991, AJ, 101, 2287

Roig F., Ferraz-Mello S., 1999, Planet. Space Sci., 47, 653

Roig F., Nesvorný D., Ferraz-Mello S., 2002, MNRAS, 335, 417

Tanga P., Cellino A., Michel P., Zappatà V., Paolicchi P., Dell'Oro A., 1998, Icarus, 141, 65

Tsiganis K., Gomes R., Morbidelli A., Levison H. F., 2005, Nat, 435, 459

Tsiganis K., Knežević Z., Varvoglis H., 2007, Icarus, 186, 484

Vokrouhlický D., Brož M., Bottke W. F., Nesvorný D., Morbidelli A., 2006a, Icarus, 182, 92

Vokrouhlický D., Brož M., Bottke W. F., Nesvorný D., Morbidelli A., 2006b, Icarus, 182, 118

Vokrouhlický D., Nesvorný D., Levison H. F., 2008, AJ, 136, 1463

Zappalà V., Cellino A., Farinella P., Milani A., 1994, AJ, 107, 772

This paper has been typeset from a TeX/LA TeX file prepared by the author.