Clarissa Family Age from the Yarkovsky Effect Chronology
Vanessa C. Lowry, David Vokrouhlicky, David Nesvorny, Humberto Campins
DDraft version September 15, 2020
Typeset using L A TEX default style in AASTeX63
Clarissa Family Age from the Yarkovsky Effect Chronology
Vanessa C. Lowry, David Vokrouhlick´y, David Nesvorn´y, and Humberto Campins University of Central Florida, 4000 Central Florida Blvd, Orlando, FL 32816, USA Astronomical Institute, Charles University, Prague, V Holeˇsoviˇck´ach 2, CZ 18000, Prague 8, Czech Republic Department of Space Studies, Southwest Research Institute, 1050 Walnut St., Suite 300, Boulder, CO 80302, USA
ABSTRACTThe Clarissa family is a small collisional family composed of primitive C-type asteroids. It is locatedin a dynamically stable zone of the inner asteroid belt. In this work we determine the formation age ofthe Clarissa family by modeling planetary perturbations as well as thermal drift of family members dueto the Yarkovsky effect. Simulations were carried out using the SWIFT-RMVS4 integrator modifiedto account for the Yarkovsky and Yarkovsky–O’Keefe–Radzievskii–Paddack (YORP) effects. We ranmultiple simulations starting with different ejection velocity fields of fragments, varying proportion ofinitially retrograde spins, and also tested different Yarkovsky/YORP models. Our goal was to matchthe observed orbital structure of the Clarissa family which is notably asymmetrical in the propersemimajor axis, a p . The best fits were obtained with the initial ejection velocities (cid:46)
20 m s − ofdiameter D (cid:39) ∼ (cid:39)
80% ofsmall family members initially had retrograde rotation. The age of the Clarissa family was found tobe t age = 56 ± ρ = 1 . − . Small variation of densityto smaller or larger value would lead to slightly younger or older age estimates. This is the first casewhere the Yarkovsky effect chronology has been successfully applied to an asteroid family youngerthan 100 Myr. INTRODUCTIONAsteroid families consist of fragments produced by catastrophic and cratering impacts on parent bodies (see Nesvorn´yet al. 2015 for a review). The fragments produced in a single collision, known as family members, share similar propersemimajor axes ( a p ), proper eccentricities ( e p ), and proper inclinations ( i p ) (Kneˇzevi´c et al. 2002; Carruba & Nesvorn´y2016). Family members are also expected to have spectra that indicate similar mineralogical composition to the parentbody (Masiero et al. 2015). After their formation, families experience collisional evolution (Marzari et al. 1995), whichmay cause them to blend into the main belt background, and evolve dynamically (Bottke et al. 2001). Since thecollisional lifetime of a 2 km sized body in the main belt is greater than 500 Myr (Bottke et al. 2005), which is nearly10 times longer than the Clarissa family’s age (Section 4), we do not need to account for collisional evolution. Instead,we only consider dynamical evolution of the Clarissa family to explain its current orbital structure and constrain itsformation conditions and age.The Clarissa family is one of small, primitive (C or B type in asteroid taxonomy) families in the inner asteroid belt(Morate et al. 2018). Of these families, Clarissa has the smallest extent in semimajor axis suggesting that it might bethe youngest. This makes it an interesting target of dynamical study. Before discussing the Clarissa family in depth,we summarize what is known about its largest member, asteroid (302) Clarissa shown in Figure 1. Analysis of denseand sparse photometric data shows that (302) Clarissa has a retrograde spin with rotation period P = 14 . − ◦ , and two possible solutions for pole ecliptic longitude, 28 ◦ or 190 ◦ (Hanuˇs et al. 2011). Thelatter solution was excluded once information from stellar occultations became available ( ˇDurech et al. 2010, 2011).The stellar occultation also provided a constraint on the volume-equivalent diameter D = 43 ± ± Corresponding author: Vanessa Lowryvanessa [email protected] a r X i v : . [ a s t r o - ph . E P ] S e p Lowry et al.
Figure 1.
Three-dimensional shape of (302) Clarissa from light-curve inversion ( ˇDurech et al. 2010, 2011; Hanuˇs et al. 2011;https://astro.troja.mff.cuni.cz/projects/damit/). Equatorial views along the principal axes of inertia are shown in the left andmiddle panels. The polar view along the shortest inertia axis is shown on the right (view from the direction of the rotationpole).
The 179 family members of the Clarissa family are tightly clustered around (302) Clarissa (Nesvorn´y et al. 2015).The synthetic proper elements of the Clarissa family members were obtained from the Planetary Data System (PDS)node (Nesvorn´y 2015). Figure 2 shows the projection of the Clarissa family onto the ( a p , e p ) and ( a p , sin i p ) planes.For comparison, we also indicate in Fig. 2 the initial distribution of fragments if they were ejected at speeds equal tothe escape speed from (302) Clarissa.To generate these initial distributions we adopted the argument of perihelion ω (cid:39) ◦ and true anomaly f (cid:39) ◦ ,both given at the moment of the parent body breakup (Zappal`a et al. 1984; Nesvorn´y et al. 2006a; Vokrouhlick´yet al. 2006a) (see Appendix A). Other choices of these (free) parameters would lead to ellipses in Figure 2 that wouldbe tilted in the ( a p , e p ) projection and/or vertically squashed in the ( a p , i p ) projection (Vokrouhlick´y et al. 2017a).Interestingly, the areas delimited by the green ellipses in Fig. 2 contain only a few known Clarissa family members.We interpret this as a consequence of the dynamical spreading of the Clarissa family by the Yarkovsky effect.Immediately following the impact on (302) Clarissa, the initial spread of fragments reflects their ejection velocities.We assume that the Clarissa family was initially much more compact than it is now (e.g., the green ellipses in Fig. 2).As the family members drifted by the Yarkovsky effect, the overall size of the family in a p slowly expanded. It isapparent that the Clarissa family has undergone Yarkovsky drift since there is a depletion of asteroids in the centralregion of the family in Figure 3.There are no major resonances in the immediate orbital neighborhood of (302) Clarissa. The e p and i p values offamily members therefore remained initially unchanged. Eventually, the family members reached the principal meanmotion resonances, most notably the J4-S2-1 three-body resonance at (cid:39) .
398 au (Nesvorn´y & Morbidelli 1998) whichcan change e p and i p . This presumably contributed to the present orbital structure of the Clarissa family, wheremembers with a p < .
398 au have significantly larger spread in e p and i p than those with a p > .
398 au. Note, inaddition, that there are many more family members sunward from (302) Clarissa, relative to those on the other side(Fig. 3). We discuss this issue in more detail in Section 4.Figure 3 shows the absolute magnitudes H of family members as a function of a p . We use the mean WISE albedo ofthe Clarissa family, p V = 0 . ± .
017 (Masiero et al. 2013, 2015) to convert H to D (shown on the right ordinate).As often seen in asteroid families, the small members of the Clarissa family are more dispersed in a p than the largemembers. The envelope of the distribution in ( a p , H ) is consequently “V” shaped (Vokrouhlick´y et al. 2006a). Thesmall family members also concentrate toward the extreme a p values, while there is a lack of asteroids in the familycenter, giving the family the appearance of ears. This feature has been shown to be a consequence of the YORP effect,which produces a perpendicular orientation of spin axis relative to the orbital plane and maximizes the Yarkovskydrift (e.g., Vokrouhlick´y et al. 2006a). Notably, the observed spread of small family members exceeds, by at least afactor of two, the spread due to the escape velocity from (302) Clarissa, and the left side of the family in Fig. 3 isoverpopulated by a factor of ∼ Re-running the hierarchical clustering method (Zappal`a et al. 1990; 1994) on the most recent proper element catalog results in only aslightly larger membership of the Clarissa family. As the orbital structure of the family remains the same, we opt to use the original PDSidentification.
Figure 2.
Orbital structure of the Clarissa family: e p vs. a p (top panel), sine of i p vs. a p (bottom panel). (302) Clarissa ishighlighted by the green diamond. The locations of principal mean motion resonances are indicated by the gray shaded regionswith the three-body J4-S2-1 resonance on the left and the exterior M1/2 resonance with Mars on the right. The resonancewidths were computed using software available from Gallardo (2014) (see also Nesvorn´y & Morbidelli 1998). It can be notedthat the dispersion of the fragments surrounding (302) Clarissa is narrow in e p , but sunward of the J4-S2-1 resonance, fragmentsare more dispersed in e p . The red symbols mark potential interlopers in the family based on their location in the ( a p , H ) plane(see Figure 3). The green ellipses indicate a region in proper element space where the initial fragments would land assuming:(i) isotropic ejection from the parent body, (ii) f (cid:39) ω (cid:39) ◦ (see Appendix A), and (iii) velocity of 20 m s − . This is equal tothe escape velocity from (302) Clarissa ( v esc = 19 . ± . − for D = 43 ± ρ = 1 . − ). The solid gray curves in the bottom panel of Figure 3 delimit the boundaries wherein most family members arelocated. The curves in the figure are calculated from the equation H = 5 log ( | a p − a c | /C ) , (1)where a c is the family center and C is a constant. The best fit to the envelope of the family is obtained with C = (5 ± × − au. As explained in Nesvorn´y et al. (2015), the constant C is an expression of (i) the ejectionvelocity field with velocities inversely proportional to the size, and (ii) the maximum Yarkovsky drift of fragments overthe family age. It is difficult to decouple these two effects without detailed modeling of the overall family structure(Sections 3 and 4). Ignoring (i), we can crudely estimate the Clarissa family age. For that we use t age (cid:39) (cid:18) C − au (cid:19) (cid:16) a . (cid:17) (cid:18) ρ . − (cid:19) (cid:18) . p V (cid:19) / (2)from Nesvorn´y et al. (2015), where ρ is the asteroid bulk density and p V is the visual albedo. For the Clarissa familywe adopt C = (5 ± × − au, and values typical for a C-type asteroid: ρ = 1 . − , p V = 0 .
05, and find t age (cid:39) ±
11 Myr. Using similar arguments Paolicchi et al. (2019) estimated that the Clarissa family is ∼ t age (cid:39)
60 Myr for the Clarissa family, but did not attach an errorbar to this value.
Lowry et al.
Figure 3.
Clarissa family distribution in a p and H (bottom panel). (302) Clarissa is highlighted by the black diamond( D = 43 ± H >
17) are missing in the center of the family and they arepushed toward the borders. The curves are calculated by Eq. (1), where a c = 2 . C = (5 ± × − au (Nesvorn´y et al. 2015). Asteroids located far outside this envelope, shown in red, are suspected interlopers.The top panel shows the a p distribution of family members (suspected interlopers excluded) with 16 . < H <
18 (correspondingto D (cid:39) a p values that would be initially populated by Clarissa family members(assumed ejection speeds (cid:46)
20 m s − ). Objects residing far outside of the curves given by Eq. (1) are considered interlopers (marked red in Figs. 2 and 3)and are not included in the top panel of Fig. 3. Further affirmation that these objects are interlopers could be obtainedfrom spectroscopic studies (e.g., demonstrating that they do not match the family type; Vokrouhlick´y et al. 2006c).The spectroscopic data for the Clarissa family are sparse, however, and do not allow us to confirm interlopers. Wemention asteroid (183911) 2004 CB100 (indicated by a red triangle) which was found to be a spectral type X (likely aninterloper). Asteroid (36286) 2000 EL14 (indicated by a red square in Fig. 3) has a low albedo p V (cid:39) .
06 in the WISEcatalog (Masiero et al. 2011) and Morate et al. (2018) found it to be a spectral type C. Similarly, asteroid (112414)2000 NV42 (indicated by a red circle) was found to be a spectral type C. These bodies may be background objectsalthough the background of primitive bodies in the inner main belt is not large. The absolute magnitudes H couldbe determined particularly badly (according to Pravec et al. (2012) the determination of H may have an uncertaintyaccumulated up to a magnitude), but it is not clear why only these two bodies of spectral type C would have this bias.The absolute magnitude H distribution of Clarissa family members can be approximated by a power law, N ( 75 (Fig. 4). The relatively large value of γ and large size of (302) Clarissarelative to other family members are indicative of a cratering event on (302) Clarissa (Vokrouhlick´y et al. 2006a). Thesignificant flattening of (302) Clarissa in the northern hemisphere (Fig. 1) may be related to the family-forming event(e.g., compaction after a giant cratering event). This is only speculation and we should caution the reader about theuncertainties in shape modeling. In the next section, we describe numerical simulations that were used to explain thepresent orbital structure of the Clarissa family and determine its formation conditions and age. We take particularcare to demonstrate the strength of the J4-S2-1 resonance and its effect on the family structure inside of 2.398 au. Figure 4. Cumulative absolute magnitude H distribution of 169 Clarissa family members (10 suspected interlopers removed;Fig. 3). The large dot indicates (302) Clarissa. The gray reference line is N ( 75. The diameters arelabeled on the top. 2. NUMERICAL MODELAs a first step toward setting up the initial conditions, we integrated the orbit of (302) Clarissa over 1 Myr. Wedetermined the moment in time when the argument of perihelion of (302) Clarissa reached 90 ◦ (note that the argumentis currently near 54 ◦ ). Near that epoch we followed the asteroid along its orbit until the true anomaly reached 90 ◦ (Figure 5). This was done to have an orbital configuration compatible with ω (cid:39) f (cid:39) ◦ (see Appendix A for morecomments) that we used to plot the ellipses in Fig. 2. At that epoch we recorded (302) Clarissa’s heliocentric positionvector R and velocity V , and added a small velocity change δ V to the latter. This represents the initial ejection speedof individual members. From that we determined the orbital elements using the Gauss equations with f = 90 ◦ and ω = 90 ◦ (Zappal`a et al. 1996). Figure 5. Osculating elements of (302) Clarissa over 1 Myr into the future. At time t (cid:39) . 35 ky (red dot), ω = 90 ◦ and f = 90 ◦ . Lowry et al. We generated three distributions with δ V = 10, 20, and 30 m s − to probe the dependence of our results onthis parameter (ejection directions were selected isotropically). Note that δ V = 20 m s − best matches the escapespeed from (302) Clarissa. The assumption of a constant ejection speed of simulated fragments is not a significantapproximation, because we restrict our modeling to D (cid:39) N massive (the sun and planets) and n massless bodies (asteroids in our syntheticClarissa family). For all of our simulations we included eight planets (Mercury to Neptune), thus N = 9, and n = 500test family members. We used a time step of two days and simulated all orbits over the time span of 150 Myr. Thisis comfortably longer than any of the age estimates discussed in Section 1. The integrator eliminates test bodies thatget too far from the sun or impact a massive body. Since (302) Clarissa is located in a dynamically stable zone andthe simulation is not too long, we did not see many particles being eliminated. Only a few particles leaked from theJ4-S2-1 resonance onto planet-crossing orbits. The Clarissa family should thus not be a significant source of near-Earthasteroids.The integrator mentioned above takes into account gravitational forces among all massive bodies and their effect onthe massless bodies. Planetary perturbations of Clarissa family members include the J4-S2-1 and M1/2 resonancesshown in Fig. 3 and the secular perturbations which cause oscillations of osculating orbital elements (Fig. 5). Theprincipal driver of the long-term family evolution, however, is the Yarkovsky effect, which causes semimajor axis driftof family members to larger or smaller values. Thus, we extended the original version of the SWIFT code to allow usto model these thermal accelerations. Details of the code extension can be found in Vokrouhlick´y et al. (2017a) andBottke et al. (2015) (also see Appendix B). Here we restrict ourselves to describing the main features and parametersrelevant to this work.The linear heat diffusion model for spherical bodies is used to determine the thermal accelerations (Vokrouhlick´y1999). We only account for the diurnal component since the seasonal component is smaller and its long-term effect onsemimajor axis vanishes when the obliquity is 0 ◦ or 180 ◦ . We use D = 2 km and assume a bulk density of 1.5 g cm − which should be appropriate for primitive C/B-type asteroids (Nesvorn´y et al. 2015; Scheeres et al. 2015). Consideringthe spectral class and size, we set the surface thermal inertia equal to 250 in the SI units (Delbo et al. 2015). To modelthermal drift in the semimajor axis we also need to know the rotation state of the asteroids: the rotation period P and orientation of the spin vector s . The current rotational states of the Clarissa family members, except for (302)Clarissa itself (see Section 1), are unknown. This introduces another degree of freedom into our model, because wemust adopt some initial distribution for both of these parameters. For the rotation period P we assume a Maxwelliandistribution between 3 and 20 hr with a peak at 6 hr based on Eq. (4) from Pravec et al. (2002). The orientation ofthe spin vectors was initially set to be isotropic but, as we will show, this choice turned out to be a principal obstaclein matching the orbital structure of the Clarissa family (e.g., the excess of members sunward from (302) Clarissa). Wetherefore performed several additional simulations with non-isotropic distributions to test different initial proportionsof prograde and retrograde spins.The final component of the SWIFT extension is modeling the evolution of the asteroids’ rotation state. For this weimplement an efficient symplectic integrator described in Breiter et al. (2005). We introduce ∆ which is the dynamicalelipticity of an asteroid. It is an important parameter since the SWIFT code includes effects of solar gravitationaltorque. We assume that ∆ = ( C − . A + B )) /C, where ( A, B, C ) are the principal moments of the inertia tensor; hasa Gaussian distribution with mean 0.25 and standard deviation 0.05. These values are representative of a populationof small asteroids for which the shape models were obtained (e.g., Vokrouhlick´y & ˇCapek 2002).The YORP effect produces a long-term evolution of the rotation period and direction of the spin vector (Bottke et al.2006; Vokrouhlick´y et al. 2015). To account for that we implemented the model of ˇCapek & Vokrouhlick´y (2004) wherethe YORP effect was evaluated for bodies with various computer-generated shapes (random Gaussian spheroids). Fora 2 km sized Clarissa family member, this model predicts that YORP should double the rotational frequency over amean timescale of (cid:39) 80 Myr.We define one YORP cycle as the evolution from a generic initial rotation state to the asymptotic state with veryfast or very slow rotation, and obliquity near 0 ◦ or 180 ◦ . Given that the previous Clarissa family age estimatesare slightly shorter than the YORP timescale quoted above, we expect that D (cid:39) Summary of Physical and Dynamical Model Assumptions Asteroid Physical Properties (C-type) Value (302) Clarissa’s diameter 43 ± . ± . . − Thermal inertia 250 J m − s − . K − Constant C (see Nesvorn´y et al. 2015) (5 ± × − au Dynamical Properties Value Initial velocity field Isotropic with 10-30 m s − Initial percentage of asteroid retrograde rotation Varying from 50 to 100% by 10% incre-mentsAsteroid rotational period Maxwellian 3-20 hr with peak at 6 hrOnly considered diurnal component of Yarkovsky Drift Dominates over seasonal componentAsteroid dynamical ellipticity ∆ Gaussian with µ = 0 . 25 and σ = 0 . Table 1. Note. The asteroid physical parameters are based on typical values for C-type asteroids (Nesvorn´y et al. 2015). Thedynamical properties stem from current predictions of YORP theory (see Vokrouhlick´y & ˇCapek 2002; ˇCapek & Vokrouhlick´y2004; Breiter et al. 2005; Bottke et al. 2006; Vokrouhlick´y et al. 2015). Lowry et al. ANALYSISWe simulated 500 test bodies over 150 Myr to model the past evolution of the synthetic Clarissa family. For eachbody, we compute the synthetic proper elements with 0.5 Myr cadence and 10 Myr Fourier window (ˇSidlichovsk´y &Nesvorn´y 1996). Our goal is to match the orbital distribution of the real Clarissa family. This is done as follows. Thetop panel of Figure 3 shows the semimajor axis distribution of 114 members of the Clarissa family with the sizes D (cid:39) N i obs , where i spans 25 bins as shown in Figure 3.We use one million trials to randomly select 114 of our synthetic family asteroids and compute their semimajor axisdistribution N i mod for the same bins. For each of these trials we compute a χ -like quantity, χ a ( T ) = (cid:88) i (cid:0) N i mod − N i obs (cid:1) N i obs , (3)where the summation goes over all 25 bins. The normalization factor of χ a ( T ), namely (cid:112) N i obs , is a formal statisticaluncertainty of the population in the i th bin. We set the denominator equal to unity if N i obs = 0 in a given bin.Another distinctive property of the Clarissa family is the distribution of eccentricities sunward of the J4-S2-1 reso-nance. We denote as N +obs the number of members with a p < . 398 au and e p > . N − obs thenumber of members with a p < . 398 au and e p < . D (cid:39) N +obs = 48 and N − obs = 5. It ispeculiar that N +obs /N − obs (cid:39) 10 because the initial family must have had a more even distribution of eccentricities. Thishas something to do with crossing of the J4-S2-1 resonance (see below). As our goal is to simultaneously match thesemimajor axis distribution and N +obs /N − obs , we define χ ( T ) = χ a ( T ) + (cid:0) N +mod − N +obs (cid:1) N +obs + (cid:0) N − mod − N − obs (cid:1) N − obs , (4)where N − mod and N +mod are computed from the model. The Clarissa family age is found by computing the minimum of χ ( T ) (Figure 6). Figure 6. Time-dependence of the χ function defined in Eq. (4): left panel for the YORP1 model, right panel for the YORP2model. These simulations used a 20 m s − initial ejection velocity field that is isotropic in space. Symbols show χ computedby the procedure described in Section 3. The color coding corresponds to models with different fractions of initially retrograderotators (see Figure 7 for a summary of the parameters): 50% (black), 60% (cyan), 70% (green), 80% (red), 90% (blue), and100% (purple). The red and blue parabolas are quadratic fits of χ near the minima of the respective data sets. The dashedlight-gray line marks the value of 27, equal to the number of effective bins. The best age solutions are 54 ± ± χ (cid:39) . 4. In both cases, simulationswith 70%-90% of initially retrograde rotators provide the best solutions. We find that χ ( T ) always reaches a single minimum in 0 < T < 150 Myr. The minimum of χ ( T ) was thendetermined by visual inspection, performing a second-order polynomial fit of the form χ ( T ) (cid:39) aT + bT + c in thevicinity of the minimum, and thus correcting the guessed value to T ∗ = − b/ a . After inspecting the behavior of χ ( T ),we opted to use a ± 15 Myr interval where we fit the second-order polynomial; for instance, between 40 and 70 Myr,if the minimum is at 55 Myr, and so on. A formal uncertainty is found by considering an increment ∆ χ resultingin δT = (cid:112) ∆ χ /a . Thus the age of the Clarissa family is t age = T ∗ ± δT . In a two-parameter model, where the twoparameters are T and the initial fraction of retrograde rotators, we need ∆ χ = 2 . σ ,and ∆ χ = 4 . 61 for a 90% confidence limit or 2 σ (Press et al. 2007, Chapter 15, Section 6). Our error estimatesare approximate. The model has many additional parameters, such as the initial velocity field, thermal inertia, bulkdensity, etc. Additionally, a Kolmogorov-Smirnov two-sample test (Press et al. 2007, Chapter 14, Section 3) wasperformed on selected models (Figures 8 and 9). This test provides an alternative way of looking at the orbitaldistribution of Clarissa family members in the semimajor axis. It has the advantage of being independent of binning. RESULTS4.1. Isotropic velocity field with 20 m s − These reference jobs used the assumption of an initially isotropic velocity field with all fragments launched at 20m s − with respect to (302) Clarissa. This set of simulations included two cases: the (i) YORP1 model (equal likelihoodof acceleration and deceleration of spin), and (ii) the YORP2 model (80% chance of acceleration versus 20% chanceof deceleration). In each case, we simulated six different scenarios with different percentages of initially retrograderotators (from 50% to 100% in increments of 10%). See Figure 7 for a summary diagram of model input parameters.In total, this effort represented 12 simulations, each following 500 test Clarissa family members. Model InputParameters Acceleration:Decelerationof Spin by YORP Initial Ejection Ve-locity of Fragments Fraction of InitialRetrograde Rotators50:50 (YORP1)80:20 (YORP2) 10 m s − 20 m s − 30 m s − Figure 7. Summary of all input parameters for the YORP model simulations as shown in Figure 6 and explained in Section4. This diagram can be read as follows: for example, in the YORP1 model (50:50 chance of acceleration:deceleration of spinby YORP) we simulated fragments with an initial ejection velocity of 20 m s − which included varying fractions of initiallyretrograde rotators from 50% to 100% in increments of 10%. The initial ejection velocity of 20 m s − represents overall 12simulations with the YORP1 and YORP2 models. These simulated cases include only the isotropic velocity field. Figure 6 summarizes the results by reporting the time dependence of χ ( T ) from Eq. (4). In all cases, χ ( T ) reachesa well defined minimum. Initially the test body distribution is very different from the orbital structure of the Clarissafamily and χ (0) is therefore large. For T ≥ 100 Myr, the simulated bodies evolve too far from the center of the family,well beyond the width of the Clarissa family in a p , and χ ( T ) is large again. The minimum of χ ( T ) occurs near 50-60Myr. For the models with equal split of prograde and retrograde rotators (Figure 6, black symbols) the minimum χ ( T ∗ ) (cid:39) 50, which is inadequately large for 27 data points (this applies to both the YORP1 and YORP2 models).This model can therefore be rejected. The main deficiency of this model is that bodies have an equal probability0 Lowry et al. to drift inward or outward in a P (left panel of Figure 10). The model therefore produces a symmetric distributionin semimajor axis, which is not observed (see the top panel of Figure 3). The simulations also show that the M1/2orbital resonance with Mars is not strong enough to produce the observed asymmetry. We thus conclude that the a p distribution asymmetry must be a consequence of the predominance of retrograde rotators in the family. Thisprediction can be tested observationally.The results shown in Figure 6 indicate that the best solutions are obtained when 70%-90% of fragments have initiallyretrograde rotation. These models lead to χ ( T ∗ ) (cid:39) . χ ( T ∗ ) (cid:39) . χ should attain or exceed this level by randomfluctuations is greater than 90% (Press et al. 2007, Chapter 15, Section 2). The inferred age of the Clarissa family is t age = 54 ± t age = 56 ± χ result, we obtain t age = 56 +7 − Myr for the YORP2 model (see Figs.8 and 9). The best fit to the observed a p distribution is shown for YORP2 in Figure 10 (right panel). The modeldistribution for T ∗ = 56 Myr indeed represents an excellent match to the present Clarissa family. Figure 8. Kolmogorov-Smirnov two-sample test: D -value vs. time (top panel), and p -value vs. time (bottom panel), bothcomputed from the cumulative distribution functions F ( a p ) (model) and F ( a p ) (observed). The test was applied to thepreferred YORP2 model (80% preference for acceleration by YORP with 80% of retrograde rotators; Fig. 6). The dashed redline refers to the critical D -value, D α ( α = 0 . t age = 56 +7 − Myr using this test. This result closely matches that obtained with the χ method describedin Section 3. The orbital distribution produced by our preferred YORP2 model is compared with observations in Figure 11. Wenote that the test bodies crossing the J4-S2-1 resonance often have their orbital eccentricity increased. This leadsto the predominance of orbits with e p > . a p < . 398 au. We obtain N +mod = 45 and N − mod = 5, which isnearly identical to the values found in the real family ( N +obs = 48 and N − obs = 5). Suggestively, even the observed sin i p distribution below the J4-S2-1 resonance, which is slightly wider, is well reproduced. We also note a hint of a veryweak mean motion resonance at a P (cid:39) . 404 au which manifests itself as a slight dispersal of e P . Using tools discussedand provided by Gallardo (2014), we tentatively identified it as a three-body resonance J6-S7-1, but we did not proveit by analysis of the associated critical angle.4.2. Anisotropic velocity field with 20 m s − Here we discuss (and rule out) the possibility that the observed asymmetry in a p is related to an anisotropic ejectionvelocity field (rather than to the preference for retrograde rotation as discussed above). To approximately implement1 Figure 9. Kolmogorov-Smirnov two-sample test: cumulative distribution functions F ( a p ) (model, blue line) and F ( a p )(observed, red line). The model distribution is shown for our preferred YORP2 model at T ∗ = 56 Myr. This curve (blue)corresponds to the minimum of the D -value and the maximum p -value plotted in Figure 8. The green and black dashed linesare F ( a p ) ± D α corresponding a 90% confidence band where F ( a p ) − D α ≤ F ( a p ) ≤ F ( a p ) + D α . This interval contains thetrue cumulative distribution function F ( a p ) with a 90% probability. Figure 10. Minimum χ solution for the YORP2 model in both panels with 20 m s − ejection speed with 50% of initiallyretrograde spins (left panel) and 80% of initially retrograde spins (right panel). The solid black line represents the observedfamily corresponding to the top panel of Figure 3. The dashed line is the initial distribution of test bodies. The gray histogramis the model distribution where χ reached a minimum in the simulation corresponding to age solutions of T ∗ = 63 Myr (leftpanel) and 56 Myr (right panel). In the left panel the model distribution is quite symmetric and does not match the observedfamily distribution. There is only a slight asymmetry in the model distribution (left panel) which is due to asteroids leakingout of the family range via the M1/2 resonance. The light-gray bars highlight locations of the principal resonances. an anisotropic velocity field we select test bodies initially populating the left half of the green ellipses in Fig. 2 (i.e.,all fragments assumed to have initial a p lower than (302) Clarissa) and adopt a 50-50% split of prograde/retrograderotators. This model does not work because the evolved distribution in a p becomes roughly symmetrical (with only asmall sunward shift of the center). This happens because the effects of the Yarkovsky drift on the a p distribution aremore important than the initial distribution of fragments in a p .We also tested a model that combined the preference for retrograde rotation with the anisotropic ejection field. Asbefore, we found that the best-fitting models were obtained if there was an ∼ 80% preference for retrograde rotation.The fits were not as good, however, as those obtained for the isotropic ejection field. The minimum χ ( T ∗ ) achieved2 Lowry et al. Figure 11. Orbital evolution of family members in our preferred YORP2 model shown in Fig. 10 (right panel). The properorbital elements e p and a p are shown in the top panel, and sin i p and a p are shown in the bottom panel. The red symbols arethe 114 members of the Clarissa family with sizes D (cid:39) T ∗ = 56 Myr. The gray lines show theevolutionary tracks of the test bodies. was (cid:39) 12, which is significantly higher than the previous result with χ ( T ∗ ) (cid:39) . 4. We therefore conclude that theobserved structure of the Clarissa family can best be explained if fragments were ejected isotropically and there was (cid:39) Isotropic velocity field with 10 and 30 m s − We performed additional simulations with the isotropic ejection field and velocities of 10 and 30 m s − . The maingoal of these simulations was to determine the sensitivity of the results to this parameter. Analysis of the simulationswith 10 m s − revealed results similar to those obtained with 20 m s − (Section 4.1). For example, the best-fittingsolution for the preferred YORP2 had χ ( T ∗ ) (cid:39) . t age = 59 ± − showed that thisvalue is already too large to provide an acceptable fit. The best-fit solution of all investigated models with 30 m s − is χ (cid:39) 25. This happens because the initial spread in the semimajor axis is too large and the Yarkovsky and YORPeffects are not capable of producing Clarissa family ears (see Fig. 3). We conclude that ejection speeds (cid:38) 30 m s − canbe ruled out. Figure 12 shows results similar to Figure 11, but for ejection velocities v ej = 10 m s − (left panel) and v ej = 30 m s − (right panel); all other simulation parameters are the same as in Fig. 11 and the configuration is shownwhen χ of the corresponding run reached a minimum. The former simulation, 10 m s − ejection speed, still providesvery good results. The initial spread in proper eccentricities near (302) Clarissa (at a P (cid:39) . . 404 au suitably extends the family at smaller a P region as the members drift across. This helps to balance the e P distribution of the family members below the J4-S2-1resonance and also provides a tight e P distribution near the M1/2 resonance. On the other hand, the simulation witha 30 m s − ejection speed gives much worse results (the best we could get was χ ( T (cid:63) ) (cid:39) . e P and sin i P is large, and this implies that also the3population of fragments that crossed the J4-S2-1 resonance remains unsuitably large and contradicts the evidence ofa shift toward larger values on its sunward side. Figure 12. Orbital evolution of family members in our preferred YORP2 model (80% of initially retrograde spins and 80%preference for spin acceleration by YORP) for initial ejection velocities of v ej = 10 m s − (left panel) and v ej = 30 m s − (rightpanel). This figure is similar to Figure 11 but with two different ejection speeds. The red symbols are the 114 members of theClarissa family with sizes D (cid:39) χ of the respective simulation, T ∗ = 59 Myr(left panel) and T ∗ = 56 Myr (right panel). The gray lines show the evolutionary tracks of the test bodies in both simulations.5. DISCUSSION AND CONCLUSIONSThe Clarissa family is an interesting case. The family’s location in a dynamically quiet orbital region of the innerbelt allowed us to model its structure in detail. Its estimated age is older than any of the very young families (e.g.,Nesvorn´y et al. 2002; 2006a; 2006b) but younger than any of the families to which the Yarkovsky effect chronologywas previously applied (e.g., Vokrouhlick´y et al. 2006a; 2006b; 2006c). Specifically, we found that the Clarissa familyis 56 ± ρ = 1 . − . In the case of pure Yarkovsky driftthe age scales with ρ as t age ∝ ρ ; higher/lower densities would thus imply older/younger ages. However, in our modelthis scaling is more complicated since altering ρ changes the YORP timescale and the speed of resonance crossing.The initial ejection velocities were constrained to be smaller than (cid:39) 20 m s − , a value comparable to the escapevelocity from (302) Clarissa. We found systematically better results for the model where ∼ 80% of fragments hadrotation accelerated by YORP and the remaining ∼ 20% had rotation decelerated by YORP. This tendency is consistentwith theoretical models of YORP and actual YORP detection, which suggest the same preference (as reviewed inVokrouhlick´y et al. 2015).The most interesting result of this work is the need for asymmetry in the initial rotation direction for small fragments.We estimate that between 70% and 90% of D (cid:39) a < . 406 au (i.e., lower than the semimajor axis of (302) Clarissa) must be retrograderotators today. This prediction can be tested observationally.In fact, prior to running the test cases mentioned in Figure 7 of Section 4 we expected that simulating more retrograderotators in roughly the 80:20 proportion would match the distribution of the observed Clarissa family. We see roughly4