W F Bottke
email: bottke@boulder.swri.edu

D Vokrouhlický

R.-L Ballouz

Olivier S Barnouin

H C Connolly

C Elder

S Marchi

T J Mccoy

Patrick Michel

M C Nolan

B Rizk

D J Scheeres

S R Schwartz

K J Walsh

D S Lauretta

Interpreting the Cratering Histories of Bennu, Ryugu, and Other Spacecraft-explored Asteroids

HAL is

Introduction

The target of NASA's asteroid sample return mission OSIRIS-REx (Origins, Spectral Interpretation, Resource Identification, and Security-Regolith Explorer) is the near-Earth object (NEO) (101995) Bennu. Bennu has a diameter Dast ~ 0.5 km , a 4.4% mean albedo, and a spectral signature consistent with a composition similar to CM-or CI-type carbonaceous chondrite meteorites (e.g., [START_REF] Lauretta | The unexpected surface of asteroid (101955) Bennu[END_REF][START_REF] Hamilton | Evidence for widespread hydrated minerals on asteroid (101955) Bennu[END_REF]. The retrieval and study of primitive asteroidal materials, whose provenance may go back to the earliest times of Solar System history, may allow us to glean insights into the nature of planetesimal and planet formation. Another goal of the OSIRIS-REx mission is to determine whether samples from Bennu can inform us about its individual evolution, as well as that of its parent body. A critical part of this analysis will be to place Bennu's samples into a geologic, geochemical, and dynamical context, and that means learning as much as we can about Bennu's history from its physical and orbital properties. As part of this work, our goal in this paper is to interpret Bennu's cratering history and what it can tell us about Bennu's trek from its formation location, presumably in the main asteroid belt, to its current orbit (e.g., Bottke et al., 2015b).

To set the stage for our work, we first describe what has been inferred about Bennu's collisional and dynamical history to date. A plausible evolution scenario is that Bennu was created in the catastrophic disruption of a main belt parent body with Dast > 100-200 km approximately 1-2 billion years (Ga) ago (e.g., [START_REF] Campins | The origin of asteroid 101955[END_REF][START_REF] Walsh | Introducing the Eulalia and new Polana asteroid families: Re-assessing primitive asteroid families in the inner Main Belt[END_REF]Bottke et al., 2015b).

Using numerical simulations, and building on earlier work by [START_REF] Campins | The origin of asteroid 101955[END_REF] and [START_REF] Walsh | Introducing the Eulalia and new Polana asteroid families: Re-assessing primitive asteroid families in the inner Main Belt[END_REF], Bottke et al. (2015b) argued that Bennu most likely came from the low-albedo Eulalia asteroid family (once called the Polana family) or New Polana asteroid family (the actual family associated with (142) Polana). Both have low inclinations (i ~ 2-3°) and are located in the region adjacent to Jupiter's 3:1 mean motion resonance at ~2.5 au. The largest remnant of the Eulalia family, likely (495) Eulalia, is located at semimajor axis a = 2.487 au, whereas the largest remnant of the New Polana family, (142) Polana, is at a = 2.42 au. The estimated age of the Eulalia family as derived by its dynamical evolution is 830 [+370, -100] Ma, whereas the age of New Polana is thought to be 1400 [+150, -150] Ma, respectively. Using suites of numerical runs, Bottke et al. (2015b) also showed that the New Polana was modestly favored as a source for Bennu over Eulalia by a 70 [+8, -4]% to 30 [+4, -8]% margin, a result consistent with previous work (e.g., [START_REF] Campins | The origin of asteroid 101955[END_REF]).

Bennu's orbit and spin state is affected by the non-gravitational Yarkovsky and Yarkovsky-O'Keefe-Radzievskii-Paddack (YORP) thermal effects (e.g., [START_REF] Rubincam | Radiative spin-up and spin-down of small asteroids[END_REF]Bottke et al., 2006a;[START_REF] Vokrouhlický | The Yarkovsky and YORP effects[END_REF][START_REF] Chesley | Orbit and bulk density of the OSIRIS-REx target Asteroid (101955) Bennu[END_REF][START_REF] Nolan | Detection of Rotational Acceleration of Bennu Using HST Light Curve Observations[END_REF][START_REF] Hergenrother | The operational environment and rotational acceleration of asteroid (101955) Bennu from OSIRIS-REx observations[END_REF]. The former is a small force caused by the absorption of sunlight and re-emission of this energy as infrared photons (heat). The recoil produces a thrust that leads to steady changes in Bennu's semimajor axis over long timescales. The latter is a thermal torque that, complemented by a torque produced by scattered sunlight, modifies Bennu's rotation rate and obliquity. Modeling results indicate that the YORP effect readily modified Bennu's spin axis to a value approaching 180°, the same value it has today, and this allowed the Yarkovsky effect to drive Bennu inward across in the inner main belt (e.g., Bottke et al., 2015b). Additional consequences of the YORP effect on Bennu's shape and surface are discussed below.

After spending most of its lifetime moving inward toward the Sun across the inner main belt, Bennu entered into the ν6 secular resonance that defines the innermost boundary of the main asteroid belt. From there, Bennu was driven onto a high eccentricity (e) orbit where it underwent encounters with the terrestrial planets. One such encounter, most likely with Earth, removed it from the ν6 resonance and placed it onto an a < 2 au orbit. At that point, planetary encounters and smaller planetary resonances moved Bennu onto its current Earth-like orbit with (a, e, i) = (1.126 au, 0.204, 6.035°).

At some point along the way, Bennu achieved an orbit low enough in eccentricity to become collisionally decoupled from the main belt. At that point, sizable collisions on Bennu became far less common, with the NEO population smaller by roughly a factor of 1000 than the main belt (e.g., [START_REF] Bottke | Velocity distribution among colliding asteroids[END_REF]Bottke et al., , 2015a)). Using the population of 682 asteroids with Dast ≥ 50 km defined by [START_REF] Farinella | Collision rates and impact velocities in the main asteroid belt[END_REF], [START_REF] Bottke | Origin of the Spacewatch small Earth-approaching asteroids[END_REF] found that NEOs were largely safe from striking main belt bodies when their aphelion values Q ≤ 1.6 au (e.g., Figs. 2 and 3 of [START_REF] Bottke | Origin of the Spacewatch small Earth-approaching asteroids[END_REF]. According to dynamical runs from Bottke et al. (2015b), we found that the median timescale to go from this boundary to Bennu's current (a, e, i) orbit was 2.6 Ma. Most test bodies took < 20 Ma, though 3% of them managed to avoid it for 70-140 Ma.

Accordingly, if Bennu came from the Eulalia or New Polana families, the time spent on an orbit collisionally decoupled from the main belt was probably a tiny fraction of its entire lifetime.

Therefore, if Bennu's largest craters date back to those times, we can deduce that they were formed by main belt projectiles. For reference, comparable arguments can be made for (162173) Ryugu, the 1-km-diameter carbonaceous chondrite-like target of JAXA's Hayabusa2 sample return mission [START_REF] Watanabe | Hayabusa2 arrives at the carbonaceous asteroid 162173 Ryugu --A spinning top-shaped rubble pile[END_REF], which also likely came from the Eulalia or New Polana families (Bottke et al. 2015b). The net number of impacts, though, may only be part of the story, particularly if Bennu has experienced frequent global crater erasure events.

The origin of Bennu's top-like shape may also tell us about its history. Bennu is a gravitational aggregate made of smaller components, what is often referred to as a "rubble-pile" asteroid (e.g., [START_REF] Barnouin | Shape of (101955) Bennu indicative of a rubble pile with internal stiffness[END_REF][START_REF] Scheeres | The dynamic geophysical environment of (101955) Bennu based on OSIRIS-REx measurements[END_REF]. Michel et al. (2020) argues that Bennu's shape may have been derived from the re-accretion of fragments produced when the parent body disrupted. Alternatively, it may have been spun up by YORP torques into a top-like shape (e.g., [START_REF] Walsh | Formation and evolution of binary asteroids[END_REF].

The YORP effect is also active today. An analysis of rotation data spanning the years 1999-2019 indicate that Bennu is currently spinning up at a rate of (3.63 ± 0.52) × 10 -6 deg day -2 [START_REF] Lauretta | The unexpected surface of asteroid (101955) Bennu[END_REF][START_REF] Nolan | Detection of Rotational Acceleration of Bennu Using HST Light Curve Observations[END_REF][START_REF] Hergenrother | The operational environment and rotational acceleration of asteroid (101955) Bennu from OSIRIS-REx observations[END_REF]. If these kinds of accelerations were common in the past, it seems reasonable that Bennu's shape has been heavily influenced by YORP spin-up processes [START_REF] Scheeres | The dynamic geophysical environment of (101955) Bennu based on OSIRIS-REx measurements[END_REF].

The invocation of the YORP spin-up mechanism to explain the shape of Bennu and other top-shaped asteroids, however, presents us with a paradox. If YORP is actively affecting the shape and surface properties of small asteroids, creating a dynamic environment where landslides, mass shedding events, and satellite formation are common (e.g., [START_REF] Barnouin | Shape of (101955) Bennu indicative of a rubble pile with internal stiffness[END_REF][START_REF] Scheeres | The dynamic geophysical environment of (101955) Bennu based on OSIRIS-REx measurements[END_REF], one would expect to see few if any craters on that surface. Instead, an analysis of images from Bennu indicates that it has several tens of craters of diameters 10 m < Dcrater < 150 m [START_REF] Walsh | Craters, boulders and regolith of (101955) Bennu indicative of an old and dynamic surface[END_REF]. The largest craters are perhaps the most unexpected, because they are likely to be the oldest and the least susceptible to erasure via impact-induced seismic shaking (e.g., [START_REF] Richardson | The global effects of impactinduced seismic activity on fractured asteroid surface morphology[END_REF] or some other process. Comparable crater signatures were also found Ryugu (Sugita et al., 2019). Like Bennu, Ryugu is top-shaped and shows evidence of mass movement. Even small potato-shaped asteroids imaged by spacecraft, such (4179) Toutatis, and (25143) Itokawa, which have mean diameters of ~2.5 and ~0.3 km, show a plethora of craters, with several having diameters Dcrater > 100 m [START_REF] Jiang | Boulders on asteroid Toutatis as observed by Chang'e-2[END_REF][START_REF] Marchi | Cratering on asteroids[END_REF].

When considered together, we are left with only a few options to explain the craters on these small asteroids.

Option 1 is that the surfaces of many small asteroids are in fact ancient. This implies that some process is regulating YORP-driven mass shedding. As discussed in Bottke et al. (2015b), a possible mechanism for this would be "stochastic YORP". [START_REF] Statler | Extreme sensitivity of the YORP effect to small-scale topography[END_REF] showed that modest shape changes to asteroids, produced by a variety of processes (e.g., crater formation, changes to asteroid rotational angular momentum by YORP), caused asteroids' spin rates, but not their obliquities, to undergo a random walk. This mechanism could slow down how often asteroids achieve YORP-driven mass shedding events. In fact, Bottke et al. (2015b) found that some stochastic YORP-like process was needed to explain the orbital distribution of asteroid families such as Eulalia and New Polana. Another possible process with approximately the same effect would be that small asteroids achieve YORP equilibrium states from time to time, where further spin-up or spin-down is minimized until some shape change takes place (e.g., cratering, boulder movement; [START_REF] Golubov | Systematic structure and sinks in the YORP effect[END_REF].

Option 2 would be that Bennu's surface, and the surfaces of Ryugu, Itokawa, and Toutatis, are instead relatively young. The craters found on these worlds would then need to form at a much higher rate than in Option 1. One way to achieve this would be to assume that the crater-projectile scaling laws for small asteroids (hereafter crater scaling laws) allow relatively small impactors to make large craters on the surface of these Dast < 2.5 km bodies. The crater scaling laws for Option 1 would instead predict that larger asteroids are needed to make the observed craters.

At this time, we argue that the crater scaling laws for small asteroids are not well enough constrained to rule out Options 1 or 2 for Bennu. If we treat crater scaling laws as a free parameter, both scenarios appear to be consistent with the observational evidence we have for Bennu thus far, namely that substantial YORP accelerations have been measured [START_REF] Hergenrother | The operational environment and rotational acceleration of asteroid (101955) Bennu from OSIRIS-REx observations[END_REF][START_REF] Nolan | Detection of Rotational Acceleration of Bennu Using HST Light Curve Observations[END_REF], evidence for landslides exist [START_REF] Barnouin | Shape of (101955) Bennu indicative of a rubble pile with internal stiffness[END_REF], yet numerous craters have been identified [START_REF] Walsh | Craters, boulders and regolith of (101955) Bennu indicative of an old and dynamic surface[END_REF]. Comparable arguments can be made for Ryugu, Itokawa, and Toutatis.

What is needed is additional evidence that can tip the balance between Options 1 and 2.

Ultimately, this comes down to finding a way to assess crater scaling laws for asteroids.

Methodology for our Crater Production Models

In this paper, we attempt to glean insights into crater scaling laws for small asteroids like Bennu, Ryugu, Itokawa, and Toutatis by first modeling and interpreting the crater records of large main belt asteroids observed by spacecraft (i.e., diameter Dast > 10 km). Our procedure is to create a crater production model specific to each target asteroid. This involves the calculation of several components:

1. An assessment of the size-frequency distribution (SFD) of the main asteroid belt.

2. A crater scaling law that can transform asteroid impactors from Component 1 into craters on the target asteroid.

3. A calculation of the estimated collision probabilities and impact velocities between objects in the main belt population and the target asteroid.

4. The time that a stable surface on the target asteroid (or possibly the entire target asteroid itself) has been recording craters above a threshold crater diameter. This time will be referred to in the paper as the crater retention age or surface age.

Components 1 through 3 come from models whose accuracy depends on constraints and issues that are discussed in more detail below. Component 4, the crater retention age, is an output value that is calculated from a fit between the observed crater SFD found on the surface of the target asteroid and that target's crater production model (e.g., [START_REF] Marchi | Cratering on asteroids[END_REF].

For each crater production model, we intend to test a range of formulations for Components 1 and 2 against the crater SFD found regionally or globally on the target asteroid. This means that for every target asteroid discussed below, there will be an envelope of model main belt SFDs, possible crater scaling laws, and estimated crater retention ages that provide good fits to the data as measured using chi-squared tests. Our preference is to let these fits tell us which combinations of components yield superior results. At the completion of our runs, a confluence of similar components across many different target asteroids, each with different physical parameters, will allow us to predict those that nature prefers.

We purposely avoid terrains that have reached saturation equilibrium or have experienced substantial crater erasure. This leads us to exclude small craters below some threshold diameter from our analysis, with the definition of "small" defined on a case-by-case basis.

Our method also makes use of a number of assumptions that the reader should understand prior to a more in-depth discussion of the components within each crater production model:

• Assumption 1. The size and shape of main belt SFD has been in steady state for billions of years (within a factor of 2 or so) for projectile sizes that make observable craters on our target asteroids.

As discussed in a review by Bottke et al. (2015a), the main belt is the primary source for the near-Earth asteroid population, which in turn provides impactors to the Moon and other terrestrial planets. The evidence suggests the lunar impact flux over the last 3 to 3.5 Ga has been fairly constant (within a factor of 2 or so) over this time (e.g., [START_REF] Ivanov | The comparison of sizefrequency distributions of impact craters and asteroids and the planetary cratering rate[END_REF][START_REF] Marchi | A New Chronology for the Moon and Mercury[END_REF][START_REF] Hiesinger | How old are young lunar craters[END_REF]; but see also [START_REF] Robbins | New crater calibrations for the lunar crater-age chronology[END_REF][START_REF] Mazrouei | Earth and Moon impact flux increased at the end of the Paleozoic[END_REF]. This constraint suggests that the main belt SFD for asteroids smaller than 10 km or so has largely been in a steady state over this time (within a factor of 2). A strongly decaying main belt SFD would produce a very different lunar impact rate.

Results of collisional evolution models also suggest that a steady state emerged in the main belt SFD over the past several billion years (e.g., Bottke et al., 2015a). Asteroid families are produced from time to time in the main belt, but their fragment SFDs are much smaller than the main belt background SFD, at least for impactor sizes of interest in this paper. Once a family is created, the SFD begins to undergo collisional evolution via the same asteroid disruption laws that affect all other asteroids. This slowly grinds the new family's SFD into the same shape as the background SFD. The consequence is that the main belt is constantly replenished by new breakup events, but these events are rarely substantial enough to strongly modify the overall main belt SFD.

A potential test of Assumption 1 is to compare the crater retention ages of target asteroids in asteroid families (or the surface of a target asteroid that can be connected with the origin of an asteroid family) with independent measures of the family's age. If the surface of the target asteroid in question has been recording impact craters from a time almost immediately after the familyforming event, we would expect all of these age constraints to be similar to one another. Examples of independent chronometers are (i) estimates of asteroid family age from models that track the dynamical evolution of family members and (ii) shock degassing ages of meteorite samples that were reset by impact heating caused by the family-forming event.

A concurrence of ages may represent potential evidence that the components applied in the crater production model are reasonably accurate. We will explore this issue below using data from the asteroids (4) Vesta, (243) Ida, and (951) Gaspra.

• Assumption 2. Most main belt asteroids with diameter Dast > 10 km are on reasonably stable orbits and commonly have been on such orbits for billions of years.

Although the main belt was potentially affected by giant planet migration early in its history (e.g., [START_REF] Morbidelli | The Dynamical Evolution of the Asteroid Belt[END_REF][START_REF] Vokrouhlický | Capture of trans-Neptunian planetesimals in the main asteroid belt[END_REF][START_REF] Nesvorný | Modeling the historical flux of planetary impactors[END_REF], the conclusion of these titanic dynamical events left the majority of Dast > 10 km asteroids on fairly stable orbits within the main belt region. From there, new Dast > 10 km asteroids are created from time to time by family-forming events, but they are unlikely to move far from the orbits on which they were placed by the ejection event itself (e.g., [START_REF] Nesvorný | Identification and dynamical properties of asteroid families[END_REF].

Evidence for this comes in a variety of forms, ranging from the calculations of asteroid proper elements for larger main belt asteroids, where dynamical stability can be demonstrated (e.g., [START_REF] Knežević | The Determination of Asteroid Proper Elements[END_REF], to billion-year integrations of the future dynamical evolution of Dast > 10 km asteroids, where only a small fraction can escape the main belt [START_REF] Nesvorný | Dynamical origin and terrestrial impact flux of large near-Earth asteroids[END_REF]. All main belt asteroids undergo modest oscillations in their eccentricities and inclinations from secular perturbations, but the forced components of this oscillation do not modify the free components.

These results indicate that nearly all of our large target asteroids have been in the same approximate orbits for a long enough period of time that our crater production models can be based on their present-day orbits.

• Assumption 3. Models are currently the best option to estimate the main belt SFD at sub-km sizes.

As discussed in more detail below, existing asteroid surveys are unable to detect large numbers of sub-km main belt asteroids, and those sub-km bodies that have been detected have to be carefully debiased to avoid selection effects (i.e., for a given absolute magnitude, a survey will find more high-albedo bodies than low-albedo bodies; [START_REF] Morbidelli | The shallow magnitude distribution of asteroid families[END_REF]Maseiro et al. 2011).

To sidestep this limitation, we will use model main belt SFDs calculated from collision evolution models as input for our crater production models. These model main belt SFDs are constructed to fit existing main belt constraints (as we understand them) and therefore are probably the best we can do with what is available at this time (e.g., Bottke et al. 2015a).

In the next few sections, we discuss our calculations of the components discussed above, starting with Component 1, the predicted main belt SFD.

Deriving a Model Main Belt Size Frequency Distribution (Component 1)

Understanding Collisional Evolution in the Main Belt

To understand cratering on Bennu and other main belt asteroids, our first task is to assess the main belt SFD (i.e., Component 1 from Sec. 1.2). This entails modeling how the main asteroid belt undergoes collisional evolution.

First, although the main belt has a diverse population, nearly all asteroids have orbits that cross one another, especially when secular perturbations are included [START_REF] Bottke | Velocity distribution among colliding asteroids[END_REF][START_REF] Bottke | Origin of the Spacewatch small Earth-approaching asteroids[END_REF].

For example, using the 682 asteroids with Dast > 50 km located between 2 and 3.2 au (Bottke et al., 1994), we find 90% and 71% of individual asteroids cross 80% and 90% of the population, respectively. Even those located along the innermost edge of the main belt near 2.2 au can still be struck by nearly half of all main belt asteroids. Effectively, this means there are no hiding places.

Accordingly, one would expect the shape of the impactor size frequency distribution (SFD) hitting most target bodies should largely represent an amalgam of the main belt SFD as a whole.

Second, collisional evolution models indicate that the main belt SFD is in a quasi-steady state with a wave-like shape driven by the shape of the asteroid disruption law (e.g., Bottke et al., 2005aBottke et al., ,b, 2015a)). Assuming all asteroids disrupt in a similar manner, which impact modeling work suggests is a fairly reasonable approximation (e.g., [START_REF] Jutzi | The structure of the asteroid 4 Vesta as revealed by models of planet-scale collisions[END_REF], simulations that produce the best match with both the main belt SFD and constraints provided by asteroid families indicate that asteroid disruption scaling laws undergo a transition between strength and gravity-scaling near D ~ 0.2 km [START_REF] Benz | Catastrophic disruptions revisited[END_REF]Bottke et al., 2005aBottke et al., ,b, 2015a)). Asteroids near this transition are relatively easy to disrupt, leading to a relative deficit of bodies with Dast ~ 0.1-0.5 km. This "valley" in the SFD leads to an overabundance, or a "peak", of multi-kilometer bodies that would be destroyed by such projectiles. Collisional models suggest this peak in the main belt SFD is near Dast ~2-3 km (e.g., O'Brien and [START_REF] O'brien | Steady-state size distributions for collisional populations: Analytical solution with size dependent strength[END_REF]reviewed in Bottke et al., 2015a).

As new-fragment SFDs are input into the asteroid belt from cratering or catastrophic disruption events, the individual bodies in the SFD undergo collisional evolution. As this grinding proceeds, for asteroids with Dast < 10 km, the shapes of the new-fragment SFDs take on the same wavy profile as the background main belt SFD over tens to hundreds of million years (e.g., Bottke et al. 2005aBottke et al. ,b, 2015a)). In this manner, the wavy shape of main belt SFD can considered to be in a quasi-steady-state.

Collisions may not be the only mechanism affecting asteroid sizes and the wavy shape of the main belt SFD. Asteroids with diameters smaller than a few km may also be affected by mass shedding events produced by YORP thermal torques, the same processes that can modify the spin vectors of small asteroids (e.g., [START_REF] Marzari | Combined effect of YORP and collisions on the rotation rate of small Main Belt asteroids[END_REF][START_REF] Jacobson | Effect of rotational disruption on the size frequency distribution of the main belt asteroid population[END_REF]. The influence of YORP torques on asteroid sizes and the main belt SFD itself depends on the frequency of these mass shedding events (e.g., Bottke et al. 2015a). Any changes to the main belt SFD produced by YORP mass shedding, however, would drive new-fragment SFDs to the same shape as the background main belt SFD.

The consequence is that the main belt SFD likely maintains a wavy profile that stays relatively constant over billions of years. The absolute number of asteroids in the inner, central, and outer main belt SFDs may change as asteroids are dynamically lost or as new families are formed, but modeling work suggests that these effects rarely modify the overall shape of the main belt SFD as a combined whole for very long.

Motivation for Generating a Different Main Belt Size Distribution

With that said, there are several reasons to consider different formulations of the main belt SFD than those discussed in Bottke et al. (2005b). The changes we suggest below have been driven by substantial progress in small body studies over the last two decades. In that time, a plethora of new data has been obtained on the shape of the present-day NEO SFD from a wide variety of surveys (e.g., the Catalina Sky Survey, Lincoln Near-Earth Asteroid Research (LINEAR), Panoramic Survey Telescope and Rapid Response System (Pan-STARRS), Spacewatch, the Near-Earth Object Wide-field Infrared Survey Explorer (NEOWISE) survey) and on asteroid crater SFDs from various missions (see Secs 3 and 4). Both components -revised asteroid SFDs and new asteroid crater SFDs -suggest the SFD presented in Bottke et al. (2005a,b) may be modestly inaccurate for Dast < 1 km.

For example, consider Figure 2 of Bottke et al. (2015a). It shows the cumulative model main belt and NEO SFDs of Bottke et al. (2005b) against recent formulations of the NEO SFD by [START_REF] Harris | The population of near-Earth asteroids[END_REF] (see also [START_REF] Stokes | Update to Determine the Feasibility of Enhancing the Search and Characterization of NEOs. Report of the Near-Earth Object Science Definition Team[END_REF]. We find the shape of [START_REF] Harris | The population of near-Earth asteroids[END_REF] NEO SFD is fairly wavy, with substantial slope changes taking place near Dast ~ 0.1-0.2 km and 2-3 km. Bottke et al. (2005b) instead predicted that (i) the smaller of the two inflection points should occur at Dast ~ 0.5 km, (ii) that a less shallow slope should occur between 0.1-0.2 < Dast < 2-3 km, and (iii) a less steep SFD should occur between 0.01 < Dast < 0.1 km.

While features (ii) and (iii) are somewhat dependent on the removal rates of small asteroids from the main belt via the Yarkovsky effect, feature (i) cannot be explained in such a manner. In general, an inflection point in the SFD of a source population should also be reflected in the daughter population unless removal rates are highly variable.

From the crater perspective, Figure 4 of [START_REF] Marchi | Cratering on asteroids[END_REF] showed a fit between the Bottke et al. (2005b) formulation of the main belt SFD and the SFD of 0.1 < Dcrater < 10 km craters found on or near Vesta's Rheasilvia basin. It indicated that the [START_REF] Marchi | Cratering on asteroids[END_REF] fit, while tolerable, seemed to miss a key feature and inflection point between 0.7 < Dcrater < 2 km. A mismatch in this size range would be consistent with Bottke et al. (2005b) predicting a key inflection point is at Dast ~ 0.5 km rather than 0.1-0.2 km.

There are several plausible ways we could modify the main belt SFD of Bottke et al. (2005b) in the size range of interest:

1. Modifying the Yarkovsky depletion rates of asteroids from the main belt.

2. Allowing YORP-driven mass shedding to strongly affect the diameters of sub-km asteroids, as suggested by [START_REF] Marzari | Combined effect of YORP and collisions on the rotation rate of small Main Belt asteroids[END_REF] and [START_REF] Jacobson | Effect of rotational disruption on the size frequency distribution of the main belt asteroid population[END_REF].

3. Modify the disruption scaling law for main belt asteroids.

We do not favor Scenario 1. Our tests using the Bottke et al. (2005b) model indicate that to move the position of a main belt inflection point from 0.5 to 0.1-0.2 km, we would need to assume (i) much larger Yarkovsky-driven removal rates than in Bottke et al. (2005b), which would require even more main belt disruptions to keep the NEO population resupplied, and (ii) the removal process has a strong size dependence between 0.1-0.2 and 0.5 km. Such changes produce strong modifications to the model NEO SFD, giving it a shape inconsistent with the observed SFD.

We find Scenario 2 to be more intriguing, with modeling work from [START_REF] Marzari | Combined effect of YORP and collisions on the rotation rate of small Main Belt asteroids[END_REF] and [START_REF] Jacobson | Effect of rotational disruption on the size frequency distribution of the main belt asteroid population[END_REF] suggesting YORP-driven mass shedding could be a major factor in decreasing the diameter of sub-km bodies and thereby changing the main belt SFD. A potential concern with this hypothesis, however, is that small asteroids observed by spacecraft have a number of Dcrater > 0.1 km craters (e.g., Bennu, Ryugu, Itokawa, Toutatis). If the Option 1 interpretation turns out to be true, and these asteroids have long crater retention ages for the largest craters, it would rule out substantial YORP-driven mass shedding from these worlds. Note that this does not mean that YORP is unimportant; it still provides an easy way to explain the obliquities, top-like shapes, the existence of satellites, and the mass shedding events seen for many small asteroids (e.g., [START_REF] Jewitt | The active asteroids[END_REF]. Nevertheless, it would imply that YORP's ability to influence the main belt SFD may be more limited than suggested by these models.

In this paper, we focus our investigation on Scenario 3. Our work indicates that it is possible to modify the asteroid disruption scaling law in a manner that yields a main belt SFD consistent with constraints (e.g., shape of the observed main belt SFD, number of asteroid families, asteroid craters, NEO SFD, laboratory impact experiments). With that said, though, Scenario 2 might still be a major player in explaining the shape of the main belt SFD in this size range.

Modeling Collisional Evolution in the Main Asteroid Belt

Collisional and dynamical depletion evolution code (CoDDEM)

Most of the asteroids that hit Bennu-sized bodies are a few tens of meters or smaller in diameter, well below the observational limit of the asteroid belt. For reference, current surveys only able detect large numbers of ~1-2 km diameter bodies (e.g., [START_REF] Jedicke | Observational selection effects in asteroid surveys[END_REF][START_REF] Gladman | On the asteroid belt's orbital and size distribution[END_REF]. Accordingly, the precise nature of the impactor population making craters on most asteroids observed by spacecraft is not yet known. Progress is being made, with digital tracking on ground-based telescopes having great potential (e.g., [START_REF] Heinze | The Flux Distribution and Sky Density of 25th Magnitude Main Belt Asteroids[END_REF]. New data may also become available within the 2020's from both the Large Synoptic Survey Telescope and spacebased infrared surveys like the Near-Earth Object Surveillance Mission (formally NEOCam). Still, a full observational assessment of the sub-km main belt SFD will not be available for some time.

Until that time arrives, it makes sense to use collisional evolution models to estimate the unknown nature of the small body main belt SFD. To this end, we model the main belt SFD using the self-consistent one-dimensional collisional evolution code CoDDEM (Collisional and Dynamical Depletion Evolution Model). Model details and the testing procedure for CoDDEM are discussed in Bottke et al. (2005a,b; see also the review in Bottke et al., 2015a). Here we provide the essentials needed to understand our new results.

We run CoDDEM by entering an initial main belt SFD where the population (N) has been binned between 0.0001 km < D < 1000 km in logarithmic intervals dLogD = 0.1. The particles in the bins are assumed to be spherical and are set to a bulk density of 2.7 g cm -3 , a common asteroid bulk density value. CoDDEM then computes the time rate of change in the differential population N per unit volume of space over a size range between diameter D and D + dD [START_REF] Dohnanyi | Collisional model of asteroids and their debris[END_REF][START_REF] Williams | Size distribution of collisionally-evolved asteroidal populations. Analytical solution for self-similar collision cascades[END_REF]:

∂N ∂t (𝐷, 𝑡) = -𝐼 !"#$%&' + 𝐼 ($)* -𝐼 !+, (1)
Here IDISRUPT is the net number of bodies that leave between D and D + dD per unit time from catastrophic disruptions. The collisional lifetime of a given target body in a bin in the current main belt is computed using estimates of the intrinsic collision probability and mean velocities between asteroids in the main belt, defined as Pi = 2.86 × 10 -18 km -2 yr -1 and Vimp = 5.3 km s -1 , respectively [START_REF] Bottke | Velocity distribution among colliding asteroids[END_REF]; see also Bottke et al., 2015a).

The projectile capable of disrupting Dtarget is defined as ddisrupt:

𝑑 -./0123 = / 2𝑄 ! * 𝑉 .52 6 3 7/9 𝐷 3:0;<3 (2)
We set 𝑄 ! * as the critical impact specific energy-the energy per unit target mass needed to disrupt the target and send 50% of its mass away at escape velocity. Our functions for 𝑄 ! * at different asteroid sizes is tested below.

When a body breaks up, the results go into the IFRAG parameter, which describes the number of bodies entering a given size bin per unit time that were produced by a given disruption event.

CoDDEM uses fragment SFDs as discussed in Bottke et al. (2005a,b). The implication is that the destruction of large asteroids serves as a source to replenish the small body population via a "collisional cascade".

The IDYN parameter account for the number of bodies lost from a given size bin via dynamical processes, such as asteroids being removed by planetary perturbations or an object entering into a dynamical resonance via the Yarkovsky effect and escaping into planet-crossing orbits. This component is used to create our synthetic NEO SFD from the main belt population, as described in Bottke et al. (2005b).

In those runs, which we exactly duplicate in our new simulations, it was assumed that the primordial main belt contained on the order of 200 times the number of objects in the existing main belt, with the vast majority of the material ejected by interactions with planetary embryos within 1-2 Myr of the formation of the first solids. The dynamical removal mechanism used in Bottke et al. (2005b) may or may not end up reflecting reality, but that is not the salient point.

Their model results instead serve as a reasonable proxy for scenarios where a large population of small bodies on planet-crossing orbits early in Solar System history batters the surviving main belt population. This may include the removal of primordial main belt asteroids onto planet-crossing orbits via interactions with migrating giant planets (e.g., [START_REF] Walsh | A low mass for Mars from Jupiter's early gas-driven migration[END_REF] or early giant planet instabilities (e.g., [START_REF] Clement | Mars' growth stunted by an early giant planet instability[END_REF]Nesvorný et al., 2018). It is even possible the Bottke et al. (2005b) model results are fairly consistent with a primordial low-mass asteroid belt bombarded by populations introduced into the terrestrial planet region by planet formation processes. In terms of our model results, all these small body sources provide an additional source of early collisional evolution that sets the stage to explain the current main belt SFD.

Initial conditions and model constraints for CoDDEM

The initial main belt population entered into CoDDEM is divided into two components that are tracked simultaneously: a small component of main belt asteroids that will survive the dynamical excitation event (Nrem) and a much larger component that will be excited and ejected from the main belt (Ndep). Thus, our initial population is N = Nrem + Ndep. We can use this procedure because we know in advance the dynamical fate of each population via the dynamics simulations described in Bottke et al. (2015b). The two populations undergo comminution with themselves and with each other. When Ndep = 0, CoDDEM tracks the collisional and dynamical evolution of Nrem alone for the remaining simulation time.

The size and shape of our initial size distribution was determined by running different initial populations through CoDDEM-like codes, then testing the results against the constraints described in Sec. 4 (Bottke et al., 2005a,b). The size distribution that provided the best fit for Nrem followed the observed main belt for bodies with Dast > 200 km, an incremental power law index of -4.5 for bodies with 110 < Dast < 200 km, and an incremental power law index of -1.2 for bodies with Dast < 110 km (Bottke et al., 2005b). The initial shape of the Ndep population is always the same as Nrem, and its size is set to Ndep = 200 Nrem. Additional starting condition details can be found in Bottke et al. (2005a,b).

For constraints in Bottke et al. (2005b), our model main belt SFD at the end of 4.6 Ga of evolution had to reproduce the wavy-shaped main belt SFD for Dast > 1 km. To determine its value, we converted the absolute magnitude H distribution of the main belt described by [START_REF] Jedicke | Observational selection effects in asteroid surveys[END_REF], who combined observations of bright main belt asteroids with renormalized results taken from the Sloan Digital Sky Survey [START_REF] Ivezić | Solar system objects observed in the Sloan Digital Sky Survey commissioning data[END_REF], into a size distribution. This was accomplished using the relationship (Fowler and Chilemi, 1992;Appendix in Pravec and Harris, 2007):

𝐷 :/3 (𝑘𝑚) = 1329 × 10 =>/? 𝑝 @ =7/6
(3) and a representative visual geometric albedo 𝑝 @ = 0.092. The shape of the main belt SFD is shown as the large dots in Fig. 1. This SFD is in general agreement with the diameter-limited survey produced by WISE [START_REF] Masiero | Main belt asteroids with WISE/NEOWISE I: Preliminary albedos and diameters[END_REF], though their study is only complete in the outer main belt to asteroids larger than Dast > 5 km. It also has had some success matching crater SFDs on asteroids (e.g., [START_REF] Marchi | Cratering on asteroids[END_REF], though we will return to this issue below.

INCLUDE FIGURE 1 HERE

There have been many additional attempts to estimate the shape of the main belt SFD since Bottke et al. (2005a,b). We only mention a few of these examples here. [START_REF] Gladman | On the asteroid belt's orbital and size distribution[END_REF] used a pencil beam survey of main belt asteroids and their likely colors to generate their SFD. Test fits of their SFD against crater SFDs on Vesta, however, have not been as successful as those derived from Bottke et al. (2005a,b) (Marchi et al., 2012a). [START_REF] Ryan | The kilometer-sized Main Belt asteroid population revealed by Spitzer[END_REF] used the Spitzer space telescope to target known objects, find their diameters, and eventually generate a main belt SFD.

Their results are similar to [START_REF] Gladman | On the asteroid belt's orbital and size distribution[END_REF] in many respects. Accordingly, their SFD would likely also have similar problems matching constraints.

A potential issue with results from [START_REF] Gladman | On the asteroid belt's orbital and size distribution[END_REF] and [START_REF] Ryan | The kilometer-sized Main Belt asteroid population revealed by Spitzer[END_REF] is how their methods treat observational selection effects near the detection limit of main belt surveys. In an absolute magnitude-limited survey, it is easier to detect high-albedo S-type asteroids than lowalbedo C-type asteroids. This bias is pervasive through the catalog of known main belt objects.

Studies employing this catalog may be overemphasizing S-types at the expense of C-types, which are numerous in the outer main belt. This effect was demonstrated by [START_REF] Masiero | Main belt asteroids with WISE/NEOWISE I: Preliminary albedos and diameters[END_REF], who

showed that nearly all of the main belt asteroids discovered by WISE were low-albedo. Up to that point, these bodies had been missed by telescopes looking in visual wavelengths.

Finally, a formulation of the main belt SFD by [START_REF] Minton | Re-examining the main asteroid belt as the primary source of ancient lunar craters[END_REF] indicated that it should change slope close to Dast ~ 3.5 km. We point out that this break is discordant with the shape of the inner main belt SFD determined by the diameter-limited WISE survey, which shows no change in slope at that size [START_REF] Masiero | Main belt asteroids with WISE/NEOWISE I: Preliminary albedos and diameters[END_REF]. We also see no change in slope at Dast ~ 3.5 km in the observed NEO SFD [START_REF] Harris | The population of near-Earth asteroids[END_REF][START_REF] Stokes | Update to Determine the Feasibility of Enhancing the Search and Characterization of NEOs. Report of the Near-Earth Object Science Definition Team[END_REF] (Fig. 1); recall that the main belt is the primary source for NEOs, so a change in slope in the parent size distribution should probably be seen in the daughter size distribution as well.

A second set of constraints for Bottke et al. (2005a,b) was provided by asteroids families, particularly those that are potentially too large to be dispersed by the Yarkovsky effect over the age of the Solar System. Using hydrocode simulations from [START_REF] Durda | Size-frequency distributions of fragments from SPH/ N-body simulations of asteroid impacts: Comparison with observed asteroid families[END_REF] to estimate the amount of material in families located below the observational detection limit, Bottke et al. (2005a, b) suggested that ∼ 20 families have been produced by the breakup of Dast > 100 km asteroids over the past ∼ 3.5 Ga. Although there have been recent attempts to update this number (e.g., see the review of this issue in Bottke et al., 2015a), we believe the distribution used by Bottke et al. (2005b) is still reasonable. Here we adopt the same constraint; we assume that the size distribution bins centered on Dast = 123.5,155.5,195.7,246.4,310.2,and 390.5 km experienced 5, 5, 5, 1, 1, and 1 breakups over the past 3.5 Ga, respectively. Our testing procedure also gives us some margin, so assuming that additional large asteroids disrupted over the past 3.5 Ga can be considered reasonable as well.

To quantify the fit between the model and observed population, we follow the methods described in Bottke et al. (2005a,b). Our first metric compares the shape of the model main belt SFD to a small envelope of values surrounding the observed main belt SFD (defined as NMB):

𝜓 #(! 6 = > / 𝑁 $AB (𝐷) -𝑁 BC (𝐷) 0.2 𝑁 BC (𝐷) 3 6 ! (4)
We assume that our model is a good fit if lies within 20% of the observed main belt between 0.98 km and 390.5 km (across 27 incremental bins) as defined by Bottke et al. (2005a) (see also [START_REF] Jedicke | Observational selection effects in asteroid surveys[END_REF]. As discussed in Bottke et al. (2005a) , is better than 2σ (i.e., probability >5%).

𝑄 ! * (𝑅) = 𝑎𝑅 F + 𝑏𝑅 G (5) 𝑎 = 𝑄 ! * D)C 𝑅 D)C F 1 1 - 𝛼 𝛽 H 𝑅 D)C 𝑅 5.E I G=F (6) 𝑏 = - 𝛼 𝛽 𝑎𝑅 5.E F=G (7)
Here R = D / 2 and 𝑅 D)C = 𝐷 D)C / 2. The parameters for our different 𝑄 ! * functions, and their rate of success against our main belt testing metrics in Sec. 2.2, are given in Table 1.

INSERT TABLE 1 HERE

We show our model SFDs #1-8 in Fig. 1. These test runs, corresponding to 𝑄 ! * #1-8 (Fig The power law slopes of the SFDs for Dast < 0.1 km in Fig. 1 range from q = -2.6 cumulative for SFD #1 to q = -2.7 for SFD #8. These outcomes match predictions from O'Brien and [START_REF] O'brien | Steady-state size distributions for collisional populations: Analytical solution with size dependent strength[END_REF], who show that the slope of the 𝑄 ! * function in the strength regime, defined using the α parameter in Table 1, yield these approximate values for the α range shown there (i.e., -0.35 to -0.63). Our results also match observational constraints of the main belt SFD from [START_REF] Heinze | The Flux Distribution and Sky Density of 25th Magnitude Main Belt Asteroids[END_REF], who used Dark Energy Camera observations of main belt asteroids and digital tracking methods to find a slope of -2.575 < q < -2.825.

The cumulative power law slope between the inflection points in Fig. 1, located between Dast ~ 0.2-0.5 km and 2-3 km, is shallower than the q values above. If we measure the slope for all of our model SFDs between 0.5 km and 1.5 km, we find values that go from q = -1.5 for SFD #1 to q = -1.2 for SFD #8. [START_REF] Heinze | The Flux Distribution and Sky Density of 25th Magnitude Main Belt Asteroids[END_REF] report a cumulative slope in this range of q = -1.31 ± 0.01, a value that matches Yoshida et al. ( 2007) (q = -1.29 ± 0.02) but disagrees with Yoshida et al. ( 2003) (q = -1.2). If we assume the preferred slope in this part of the main belt SFD is indeed q = -1.3, the best match comes from SFD #6, with q = -1.3.

The intriguing matches between our SFDs and observational data are necessary but not sufficient proof that our 𝑄 ! * functions reflect reality. For example, [START_REF] Holsapple | The catastrophic disruptions of asteroids: History, features, new constraints and interpretations[END_REF] point out that asteroid disruption scaling laws with α parameters more negative than -0.5 are inconsistent with those inferred from materials tested to date. They instead argue that slopes in the strength regime of -0.2 to -0.3 provide the best matches with scaling law theory. Taken at face value, the best match to α parameters of -0.2 to -0.3 comes from our baseline 𝑄 ! * #1, which yields model SFD #1 within CoDDEM. As we will show below, however, this SFD does not reproduce crater SFDs on many different asteroids as well as other choices.

There are different ways to potentially resolve this paradox beyond simply assuming that our collisional evolution model is inaccurate. The first possibility would be that YORP-spin up is indeed a major factor in the disruption of small asteroids (i.e., Scenario 2 from Sec. 2.2), and that the demolition of small asteroids from this effect is needed in combination with 𝑄 ! * #1 to get the correct SFD. In other words, our steeper α parameter is compensating for the lack of YORP disruption in our model.

A second possibility is that existing impact studies have not yet accounted for the unusual material properties found on some small asteroids. For example, in Hayabusa2's Small Carry-on Impactor (SCI) experiment, a 2.-kg copper plate was shot into the surface of the km-sized carbonaceous chondrite-like asteroid Ryugu at 2 km/s, where it made a semi-circular crater with a rim to rim diameter of 17.6 ± 0.7 m [START_REF] Arakawa | An artificial impact on the asteroid 162173 Ryugu formed a crater in the gravity-dominated regime[END_REF]. This outcome was a surprise to many impact modelers, in that Ryugu's surface acted like it had the same strength as cohesionless sand upon impact (i.e., the crater formed in the gravity-dominated regime). Related studies suggest that the boulders on Ryuyu have estimated porosities as large as 55% [START_REF] Grott | Low thermal conductivity boulder with high porosity identified on C-type asteroid (162173) Ryugu[END_REF]. Put together, these results may indicate that modified asteroid disruption laws are needed to accommodate how carbonaceous chondrite-like asteroids with Ryugu-like properties behave in a disruption event.

The stage is now set to test our eight bounding model asteroid SFDs (Fig. 1) against observed crater SFDs on asteroids observed by spacecraft.

Crater Scaling Laws for Asteroids (Component 2)

To determine the crater retention age of a given asteroid surface, we need to know the crater scaling law that turns projectiles into craters. Typical crater scaling laws require a range of projectile quantities (e.g., size, mass, impact velocity, impact angle, composition, internal structure) and target quantities (e.g., target gravity, surface composition, structure and density, target interior structure and density, effects of surface and internal porosity). Unfortunately, many of these quantities are unknown for observed asteroids. Our ability to calibrate crater scaling laws is also somewhat limited, given that most test data come from laboratory shot experiments, conventional explosions, or nuclear bomb detonations. The energies involved in making observed asteroid craters is typically orders of magnitude higher than the energies used to generate our crater scaling law constraints, even those from nuclear blasts.

Holsapple and Housen crater scaling law

A common crater scaling law used in asteroid studies is similar to the Holsapple and Housen (2007) formulation of the Pi-group scaling law (e.g., used by [START_REF] Marchi | Cratering on asteroids[END_REF]; see also [START_REF] Tatsumi | Cratering efficiency on coarse-grain targets: Implications for the dynamical evolution of asteroid 25143 Itokawa[END_REF]:

𝐷 3 = 𝑘𝑑 K 𝑔𝑑 2𝑉 & 6 H 𝜌 𝛿 I 6H/I + / 𝑌 𝜌𝑉 & 6 3 (6KI)/6 H 𝜌 𝛿 I H(6KI)/I P =I/(6KI) (8)
Here the transient crater diameter, defined by Dt, can be found using the impactor properties (impactor diameter d, velocity perpendicular to the surface Vp, bulk density δ) together with the target properties (density of target material ρ, strength of target material Y, surface gravity g).

For planets and large asteroids, the input of surface gravity g into such crater scaling law equations is straightforward; the combination of their largely spherical shapes and relatively slow spin rates means accelerations across their surfaces are similar. For smaller asteroids, however, the calculation of an effective surface gravity can be complicated by irregular shapes and centrifugal forces.

As an example, consider (243) Ida. Its elongated shape ( 59.8 × 25.4 × 18.6 km) and rapid spin period (P = 4.63 h) leads to a wide range of surface accelerations (0.3 to 1.1 cm s -2 ) [START_REF] Thomas | The shape of Ida[END_REF]. Therefore, when applying this scaling law to Ida, we follow the lead of [START_REF] Schmedemann | The cratering record, chronology and surface ages of (4) Vesta in comparison to smaller asteroids and the ages of HED meteorites[END_REF] and choose a single representative g value from this range (i.e., 0.7 cm s -2 ) as input into our crater scaling laws. We follow suit for the other target asteroids in this paper, whose g values, along with a corresponding reference, are given in Table 2.

PLACE TABLE 2 HERE

Additional parameters (k, ν, µ) account for the nature of the target terrain (i.e., whether it is hard rock, cohesive soil, or porous material). Common parameters for hard rocks are k = 0.93, ν = 0.4, µ = 0.55, and for cohesive soils are k = 1.03, ν = 0.4, µ = 0.41 (e.g., Marchi et al., 2011[START_REF] Marchi | Cratering on asteroids[END_REF].

The yield strength Y of different asteroid target materials is unknown, but we can bracket possibilities using reference values, which range from lunar regolith (Y = 1 × 10 5 dynes cm -2 ; 3 × 10 5 dynes cm -2 at 1 m depth) to dry soil (Y = 3 × 10 6 dynes cm -2 ) to dry desert alluvium (Y = 7 × 10 5 dynes cm -2 ) to soft dry rock/hard soils (Y = 1.3 × 10 7 dynes cm -2 ) to hard rocks and cold ice (Y = 1.5 × 10 8 dynes cm -2 ) (Holsapple and Housen, 2007). In general, when yield strength increases, the craters formed from projectiles are smaller, which translates into an older surface for a given crater SFD.

We also account for the collapse of the transient crater, such that the final crater size is Dcrater = λ Dt. The value λ is ≥ 1 and it is usually determined empirically. A common value for λ is 1.2, but smaller and larger values can also be found in the literature.

In terms of the final crater size, our tests show that larger values for λ can be counteracted by increasing Y; the two trade off of one another. To keep things simple when asteroid parameters are largely unknown, and to limit the amount of interpretation needed for our results, we decided to apply λ = 1.2 and vary Y for our results. Hence, we assume that

𝐷 M0:3<0 = 1.2 𝐷 3 (9)

Ivanov crater scaling law

Another commonly used asteroid scaling law, reformulated from [START_REF] Schmidt | Some recent advances in the scaling of impact and explosion cratering[END_REF], comes from [START_REF] Ivanov | Size-frequency distributions of planetary impact craters and asteroids[END_REF] (see corrected version in [START_REF] Schmedemann | The cratering record, chronology and surface ages of (4) Vesta in comparison to smaller asteroids and the ages of HED meteorites[END_REF]. It has the form:

𝐷 3 𝑑(𝛿 𝜌) ⁄ N.P9 (𝑉 2 sin 𝛼) N.?? = 1.21 [(𝐷 #* + 𝐷 3 )𝑔] N.6Q (10)
Here the yield strength and related parameters from [START_REF] Schmidt | Some recent advances in the scaling of impact and explosion cratering[END_REF] The impact angle of the projectile, 𝛼, is assumed to be 45°, the most probable impact angle for projectiles hitting a surface [START_REF] Shoemaker | Interpretation of lunar craters[END_REF].

For large craters on big asteroids like Ceres, it is assumed in the Ivanov scaling law that the craters undergo collapse following the equation:

𝐷 M0:3<0 = 𝐷 3 hK7 𝐷 #S h (11)
Hereh = 0.15 and DSC is defined as the final rim diameter where simple craters transition into complex craters, which is assumed to be 10 km on Ceres [START_REF] Hiesinger | Cratering on Ceres: Implications for its crust and evolution[END_REF].

To verify that our coded versions of Eq. 10 and 11 function correctly for the results presented below, we reproduced the lunar and asteroid crater production functions shown in Fig.

3 [START_REF] Schmedemann | The cratering record, chronology and surface ages of (4) Vesta in comparison to smaller asteroids and the ages of HED meteorites[END_REF] (i.e., their normalized crater production curves for the Moon, Vesta (versions 3 and 4), Lutetia, Ida, and Gaspra). In this situation, input parameters were taken from their paper.

A general comment should also be made about this scaling law versus the Holsapple and Housen (2007) formulation. Both are based on the same general Pi-scaling theory and have a similar heritage (e.g., [START_REF] Schmidt | Some recent advances in the scaling of impact and explosion cratering[END_REF]. The difference is that the Holsapple and Housen (2007) scaling law as shown in Eq. 8 has free parameters for the various strength parameters that can be selected to match our best understanding of asteroid materials, whereas the Ivanov et al.

(2001) scaling law in Eq. 10 has those parameters built in. Presumably, one could select material parameters for the Holsapple and Housen (2007) scaling law to make it closer to the Ivanov et al.

(2001) scaling law, and one could reformulate the [START_REF] Ivanov | Size-frequency distributions of planetary impact craters and asteroids[END_REF] scaling law to have additional options as well. Therefore, the differences between the scaling laws are essentially choices in how asteroids are predicted to behave.

Empirical crater scaling law

A new element in this work is to use empirical methods to derive the appropriate asteroid scaling laws. Our method can be explained using two thought experiments.

For the first one, we consider a cumulative projectile and crater SFD defined as "broken" power laws; two power laws with different slopes that either meet at an inflection point or join each other over a slow bend (often called a "knee"). For this example, we assume the projectile and crater SFDs are not congruent. If one wanted to glean insights into the nature of the crater scaling law, the first thing to do would be to compare the inflection points or knees between the projectile and crater SFDs. Assuming that the crater scaling law is not pathological, these locations must correspond to one another. Their connection yields the relationship between the diameter of the projectile Dast and the diameter of the final crater Dcrater. We call this ratio

𝑓 = 𝐷 M0:3<0 𝐷 :/3 (12)
and use it throughout the paper. It is the simplest possible crater scaling law. In this example, there is only one value for f, but it can still be a powerful constraint if one desires to test crater scaling laws and impact models.

For the second thought experiment, we again assume that we have broken power laws for projectile and crater SFDs but that their shapes are congruent. By mapping the shape of the projectile SFD onto the crater SFD, one can empirically obtain the crater scaling law f for all sizes where data exist. In our idealized situation, no other information is needed; the myriad of crater scaling parameters for projectile and target properties are folded into the factor f.

When we started this project, we assumed that the first thought experiment was most likely to be applicable. As we show below, however, the projectile and crater SFDs used here are in fact excellent matches to the second thought experiment. This suggests that we can calculate empirical crater scaling laws for a wide range of crater sizes on different asteroids, provided their crater SFDs have a knee or that we have sufficient alternative constraints to rule out other possible scaling laws. As we will show, this method leads to powerful insights about the craters formed on different asteroids.

Using f values, one could presumably constrain more sophisticated crater scaling laws that describe how a given impact outcome is affected by different projectile and target quantities. The difficulty would be to overcome the degeneracy between the variables, such as the tradeoff between impact velocity, projectile size, etc. We do not perform such work here, but it would be an interesting follow-up project.

Collision Probabilities Between Target and Main Belt Asteroids (Component 3)

Two additional components are needed to model the collisional evolution of individual asteroids and interpret their crater histories: the intrinsic collision probabilities Pi and mean impact velocities Vimp of our target asteroids against the main belt population. There are many published formalisms to calculate these parameters that yield comparable results; a short list includes Ӧpik (1951), [START_REF] Wetherill | Collisions in the asteroid belt[END_REF], [START_REF] Kessler | Derivation of the collision probability between orbiting objects: the lifetimes of Jupiter's outer moons[END_REF], [START_REF] Farinella | Collision rates and impact velocities in the main asteroid belt[END_REF], [START_REF] Bottke | Velocity distribution among colliding asteroids[END_REF], Vedder

(1998), [START_REF] Manley | An algorithm for determining collision probabilities between small solar system bodies[END_REF]), Dell'Oro et al. (2001[START_REF] Vokrouhlický | Ӧpik-type collision probability for highinclination orbits[END_REF]. In this paper, we use the methodology of [START_REF] Bottke | Velocity distribution among colliding asteroids[END_REF].

For cratering events, the Pi parameter can be defined as the likelihood that a given projectile will hit a target with a given cross-sectional area over a unit of time. In most such cases, the size of the projectile is small enough to be ignored. For each pair of bodies, it can be considered to be the product of two combined probabilities:

• The probability that two orbits, with orbit angles that uniformly precess on short timescales, are close enough to one another that a collision can take place. It is the calculation of the volume of the intersection space of the pair of orbits.

• The probability that both bodies will be at their mutual orbital crossing location at the same time.

Our first task is to identify an appropriate projectile population that can hit our target asteroids. At that point, we compute individual Pi and Vimp values for all of the bodies on crossing orbits with the target. The (a, e, i) values of each pair are entered into the collision probability code, with the integral examining and weighting all possible orientations of the orbits, defined by their longitudes of apsides and nodes. This approximation is valid because secular perturbations randomize these values over ~10 4 -yr timescales.

The most difficult part of this task is finding the appropriate impactor population. Consider that most asteroid craters observed to date have been produced by projectiles smaller than the observational limit of the main belt (roughly Dast ~ 1-2 km). Moreover, the catalog of main belt objects suffers from observational selection effects, particularly as one approaches the observation limit. This makes it difficult to find a completely non-biased sample of main belt bodies for collision probability calculations.

Most asteroids discovered to date have been found by surveys limited in absolute magnitude (H). In general, for a given H value, it is easier to detect high-albedo S-type asteroids than low-albedo C-type asteroids. Any studies employing this catalog need to worry about overemphasizing S-types at the expense of C-types, particularly in the outer main belt where Ctypes dominate the population.

To mitigate against these problems, it is common to use a complete population of main belt asteroids as a statistical proxy for the population of smaller projectiles. For example, Bottke et al.

(1994) used 682 asteroids with Dast ≥ 50 km as defined by [START_REF] Farinella | Collision rates and impact velocities in the main asteroid belt[END_REF] for their collision probability calculations. The use of this sample is imperfect because family SFDs may be important at small asteroid sizes, but it may be the most reasonable approximation that we can make at this time, as we show below.

As a test, we also experimented with using the WISE diameter-limited catalog of main belt objects. The WISE catalog is incomplete, yet using it leads to results that are interesting in many ways. [START_REF] Masiero | Main belt asteroids with WISE/NEOWISE I: Preliminary albedos and diameters[END_REF] showed that the ratio of the number of outer to inner main belt asteroids becomes larger as one goes from Dast ≥ 50 km to Dast ≥ 10 km and then decreases again as one goes to Dast ≥ 5 km. The latter effect occurs because the power law slope of the inner main belt between 5 < Dast < 10 km is slightly steeper than that of the outer main belt over the same size range. The outer main belt appears to become observationally incomplete for Dast < 5 km, so we perform no calculations beyond this point.

This change in population has little effect on the collision probabilities of asteroids residing in the outer main belt, but it can be important for those in the inner main belt. As a demonstration of this effect, we selected (951) Gaspra for a series of Pi tests against the WISE catalog.

Using Gaspra's proper (a, e, i) values of (2.20974 au, 0.1462, 4.77253°) (Table 3), we calculated a mean Pi value for WISE asteroids on Gaspra-crossing orbits of 5.67, 5.15, and 5.39 × 10 -18 km -2 yr -1 for bodies of Dast ≥ 50, 10, and 5 km. Little change is seen between the values.

If we then fold in the population not on crossing orbits, which is needed to derive the approximate impact flux on Gaspra, the values change to 2.67, 1.74, and 2.11 × 10 -18 km -2 yr -1 , with 252 out of 535, 2289 out of 6754, and 15044 out of 38437 on crossing orbits. The low value for Dast > 10 km is notable, in that it is only 65% of the value for Dast ≥ 50 km. The mean Pi value then partially recovers for Dast ≥ 5 km because the inner main belt has a steeper SFD that the outer main belt. If we assume this trend holds to Dast ≥ 2 km, it seems likely that the mean Pi value for Gaspra will once again approach Dast ≥ 50 km. New work on debiasing the WISE asteroid catalog is needed to confirm this hypothesis.

PLACE TABLE 3 HERE

Given these trends and considerations, we argue a reasonable compromise is to continue to use the 682 asteroids with Dast ≥ 50 km discussed in [START_REF] Farinella | Collision rates and impact velocities in the main asteroid belt[END_REF] and Bottke et al. (1994) for our collision probability calculations of main belt bodies. Table 3 shows our results for all main belt asteroids observed by spacecraft. We obtained their proper (a, e, i) elements from the Asteroids Dynamic Site AstDyS, which is located at http://hamilton.dm.unipi.it/astdys/.

Results for Main Belt Asteroids and Near-Earth Asteroids Larger than 10 Kilometers

In this section, we examine the crater histories of the following asteroids observed by spacecraft: Vesta, Ceres, Lutetia, Mathilde, Ida, and Gaspra, which are all main belt asteroids, and Eros, a near-Earth asteroid. All have average diameters larger than 10 km. We start with main belt asteroids that have the largest size range of craters and work down to Gaspra, the smallest main belt asteroid in this list. Eros is actually larger than Gaspra, but we address it last to discuss the prospective relationship between Eros, Gaspra, and the Flora asteroid family.

While this list is long, it is not comprehensive. We avoid modeling certain asteroid terrains where crater saturation is prevalent, such as those on main belt asteroid Steins and the northern hemisphere of Vesta [START_REF] Marchi | Cratering on asteroids[END_REF]; see also Marchi et al., 2012a). To be cautious, we also decided to bypass sub-km craters in the crater SFD of Ceres, in part because they were potentially influenced by secondary cratering. Our analysis of those terrains is left for future work.

Finally, there are many proposed crater counts and crater retention ages for the asteroids or features discussed below. We focus here on published craters and ages that are most germane to testing our scaling laws and methods. For the interested reader who wants to know more, including a list of references about the craters found on these worlds, a good place to start would be to examine these papers: [START_REF] Chapman | Cratering on asteroids from Galileo and NEAR Shoemaker[END_REF], O'Brien et al. ( 2006), Schmedemann et al.

(2014), and [START_REF] Marchi | Cratering on asteroids[END_REF].

Rheasilvia Basin on Vesta

(4) Vesta is the second largest main belt asteroid. It is located in the middle of the inner main belt with proper orbital elements of (a, e, i) = (2.36 au, 0.099, 6.4°). NASA's Dawn spacecraft imaged its surface at varying spatial resolutions and verified that Vesta had differentiated into a metallic core, silicate mantle, and basaltic crust. Some key physical parameters for Vesta include dimensions of 572.6 km × 557.2 km × 446.4 km, a bulk density of 3.456 ± 0.035 g cm -3 , and a surface gravity of 0.25 m s -2 [START_REF] Russell | Dawn at Vesta: Testing the protoplanetary paradigm[END_REF].

Here we re-examine the superposed crater SFD on or near Vesta's Rheasilvia basin, a 500km-diameter impact structure that defines the shape of Vesta's southern hemisphere (e.g., Schenk et al., 2012) (Fig 3). We choose this region for our modeling work for two reasons: Rheasilvia is young enough that crater saturation is not an issue, and it is broad enough that it is covered by a large range of crater diameters (0.15 < Dcrater < 35 km; [START_REF] Marchi | Cratering on asteroids[END_REF].

PLACE FIGURE 3 HERE

The craters identified and used here are located on Rheasilvia's floor and ejecta blanket;

their properties are reported in [START_REF] Marchi | Cratering on asteroids[END_REF] (see also Marchi et al., 2012a for earlier counts). Their work indicated that a plausible age for Rheasilvia was ~1 Ga (Marchi et al., 2012a).

Model components that went into this age include (i) the main belt SFD described by Bottke et al. (2005b) (SFD #1 in Fig. 1), (ii) an intrinsic collision probability between main belt asteroids and Vesta of Pi = 2.8 × 10 -18 km -2 yr -1 , and (iii) the Holsapple and Housen (2007) scaling law for cohesive soils (Y = 2 × 10 7 dynes cm -2 ).

Their estimated crater retention age for Rheasilvia is comparable to the 40 Ar/ 39 Ar ages of feldspar grains in the brecciated howardite Kapoeta, which were reset by a thermal event between 0.6 and 1.7 Ga ago (Lindsey et al., 2015). Lindsey et al. suggested that the source of the heating event was the formation of Rheasilvia basin 1.4 ± 0.3 Ga ago. They also pointed out that this age is similar to other 40 Ar/ 39 Ar ages found among the HED (howardite-eucrite-diogenite) meteorites.

A range of ages between 0.6 and 1.7 Ga seem plausible given these data.

Note that 40 Ar/ 39 Ar ages between 3.5-4.1 Ga have also been identified in eucrites. These ages are older than the crater retention ages found for Rheasilvia by [START_REF] Schenk | The geologically recent giant impact basins at Vesta's south pole[END_REF] and Marchi et al. (2013b). Using their own crater counts, and comparing their model to craters with a more limited dynamic range than those works, [START_REF] Schmedemann | The cratering record, chronology and surface ages of (4) Vesta in comparison to smaller asteroids and the ages of HED meteorites[END_REF] argued that Rheasilvia had a crater retention age that matched those ancient values. We will address this issue below.

The Rheasilvia basin-forming event also ejected numerous fragments onto escape trajectories, and these bodies likely comprise Vesta's color-, spectral-and albedo-distinctive asteroid family (e.g., [START_REF] Parker | The size distributions of asteroid families in the SDSS Moving Object Catalog 4[END_REF][START_REF] Nesvorný | Identification and dynamical properties of asteroid families[END_REF][START_REF] Masiero | Asteroid family physical properties[END_REF]. Using a collisional evolution model, [START_REF] Bottke | On the Origin and Evolution of Vesta and the V-Type Asteroids. Vesta in the Light of Dawn: First Exploration of a Protoplanet in the Asteroid Belt[END_REF] found that the Vesta family's steep SFD, composed of bodies of Dast < 10 km, showed no indication of a change produced by collisional grinding. On this basis, they estimated that the Vesta family has an 80% probability of being < 1 Ga old. The orbital distribution of the family members, and how they have likely been influenced by the Yarkovsky thermal forces, also suggests an age of ~1 Ga [START_REF] Spoto | Asteroid family ages[END_REF], though we caution that the high ejection velocity of the family members makes it difficult to precisely determine the family's dynamical age (e.g., [START_REF] Nesvorný | Identification and dynamical properties of asteroid families[END_REF].

The combined crater sets of Rheasilvia presented in [START_REF] Marchi | Cratering on asteroids[END_REF] yield 48 craters between 0.15 < Dcrater < 35 km (Fig. 4). Two knees are seen in the crater SFD: one near Dcrater ~ 2 km and a second near Dcrater ~ 20 km. The smaller of the two knees is likely related to the inflection points seen between 0.2 < Dast < 0.6 km in the main belt SFDs shown in Fig. 1. The origin of the larger knee will be discussed below.

PLACE FIGURE 4 HERE

Empirical scaling law derived by fitting model and observed crater SFDs (Rheasilvia)

To compare the shape of Vesta's crater SFD to the impactor SFDs shown in Fig. 4, we defined two parameters: (i) the crater scaling relationship factor f = Dcrater / Dast and (ii) the age of the Rheasilvia surface Tast. The number of model craters forming per square kilometer on the surface of the asteroid, Nmodel-crater (> Dcrater), as a function of time 𝑇 :/3 is given by the equation:

𝑁 5R-<T=M0:3<0 (> 𝐷 M0:3<0 ) = 𝑃 . 𝑇 :/3 𝑁 5R-<T=:/3 (> 𝐷 :/3 ) 4𝜋 (13)
The number of model asteroids larger than a given size Dast is given by Nmodel-ast (> Dast), which can be found in Fig. 1.

The quality of the fit between the observed crater SFD on Rheasilvia (Fig. 4) and those modeled is defined using chi-squared methods:

𝜒 6 = ] ^𝑁5R-<T=M0:3<0 (> 𝐷 . ) -𝑁 RU/=M0:3<0 (> 𝐷 . )_ 6 𝑁 RU/=M0:3<0 (> 𝐷 . ) B .V7 (14)
Here 𝐷 . = 1,..., M, stands for the diameters of observed and model craters on a given asteroid surface. To obtain normalized 𝜒 6 values, one should divide them by the value M, yielding the value we define here as 𝜒 6 ER05 . In this case, there were 48 Rheasilvia craters, so M = 48 (Table 4).

PLACE TABLE 4 HERE

By creating an array of (f, Tast) values 5 < f < 25, incremented by 0.1, and 0.01 < Tast < 5

Ga, incremented by 0.01 Ga, we were able to test all plausible fits between model and observed crater data. These values also allow us to calculate confidence limits of 68% (1s) and 95% (2s) relative to our best fit case that can be used to estimate error bars.

An additional issue with fitting a model SFD to a crater SFD is that the smallest craters in N (> Dcrater), which have the most data and the smallest error bars but also are closest to the observation limit, tend to dominate the 𝜒 6 values. To mitigate against this effect, we multiplied the error bars of N (> Dcrater) by a function 𝛾 that increases the error bars of the smaller craters according to:

𝛾 = (𝑤 -1) log 7N e 𝐷 M0:3<0 𝐷 M0:3<0 5:W f log 7N / 𝐷 M0:3<0 5.E 𝐷 M0:3<0 5:W 3 + 1 (15)
Here we set w to 3-5 for the asteroid craters in this paper, with 𝐷 M0:3<0

5.E

and 𝐷 M0:3<0 5:W defined by the minimum and maximum crater sizes in a given set, respectively.

Using the Pi value in We also examined how our main belt SFDs compared to Rheasilvia's superposed craters using the Holsapple and Housen (2007) formulation of the Pi-group crater scaling law (Eqs. 8 and 9). Hereafter we call this the HH crater scaling law.

Following the procedure used by [START_REF] Marchi | Cratering on asteroids[END_REF], and applying his chosen parameters for Eqs. 8 and 9, we assumed that Vesta's surface could be treated like it had the same material properties as cohesive soils (k = 1.03, ν = 0.4, µ = 0.41). We assumed that the projectile density was 2.5 g cm -3 and Vesta's surface density was 3.0 g cm -3 . After some trial and error, we found that the lowest 𝜒 6 BC values were generated from Y = 2 × 10 7 dynes cm -2 . We use this value for all of the asteroids discussed below. The values of Pi and Vimp are found in Table 3.

Our best fit came from SFD #7, which yielded 𝜒 6 B) = 19.34 (Table 4). This value indicates that our fit here is modestly superior to the empirical scaling law results in Sec. 3.1.1. The reason is that their f value decreases for larger projectiles, allowing SFD #7 to match Rheasilvia's craters with Dcrater > 20 km.

Our apparent success for large crater sizes, however, may be an issue. Existing numerical hydrocode simulations indicate that the 500-km Rheasilvia basin formed from the impact of a 37to 60-km-diameter projectile [START_REF] Jutzi | The structure of the asteroid 4 Vesta as revealed by models of planet-scale collisions[END_REF][START_REF] Ivanov | Two-dimensional numerical modeling of the Rheasilvia impact formation[END_REF], which corresponds to f = 8.3-13.5. These latter values are a good match to the f values predicted by our empirical scaling law results (Fig. 4, Table 4).

If f indeed decreases substantially for large craters, as shown by the red curve in the inset figure within Fig. 4, it would imply that much larger impactors-perhaps Dast > 100 km-would be needed to make Rheasilvia. We find this to be an unlikely scenario. We discuss this issue further in Sec. 3.1.4.

Our best fit crater retention age for this set of parameters is Tast = 1.24 [-0.06, + 0.06] Ga (Table 5). This value matches the ages of Marchi et al. (2012a) and 40 Ar/ 39 Ar ages constraints from Lindsey et al. ( 2015), but it is modestly older than the empirical fit results in Sec. 3.1.1. The reasons are that (i) this scaling law yields f values that are consistently lower than the empirical main belt best fit results of f = 9.9, which increases the surface age, and (ii) SFD #7 is shallower at small asteroid sizes and therefore has fewer small projectiles; few projectiles mean older ages.

Ivanov crater scaling fit (Rheasilvia)

The last scaling law investigated was that from [START_REF] Ivanov | Size-frequency distributions of planetary impact craters and asteroids[END_REF] [START_REF] Spoto | Asteroid family ages[END_REF] or collisional evolution [START_REF] Bottke | On the Origin and Evolution of Vesta and the V-Type Asteroids. Vesta in the Light of Dawn: First Exploration of a Protoplanet in the Asteroid Belt[END_REF]. As before, the reason has to do with the large f values applied here; if smaller projectiles make larger craters, the surface has to be younger.

Our crater retention ages are different than those of [START_REF] Schmedemann | The cratering record, chronology and surface ages of (4) Vesta in comparison to smaller asteroids and the ages of HED meteorites[END_REF], who use the same scaling law to get 3.5 ± 0.1 Ga (though some surfaces have reported ages of 1.7-1.8 Ga;

see their Table 6). Only a minor portion of this difference can be attributed to their use of different collision probabilities or impact velocities; their values are nearly the same as the ones we show in Table 3. Similarly, in our tests of their work, we find that their derived main belt SFD is similar to our SFD #8 in Fig. 1.

The main reason that Schmedemann et al. (2014) report a different crater retention age for

Rheasilvia than our work is that they focus on comparing their crater production function to craters sizes between several kilometers and several tens of kilometers in diameter. As our Fig. 4 shows, if one ignored all craters smaller than a few km, the best fit curve for the Ivanov scaling law would shift to substantially higher values, with f values that are closer to those of the HH scaling law and the empirical main belt fit scaling law. These effects would in turn yield substantially older crater retention ages.

There may be valid reasons why one should ignore fitting a crater production function model to small craters on a given surface, and it is possible that one can obtain reasonable results

by only looking at the largest craters on a surface. Nevertheless, where practicable, it is better to compare a crater model across an entire size range of craters rather than a subset. For this reason, we argue that this crater scaling law does not perform as well at modeling Vesta's crater SFD in Fig. 4 as the other choices.

What projectile sizes make the largest basins on Vesta?

In section 3.1.2, we asserted that we favor f ~ 10 to make the largest basins on Vesta. One reason is that these values are consistent with hydrocode simulations, where the 500-km Rheasilvia basin formed from the impact of a 37-to 60-km-diameter projectile [START_REF] Jutzi | The structure of the asteroid 4 Vesta as revealed by models of planet-scale collisions[END_REF][START_REF] Ivanov | Two-dimensional numerical modeling of the Rheasilvia impact formation[END_REF]. A second reason is that it matches f values predicted by our empirical scaling law results (Fig. 4, Table 4). A third reason comes from the following calculation (see also Bottke et al., 2015a).

The two largest basins on Vesta are Rheasilvia and Veneneia, with diameters of ~500 and ~400 km, respectively. Veneneia is partially buried by Rheasilvia, so its estimated crater retention age is > 2 Ga [START_REF] Schenk | The geologically recent giant impact basins at Vesta's south pole[END_REF]. Both formed after the emplacement of Vesta's basaltic crust, which solidified within a few million years of Solar System formation (e.g., [START_REF] Hopkins | Micrometer-scale U-Pb age domains in eucrite zircons, impact re-setting, and the thermal history of the HED parent body[END_REF].

In the following calculation, we will assume these basins were produced by f = 5, which would require impacts from Dast ≥ 90 km and 100 km bodies, or f = 10, which would require impacts from Dast ≥ 40 km and 50 km bodies. If we assume the main belt had approximately the same population of large asteroids over the last 4.5 Gyr as it has today, such that we can use the SFD in Fig. 1 over this interval, the population of Dast ≥ 40, 50, 90, and 100 km asteroids in the main belt is 860, 680, 270, 220, respectively. Using Pi = 2.8 × 10 -18 km -2 yr -1 , the probability that these projectiles will collide with Vesta over 4.5 Gyr is 0.79, 0.65, 0.30, 0.25 for Dast ≥ 40, 50, 90, and 100 km, respectively.

Accordingly, for f = 5, the probability that both Rheasilvia and Veneneia will form on Vesta is 0.30 × 0.25 = 7%, while for f = 10, the probability is 0.79 × 0.65 = 51%. The latter value is 7 times higher than the former. This does not mean the real scaling law must be f ~ 10, but it is fair to say that f ~ 10 makes it easier to match constraints.

Summary (Rheasilvia)

Our results for Vesta allow us to make some initial observations about how to interpret crater SFDs on asteroids.

• Results for the empirical fit and HH crater scaling laws indicate that the Rheasilvia formation event probably took place between ~0.6 Ga and ~1.3 Ga ago. This result is consistent with dynamical, collisional evolution, and meteorite constraints.

• If an asteroid's crater SFD has a knee near Dcrater ~ 2 km, it indicates that an f value near ~10

will allow it to match the main belt SFD. If f values are substantially smaller or larger than 10, one can only fit the small craters in the SFD at the expense of missing the larger ones, or vice versa.

• Main belt asteroid SFDs that remain shallow between ~0.2-0.3 < Dast < 2-3 km, results that are represented by SFDs #4 to #8 in Fig. 1, appear to be the most successful at matching constraints on Vesta.

• Using the Holsapple and Housen (2007) formulation of the Pi-group scaling law, we find that input parameters for cohesive soils appear to allow us to best match observations.

• Scaling laws that have f values substantially smaller than 10 for craters of Dcrater ~ 20 km on Rheasilvia are needed to match data, but an extrapolation of this trend would make it difficult to produce Rheasilvia and Veneneia basins. Our interpretation is that this may make f ~ 10 a reasonable choice for all of the observed crater data on or near Rheasilvia basin.

Ceres's Kerwan Basin

(1) Ceres, with dimensions of 965.2 ± 2.0 km × 961.2 ± 2.0 km × 891.2 ± 2.0 km, is the largest asteroid in the main belt (Park et al. 2016) The nature of the craters on Ceres suggests that its surface may be intermediate in strength between that of Vesta and Rhea, the icy satellite of Saturn [START_REF] Russell | Dawn arrives at Ceres: Exploration of a small, volatile-rich world[END_REF]. The lack of crater relaxation observed for smaller craters, however, indicates that the crust may be deficient in ice, and could be a mechanically strong mixture of rock, carbonates or phyllosilicates, ice, and salt and/or clathrate hydrates [START_REF] Fu | The interior structure of Ceres as revealed by surface topography[END_REF]. Curiously, Ceres is missing very large craters (Dcrater > 280 km) and is highly depleted in craters of diameter 100-150 km compared to expectations from the shape of the impacting main belt SFD [START_REF] Marchi | The missing large impact craters on Ceres[END_REF]. Their absence could suggest the viscous relaxation of long-wavelength topography, perhaps via a subsurface zone of low-viscosity weakness [START_REF] Fu | The interior structure of Ceres as revealed by surface topography[END_REF].

To glean insights into the nature of Ceres's crust, we examine the superposed crater SFD associated with Ceres's Kerwan basin (Fig. 5). Geologic mapping work indicates that Kerwan is the oldest, largest (undisputed) impact crater on Ceres (Dcrater ∼284 km) (Williams et al., 2017).

The derived age of the basin depends on the superposed crater counts and the crater age model used (see the crater SFDs from Williams et al., 2017), but craters counted in the smooth unit of Kerwan, which range from approximately 5 < Dcrater < 100 km, yield ages of 550 ± 90 Ma and 720 ± 100 Ma [START_REF] Hiesinger | Cratering on Ceres: Implications for its crust and evolution[END_REF]. With that said, none of the model crater SFDs shown in Fig. 8 of [START_REF] Hiesinger | Cratering on Ceres: Implications for its crust and evolution[END_REF] appear to reproduce the shape of the crater SFD, and the above ages seem to be determined by best fits to the largest craters.

PLACE FIGURE 5 HERE

Here we compare our model crater SFDs to craters counted by co-author S. Marchi. They have approximately the same crater SFD as [START_REF] Hiesinger | Cratering on Ceres: Implications for its crust and evolution[END_REF]. The crater counts have been slightly updated and are shown in Fig. 6. The observed inflection point in these crater SFD occurs near Dcrater ~ 20 km, approximately the same size as seen for Vesta's Rheasilvia basin (Fig. 4).

Accordingly, our prediction is that the Kerwan and Rheasilvia basins should have similar crater scaling laws.

PLACE FIGURE 6 HERE

Empirical scaling law derived by fitting model and observed crater SFDs (Kerwan)

Using our empirical main belt fit method, and a Pi value of 3.455 × 10 -18 km -2 yr -1 (Table 3), we found that our best fit comes from SFD #8, which yielded 𝜒 6 BC = 2.84 (

Housen and Holsapple crater scaling fit (Kerwan)

Our comparison between model and data using the Holsapple and Housen (2007) formulation of the Pi-group scaling law also produced a reasonable match with Rheasilvia craters (Fig. 4). As with our Vesta runs, we assumed that Ceres's surface had the strength of cohesive soils (k = 1.03, ν = 0.4, µ = 0.41), with projectile and target surface density set to 2.5 and 1.5 g cm - 3 , respectively. The yield strength was Y = 2 × 10 7 dynes cm -2 , and the values of Pi and Vimp are found in Table 3. This scaling law yielded a best fit using SFD #8, with 𝜒 6 B) = 4.61, modestly higher than the empirical main belt fit of 𝜒 6 BC = 2.84. Here the slight mismatch stems from the model crater SFD missing the craters with Dcrater > 40 km. The best fit crater retention age in this circumstance is Tast = 0.70 [-0.02, + 0.03] Ga, fairly close to the main belt fit result (Table 5).

Ivanov scaling law fit (Kerwan)

The fit using the Ivanov crater scaling law is only modestly better than the HH scaling law case, with 𝜒 6 "X = 9.44. The f values are larger here across the board than the other two scaling law cases, with crater sizes increasing substantially for Dast < 1 km. The best fit crater retention age here is 𝑇 :/3 = 0.47 [-0.01, + 0.02] Ga, younger than the previous two test cases (Table 5).

With this said, there are indications in other Kerwan crater databases not investigated in this paper that the slope of the observed crater SFD indeed becomes substantially steeper for craters with Dast < 1 km, as predicted by the Ivanov scaling law (e.g., Williams et al., 2017). If so, the crater SFD on Kerwan is radically different than the one observed in Vesta's Rheasilvia basin. It seems unlikely that the main belt SFD changed over the timescales in question, so this mismatch between Vesta and Ceres implies that the f value function is different for sub-kilometer craters on the two worlds.

Finally, we note that it is plausible that the observed difference between sub-kilometer craters on Kerwan and Rheasilvia is because secondary craters are pervasive across Ceres for Dcrater < 2-3 km.

Summary (Kerwan)

As with Rheasilvia basin on Vesta, the best fit matches in Table 4 come from the empirical main belt fit and the HH scaling law fit. As before, our results favor the high-number SFDs, with the best fit coming from #8, a value that is modestly disfavored from probability studies (Table 1).

If we use a more probable SFD, such as SFD #6 (Table 1), the empirical scaling law f ~ 9, with errors that overlap with our f ~ 10 solution for Rheasilvia basin (Table 5). The age of Kerwan basin favored by our results is ~0.8-0.9 Ga.

Lutetia's Achaia Region

(21) Lutetia is an M-type asteroid with dimensions of 121 ± 1 km × 101 ± 1 km × 75 ± 13 km [START_REF] Sierks | Images of Asteroid 21 Lutetia: A Remnant Planetesimal from the Early Solar System[END_REF]) (Fig. 7). It is located in the inner main belt and has proper orbital elements of (a, e, i) = (2.43 au, 0.13, 2.1°). The flyby of Lutetia by ESA's Rosetta mission yielded a bulk density of 3.4 ± 0.3 g cm -3 [START_REF] Sierks | Images of Asteroid 21 Lutetia: A Remnant Planetesimal from the Early Solar System[END_REF]). The composition of Lutetia is unknown, though it is thought to be related to enstatite chondrites or possibly the metal-rich CH carbonaceous chondrites [START_REF] Coradini | The surface composition and temperature of asteroid 21 Lutetia as observed by Rosetta/VIRTIS[END_REF][START_REF] Moyano-Cambero | A plausible link between the asteroid 21 Lutetia and CH carbonaceous chondrites[END_REF].

PLACE FIGURE 7 HERE Marchi et al. (2012b) examined craters on the oldest observed surface imaged by Rosetta, a flat and uniform region called Achaia (Fig. 7). The craters on Achaia range in size from ~1 < Dcrater < 50 km (Fig. 8), a large enough dynamic range that they potentially sample both inflection points in our main belt SFD (Fig. 1).

PLACE FIGURE 8 HERE

Empirical scaling law derived by fitting model and observed crater SFDs (Achaia)

Our main belt fit method, combined with a Pi value of 3.763 × 10 -18 km -2 yr -1 (Table 3), yields a best fit for SFD #7 with 𝜒 6 BC = 1.97 (Table 4). This value is the best of the three scaling laws tested for this case. The SFDs #3-#8 are also yield results within 1s of this best fit case. Our best fit f value was 10.30 [-4.90, +3.60], similar to results from Vesta and Ceres. Our best fit crater retention age is Tast = 2.57 [-2.02, +1.64] Ga. This mean value is lower than the ~3.6-Ga age estimate from Marchi et al. (2012b) but overlaps within errors. Visually, our main belt fit curve hits both inflection points and is a good visual fit to the data.

Housen and Holsapple crater scaling fit (Achaia)

Using the HH scaling law, the same yield strength as above, Pi and Vimp in Table 3, and projectile and target surface densities of 2.5 and 3.0 g cm -3 , respectively, we obtain a comparable but slightly worse fit, with 𝜒 6 B) = 2.24 for SFD #7. Other SFDs within 1s of the best fit case are #5-8. The crater retention age from this fit is Tast = 3.07 [-0.41, + 0.41] Ga, similar to the age found in Marchi et al. (2012b). The empirical main belt and HH crater scaling curves in Fig. 8 are fairly similar to one another. The difference in age is produced by modest differences in the f function, with the HH scaling curve having f values between 7 and 9, somewhat lower than the main belt fit of f ~ 10.

Ivanov scaling law fit (Achaia)

The Ivanov scaling law produces the poorest fit, with 𝜒 6 "X = 4.00 for SFD #8 (Table 4).

As with the Rheasilvia region on Vesta, the typical f value near 15 for smaller projectiles is higher than those found for either of the other two scaling laws (Fig. 8). This large value causes the best fit case to undershoot crater data between 3 < Dcrater < 15 km. It also leads to a younger age for Achaia, with Tast = 1.24 [-0.16, +0.16] Ga. [START_REF] Schmedemann | The cratering record, chronology and surface ages of (4) Vesta in comparison to smaller asteroids and the ages of HED meteorites[END_REF] report that their fit to Lutetia's large and degraded craters, which have diameters between 2 and 25 km, yields a crater retention age of 3.5 ± 0.1 Ga. We suspect the difference in age between the two calculations is caused by our use of craters between 1 and 3 km in Fig. 8. These craters drive our fit. If small craters were ignored, the best fit for the Ivanov scaling law curve would slide upward, which in turn would correspond to an older age.

Summary (Achaia)

In all three cases, the best fit matches seem to come from high-number SFDs, a common theme for all of the crater SFDs discussed up to this point. The main belt fit with f ~ 10 and the HH scaling law fit are preferred from the chi-squared metric over the Ivanov scaling law fit, probably because the latter's model crater SFD is lower than the observed crater data for middlesized craters in Fig. 8. The crater retention age of the Achaia region has a wide range of possible ages, but a reasonable value for this terrain is ~3-4 Ga.

Mathilde

(253) Mathilde was the first C-type asteroid imaged by spacecraft (Fig. 9). NASA's Near

Earth Asteroid Rendezvous -Shoemaker mission (NEAR) flew by it in 1997 en route to Eros. It is located in the middle of the central main belt, with proper orbital elements of (a, e, i) = (2.65 au, 0.22, 6.5°). Its physical dimensions are 66 km × 48 km × 46 km, and it has an estimated bulk density of 1.3 ± 0.2 g cm -3 [START_REF] Veverka | NEAR's flyby of 253 Mathilde: Images of a C asteroid[END_REF]1999). Mathilde is not in an asteroid family [START_REF] Nesvorný | Identification and dynamical properties of asteroid families[END_REF] and thus its surface may have borne witness to the early days of main belt history.

PLACE FIGURE 9 HERE

Mathilde, Bennu, and Ryugu are all C-complex bodies and thus are thought to have physical properties similar to carbonaceous chondrite meteorites. Accordingly, the crater history of Mathilde may provide us with the most direct insights into the crater scaling laws that govern Bennu and Ryugu. There are different schools of thought about how cratering should work on Mathilde. The flyby images that we have of Mathilde are dominated by its two largest craters, Ishikari (29.3 km) and Karoo (33.4 km) (Fig. 9). There appears to be little ejecta surrounding these craters [START_REF] Veverka | NEAR encounter with asteroid 253 Mathilde: Overview[END_REF][START_REF] Chapman | Cratering on Mathilde[END_REF]. [START_REF] Housen | Compaction as the origin of the unusual craters on the asteroid Mathilde[END_REF] suggested that the apparent absence of deep ejecta blankets indicates the cratering process is dominated by the effects of porous materials. In such media, craters form more by compaction than excavation, and what little ejecta is produced goes back into the central cavity [START_REF] Housen | Impact cratering on porous asteroids[END_REF][START_REF] Housen | Impacts into porous asteroids[END_REF].

The final size of the crater made into such a target, however, is unclear. Kinetic impact energy transfer can be inefficient in a porous target, and this may result in craters that are not much larger than those formed in targets with large yield stresses. Alternatively, the low strength nature of the target may make it easy for an impact to push material out of the way (e.g., comparable to impacts into sand targets; O' Brien et al., 2006), and this could result in a larger craters per projectile diameter than those found on S-type asteroids.

To glean insights into this issue, it is useful to examine Mathilde's crater SFD. Here we adopt the cumulative crater counts and errors of O'Brien et al. ( 2006), who converted their data from Fig. 3 of [START_REF] Chapman | Cratering on Mathilde[END_REF] (Fig. 10). Given how far crater counting tools have advanced in two decades, Mathilde would seem to be a ripe target for a re-examination.

It is plausible that the Mathilde data are in saturation at smaller crater sizes, as suggested by [START_REF] Chapman | Cratering on Mathilde[END_REF], but the cumulative power law slope of the data for Dcrater < 2 km is close to q = -2.6, like those of Rheasilvia basin (Vesta) and the Achaia region (Lutetia) (Figs. 4,8). Our expectation is that a crater SFD in saturation would instead have a cumulative slope of q = -2 (e.g., Melosh 1989), a value that is possible within error bars but does not appear to be the true solution.

Here we will assume that Mathilde's observed crater SFD is not in saturation.

PLACE FIGURE 10 HERE

A relatively recent calculation of the crater retention age of Mathilde was made by [START_REF] O'brien | Craters on asteroids: Reconciling diverse impact records with a common impacting population[END_REF]. They made the following assumptions: (i) the intrinsic collision probability Pi for Mathilde was 2.86 × 10 -18 km -2 yr -1 (a factor of 0.77 lower than our value from Table 3), (ii) the main belt SFD followed the estimate made by O'Brien and Greenberg (2005), (iii) Mathilde's craters could be produced by the Pi-group crater scaling relationship, provided that the target acted like loose sand, and (iv) sandblasting by small impactors, a potential crater erasure mechanism, was active on Mathilde (e.g., [START_REF] Greenberg | Collisional history of Gaspra[END_REF]1996). Their crater scaling law for Mathilde is shown in their Fig. 2. For their impactors of 0.1 to 1 km diameter, their value for f was approximately 20-40.

O 'Brien et al. (2006) report that the population of km-scale and smaller craters becomes saturated at a constant level near 1 Ga ago, but that the best crater retention age for the large craters is ∼4 Ga. We will instead assume below that Mathilde's crater SFD follows a production population and that crater erasure mechanisms are not needed to explain observations.

Empirical scaling law derived by fitting model and observed crater SFDs (Mathilde)

Assuming Pi = 3.723 × 10 -18 km -2 yr -1 (Table 3), we obtain a best fit for SFD #2 with 𝜒 6 BC = 2.67 (Table 4). We find that SFDs #1-6 all have 𝜒 6 BC within 1s of this best fit case. Our best fit value for f is 10.00 [-3.40, +1.90], and best fit crater retention age was Tast = 3.70 [-1.30, +1.52].

If we were to adopt SFD #6, these values change to f = 12.2 [-2.0, +2.3], with Tast becoming 3.61

[-1.31, +0.94] Ga. Here the smaller number of projectiles in SFD #6 is compensated by having each projectile make a modestly larger crater.

Housen and Holsapple crater scaling fit (Mathilde)

The parameters listed in Holsapple and Housen (2007) for the highly porous case are essentially the same as those used for cohesive soils. It is argued in their paper that their experimental craters were governed predominantly by some compressive strength, so strengthscaled laws would apply. Accordingly, their scaling is almost identical to those used for other asteroids and assumes that compressive strength is the relevant measure.

Assuming a target and projectile density of 1.3 g cm -3 , we obtain a best fit of 𝜒 6 B) = 2.60

for SFD #4. The f value of the scaling law is close to f = 12 for the projectile sizes used here (Fig. 10). The SFDs #1-#5 are 1s of the best fit case, and they yield ages of 2.2 < Tast < 3.4 Ga. If we include error bars, the age range expands to 1.8 < Tast < 4 Ga.

Ivanov scaling law fit (Mathilde)

The Ivanov scaling law leads to a best fit of 𝜒 6 "X = 4.98 for SFD #8 (Table 4). Its f values are substantially higher than the other two cases, and they yield a crater retention age of Tast = 0.85 [-0.15, +0.15] Ga. This outcome shows the implications of having projectiles make much larger craters on Mathilde relative to other scaling laws; the surface age of the body becomes young enough that one needs to invoke a special event to explain why the surface is young. As stated above, Mathilde has no asteroid family, and its big craters formed apparently without damaging one another. This makes it more difficult to argue that impacts have reset the surface relatively recently (i.e., over the past billion years).

Summary (Mathidle)

The crater SFD for Mathilde is murky enough that it is difficult to identify which fits are best. Until stronger evidence becomes available, it seems reasonable to adopt the simplest solutions.

Accordingly, if we apply a high-number SFD to our crater model, one that also provides a good fit to the other asteroid terrains investigated so far (i.e., SFD #6, which for Mathilde yields results that are within 1s of the best case fit for both the empirical main belt fit case and the HH scaling law case), our results suggest that f ~ 10-12 provides a good solution to the entire crater we argue that a young crater retention age for Mathilde seems unlikely.

Ida

The second asteroid observed by the Galileo spacecraft was (243) Ida, an S-type asteroid that is also a member of the Koronis asteroid family (Fig. 11). It is located in the outer main belt, with proper orbital elements of (a, e, i) = (2.86 au, 0.05, 2.1°). Its physical dimensions are 59.8 km × 25.4 km × 18.6 km [START_REF] Belton | Galileo's encounter with 243 Ida: An overview of the imaging experiment[END_REF]. The best estimate of Ida's bulk density comes from

PLACE FIGURE 11 HERE

The age of the Koronis family, and likely that of Ida, can be computed using dynamical models. By tracking how Koronis family members have drifted in semimajor axis by the combined Yarkovsky and YORP (Yarkovsky-O'Keefe-Radzievskii-Paddack) effects, it has been estimated

that the family's age is ~2-3 Ga [START_REF] Bottke | Dynamical spreading of asteroid families via the Yarkovsky effect: The Koronis family and beyond[END_REF][START_REF] Brož | Constraining the cometary flux through the asteroid belt during the Late Heavy Bombardment[END_REF][START_REF] Spoto | Asteroid family ages[END_REF]; see also [START_REF] Nesvorný | Identification and dynamical properties of asteroid families[END_REF]. Error bars on this value could take it to 4 Gyr as well.

There are also multiple Koronis family members (including Ida) whose spin vectors have been modified enough by the YORP effect, some to be caught in spin-orbit resonances [START_REF] Vokrouhlický | The vector alignments of asteroid spins by thermal torques[END_REF]. Objects with prograde spins in so-called "Slivan states" have nearly identical periods (7.5-9.5 h), obliquities between 42° and 50°, and pole longitudes confined in a tight interval between 25° and 75°. [START_REF] Vokrouhlický | The vector alignments of asteroid spins by thermal torques[END_REF] estimated that the time needed for the observed bodies with 20 < Dast < 40 km, starting with spin periods P = 5 hours, to reach their Slivan state status was ~2-3 Ga. Accordingly, our expectation is that Ida's crater retention age should be close to this value.

Modeling Ida's crater history is challenging because many smaller crater sizes appear to be near or in saturation equilibrium (Chapman et al., 1996b;reviewed in Chapman et al., 2002 and[START_REF] Marchi | Cratering on asteroids[END_REF]. There are ways to deal with crater saturation using specialized codes (e.g., Marchi et al. 2012a), but they also involve making assumptions about the nature of the saturation process; this issue will be discussed in a follow-up paper.

As a work-around, we examine the crater counts provided by Fig. 5 of Chapman et al.

(1996b) (Fig. 12). Their crater SFD between 0.6 < Dcrater < 10 km has a similar shape to those seen on Vesta and Lutetia for the same size range (Figs. 4 and8). Only the craters with Dcrater < 0.6 km seem to have the -2 cumulative power law diagnostic of saturation (e.g., Melosh 1996; see also Marchi et al. 2012a). Our main belt model will be applied to Ida's craters with Dcrater > 0.6 km.

PLACE FIGURE 12 HERE

Empirical scaling law derived by fitting model and observed crater SFDs (Ida)

Using a Pi value of 4.037 × 10 -18 km -2 yr -1 (Table 3), we obtain a best fit for SFD #8 with 𝜒 6 BC = 1.17 (Table 4). All of the other SFDs also fit within 1s of this best fit case, a byproduct of the limited number of large craters in our crater SFD. As in the case for Lutetia's Achaia region, the empirical main belt fit case is the best of the three scaling laws tested.

Our should not be taken literally. Both crater retention ages for Ida are consistent with the estimated dynamical ages of the Koronis family. Visually, our main belt fit curve hits both inflection points and provide a good visual match to the data.

Housen and Holsapple crater scaling fit (Ida)

For the HH scaling law case, we applied the same yield strength as before, Pi and Vimp values from Table 3, and projectile and target surface densities of 2.5 and 2.6 g cm -3 , respectively.

Our best fit case was for SFD #6, which yielded 𝜒 6 B) = 1.39 (Table 4). As in the empirical main belt fit case, all SFDs are within 1s of the best fit case. The crater retention age from our best fit is Tast = 2.91 [-0.43, +0.43] Ga, which again is a reasonable match with the dynamical age of the Koronis family. The f functions for this scaling law are slightly lower than the main belt fit, but this is balanced against a different best fit choice for the main belt SFD.

Ivanov scaling law fit (Ida)

The Ivanov scaling law yields a best fit of 𝜒 6 "X = 1.93 for SFD #8, which is within 1 σ of the other two scaling laws (Table 4). It produces f values near 15 for smaller projectiles, though, and this leads to a crater retention age of Tast =1.17 [-0.17, +0.16] Ga, much younger than from the other scaling laws. This value is outside the range estimated from dynamical models for the age of the Koronis family, though see the caveats in Sec. 3.5.1.

Schmedemann et al. (2014) applied the Ivanov scaling law model and a main belt SFD

similar to SFD #8 to crater counts found over a limited region of Ida. They assumed the Pi value for the impactors was 3.6 × 10 -18 km -2 yr -1 , while their impact velocities were 3.3 km s -1 . Both values are about 90% of our values in Table 3. In this region, crater spatial densities were found to be roughly a factor of 2 higher than those in Chapman et al. (1996b) (Fig. 12). Their best fit to these craters, most of which were between 1 and 2 km in diameter, yielded an age of 3.4 to 3.6 Ga.

Other surfaces with fresher craters suggested ages of ~2.1 Ga when their model was fit to data from craters between 0.5 and 1 km in diameter.

Summary (Ida)

The main belt fit scaling law with f ~10, along with the HH scaling law, are our preferred solutions, with main belt SFD solutions within high numbers favored from our fits. They produce mean crater retention ages of Ida near 2.5-2.9 Ga. These values are consistent with the modelderived dynamical age of the Koronis family, with Ida as a family member. The saturated crater regions on Ida do exist, however, and they could suggest an older age for Ida than 3 Ga. This would be tolerable given errors on existing family and dynamical constraints.

Gaspra

(951) Gaspra was the first asteroid ever observed by spacecraft (Fig. 11). The data returned by NASA's Galileo flyby in 1991 provided scientists with their first glimpse of what an S-type asteroid look likes up close. Gaspra has dimensions of 18.2 km × 10.5 km × 8.9 km and a mean diameter of 12.2 km [START_REF] Thomas | The shape of Gaspra[END_REF]. It is a prominent member of the Flora family that dominates the innermost region of the main asteroid belt [START_REF] Nesvorný | Identification and dynamical properties of asteroid families[END_REF]. The proper orbital elements of Gaspra are (a, e, i) = (2.21 au, 0.15, 5.1°).

Constraints on the age of Gaspra can be inferred from the evolution of the Flora family, which was formed from the catastrophic disruption of an asteroid of Dast > 150 km in the orbital region adjacent to the ν6 secular resonance [START_REF] Durda | Size-frequency distributions of fragments from SPH/ N-body simulations of asteroid impacts: Comparison with observed asteroid families[END_REF][START_REF] Vokrouhlický | Forming the Flora family: Implications for the near-Earth asteroid population and large terrestrial planet impactors[END_REF]. [START_REF] Granvik | Super-catastrophic disruption of asteroids at small perihelion distances[END_REF]2018). We will discuss Itokawa in more detail below, but the correspondence of these ages suggests that Gaspra should likely have a crater retention age comparable to 1.3 ± 0.3 Ga. In our modeling work below, we assign Gaspra a bulk density of 2.7 g cm -3 .

The population of Gaspra's craters has been reported and modeled by several groups (e.g., [START_REF] Belton | Galileo encounter with 951 Gaspra: First pictures of an asteroid[END_REF][START_REF] Greenberg | Collisional history of Gaspra[END_REF]Chapman et al., 1996a;[START_REF] Chapman | Cratering on asteroids from Galileo and NEAR Shoemaker[END_REF][START_REF] O'brien | Craters on asteroids: Reconciling diverse impact records with a common impacting population[END_REF][START_REF] Marchi | Cratering on asteroids[END_REF]. Numerous craters were identified between 0.16 < Dcrater < 1.9 km that followed a cumulative power law slope of -2.6. This size limit means Gaspra does not sample the knee in the crater SFDs observed near Dcrater ~ 2 km on Vesta, Lutetia, Ida, and others. A possible exception may be the mysterious facets on Gaspra, one or more of which may be ancient craters with Dcrater > 2 km (e.g., [START_REF] Greenberg | Collisional history of Gaspra[END_REF][START_REF] Thomas | The shape of Gaspra[END_REF][START_REF] O'brien | Craters on asteroids: Reconciling diverse impact records with a common impacting population[END_REF].

None of these pseudo-craters has been verified, so we do not include them in our analysis.

More recently, Gaspra's crater SFD was re-assessed using the Small Body Mapping Tool [START_REF] Runyon | Gaspra's craters: Implications for production functions and surface processes[END_REF][START_REF] Ernst | The Small Body Mapping Tool (SBMT) for Accessing, Visualizing, and Analyzing Spacecraft Data in Three Dimensions[END_REF]. The SBMT allows images to be wrapped onto an asteroid shape model, which is helpful for calculating surface areas when the body is irregular.

They counted 712 craters of 0.05 < Dcrater < 1.3 km within an area of 119.6 km 2 (Fig. 13). A rollover in their counts occurs for Dcrater < 0.17 km, which they attribute to limitations of image resolution.

Overall, the shape of their crater SFD was similar to that of Chapman et al. (1996a), with the power law slope of -2.6 reproduced. The normalization of the counts, however, was lower, probably because SBMT can more easily derive the irregular area of Gaspra's observed surface. In our modeling work, we use these new counts for Gaspra.

PLACE FIGURE 13 HERE

Empirical scaling law derived by fitting model and observed crater SFDs (Gaspra)

Assuming Pi = 2.635 × 10 -18 km -2 yr -1 (Table 3), we obtain a best fit for SFD #1 with 𝜒 6 BC = 1.94 (Table 4). Although this fit appears to favor lower-number SFDs, the other SFDs all have 𝜒 6 BC within 1s of this best fit case. The reason for this behavior is that Gaspra's observed crater SFD is effectively a power law, which is relatively easy for most models to fit. All of the main belt SFDs from Fig. 1 show power law slopes of q ~ -2.6 cumulative for their smallest craters, which matches the slope of Gaspra's craters (Chapman et al. 1996a).

Our best fit crater retention age for SFD #1 was Tast = 0.74 [-1.71, +0.41] Ga, nearly a factor of 2 lower than the Gaspra's expected age from constraints. Negative ages when the error bars are included are not meant to be taken literally. Given the similarity of our results for the different SFDs, the fact that the other asteroids favor high-number SFDs, and that Table 1 favors SFDs #5 and #6, we find it interesting that our model results for SFDs #5 and #6 yield mean ages between 1.18 < Tast < 1.38 Ga, all very close to our constrained age for Gaspra of 1.35 ± 0.3 Ga.

For these latter runs, our preferred value of f is 10.

Housen and Holsapple crater scaling fit (Gaspra)

Using the Holsapple and Housen (2007) scaling law and the same input parameters as Ida, expect for gravity, we obtain a best fit of 𝜒 6 B) = 1.93 for SFD #1. As before, all SFDs are within 1s of the best fit case, with SFDs #5 and #6 yielding mean ages of 1.12 < Tast < 1.32 Ga and a spread of 1.01 < Tast < 1.46 Ga when errors are included. The f value of the scaling law is near identical to the main belt fit's estimate of f = 10.

Ivanov scaling law fit (Gaspra)

The Ivanov scaling law provides the best fit of the three cases, with 𝜒 6 "X = 1.85 for SFD #1 (Table 4). Here SFDs #2-#8 are within 1 σ of the best fit. Its higher f values for smaller projectiles, however, yield an age for SFDs #5 and 6 between 0.25 < Tast < 0.29 Ga. These values are outside our age range derived from Flora and Itokawa constraints. [START_REF] Schmedemann | The cratering record, chronology and surface ages of (4) Vesta in comparison to smaller asteroids and the ages of HED meteorites[END_REF] applied the Ivanov scaling law model and a main belt SFD similar to SFD #8 to crater counts found over a limited region of Gaspra. They assumed a Pi value for their impactors of 3.54 × 10 -18 km -2 yr -1 , and their impact velocities were 4.69 km s -1 . Both values are modestly higher than our values in Table 3. They identified a steep crater SFD similar to that reported in Chapman et al. (1996a) that yielded a crater retention age of 0.27 ± 0.068 Ga.

This value is similar to our prediction for the Ivanov scaling law of 0.25 < Tast < 0.29 Ga.

Summary (Gaspra)

All of our best fit models favor SFD #1, but the crater retention ages produced by our fits are substantially different than age constraints from Itokawa samples and our estimated dynamical age of the Flora family. We also argue that SFD #1 is largely disfavored according to the crater SFDs of other main belt asteroids. The higher-number SFDs, however, produce results within 1s of this best fit case, with our preferred SFDs, namely #5 or #6, producing ages that are a good match to the Gaspra's additional age constraints (approximately 1.3 Ga). If we use those runs, our main belt scaling law fit results yield f ~ 10, the same values as for the asteroids discussed above.

Eros

We are now ready to consider the crater retention age of (433) Eros, the largest NEO observed by spacecraft. It was the primary target of NASA's NEAR mission and has a record of large craters that share commonalities with those of the main belt asteroids investigated above.

Eros has dimensions of 34.4 km × 11.2 km × 11.2 km and a bulk density of 2.67 ± 0.03 g cm -3 [START_REF] Thomas | Eros: Shape, topography, and slope processes[END_REF]) (Fig. 14). It is classified as an S-type asteroid, and its spectral characteristics suggest that it is similar to L or LL-type ordinary chondrites (e.g., [START_REF] Trombka | The elemental composition of asteroid 433 Eros: Results of the NEAR-Shoemaker X-ray spectrometer[END_REF][START_REF] Foley | Minor element evidence that Asteroid 433 Eros is a space-weathered ordinary chondrite parent body[END_REF][START_REF] Dunn | Mineralogies and source regions of near-Earth asteroids[END_REF][START_REF] Peplowski | Hydrogen and major element concentrations on 433 Eros: Evidence for an L-or LL-chondrite-like surface composition[END_REF][START_REF] Peplowski | The global elemental composition of 433 Eros: First results from the NEAR gamma-ray spectrometer orbital dataset[END_REF]. Recent spectral modeling work by [START_REF] Binzel | Compositional distributions and evolutionary processes for the near-Earth object population: Results from the MIT-Hawaii Near-Earth Object Spectroscopic Survey (MITHNEOS)[END_REF] agrees with this assessment; they find probabilities of 2%, 24%, and 74% that Eros is an H, L, and LL chondrite, respectively. Their interpretation is that Eros is most likely a LL chondrite.

PLACE FIGURE 14 HERE

All of these characteristics make Eros something of a mini-Ida. Ida is an S-type with a similar bulk density and shape, but it is only half the size (Sec. 3.5). A difference between the two is Eros's orbit; it currently crosses the orbit of Mars with (a, e, i) = (1.458 au, 0.223, 10.83°). This classifies it as an Amor-type NEO; it cannot currently strike the Earth, but its unstable orbit may put it in position to do so in the future [START_REF] Michel | The orbital evolution of the asteroid Eros and implications for collision with the Earth[END_REF].

Objects as large as Eros have little Yarkovsky mobility (e.g., Bottke et al., 2006a). For Eros-like bodies to escape the main belt and reach an Eros-like orbit, it is helpful if they are initially located in the innermost main belt region. Here a forest of overlapping Mars and three-body resonances creates a diffusive environment for asteroids (e.g., [START_REF] Morbidelli | Numerous weak resonances drive asteroids toward terrestrial planets orbits[END_REF][START_REF] Nesvorný | The Flora family: a case of the dynamically dispersed collisional swarm?[END_REF][START_REF] Nesvorný | Dynamical origin and terrestrial impact flux of large near-Earth asteroids[END_REF]. Alternatively, they need to have been created by a large asteroid disruption event occurring on the brink of a prominent resonance (e.g., [START_REF] Zappalà | Maria's family: Physical structure and possible implications for the origin of giant NEAs[END_REF]. Either way, dynamical models suggest that Eros likely spent hundreds of millions of years to many billions of years to escape the main belt. Only the last few millions to tens of millions of Eros's lifetime have been spent on planet-crossing orbits collisionally decoupled from the main belt population (e.g., [START_REF] Bottke | Origin of the Spacewatch small Earth-approaching asteroids[END_REF]2002). Based on this, we predict that most craters on Eros were derived from main belt impactors, and that we can model Eros's cratering history in the same manner as the main belt asteroids discussed above.

According to the NEO model of [START_REF] Granvik | Debiased orbit and absolute-magnitude distributions for near-Earth objects[END_REF], Eros was derived from the Hungaria region (44%), the innermost region of the main belt (47%), or the 3:1 resonance (9%). Given that the Hungarias currently have a paucity of Eros-sized S-type asteroids, we can probably reject that region as a source for Eros. Doing so increases the probabilities that Eros came from the innermost main belt and the 3:1 resonance to 84% and 16%, respectively. These results match those of [START_REF] Bottke | Debiased orbital and size distributions of the near-Earth objects[END_REF], whose NEO model did not include the Hungarias; they found that Eros has an 80% chance of coming from the inner main belt and 20% from the 3:1 resonance.

Given that Eros likely has an LL chondrite-like composition, it is natural to once again consider the Flora family as a possible source, especially given our results for Gaspra in Sec. 3.6.

Some have also postulated that Eros may be derived from the Maria asteroid family, which disrupted on the brink of the 3:1 resonance with Jupiter [START_REF] Zappalà | Maria's family: Physical structure and possible implications for the origin of giant NEAs[END_REF][START_REF] Zappalà | A search for the collisional parent bodies of large NEAs[END_REF].

Many regions on Eros appear to be close to crater saturation, with only the largest craters escaping this fate [START_REF] Chapman | Cratering on asteroids from Galileo and NEAR Shoemaker[END_REF][START_REF] Robinson | The geology of 433 Eros[END_REF]. The full database of Eros craters was provided to us by P. Thomas, who did the original mapping of Eros with M. Berthoud. The database is only likely to be complete for Dcrater > 0.2 km (P. Thomas, personal communication).

We estimate that the onset of saturation equilibrium takes place for Dcrater < 0.6 km, as determined when the power law slope of the craters at small sizes moves to a q = -2 cumulative index value.

For this reason, as with Ida, we will only examine the craters of Dcrater > 0.6 km (Fig. 15).

PLACE FIGURE 15 HERE

As an aside, we point out that seismic shaking appears to have reduced the spatial density of craters with 0.2 < Dcrater < 0.5 km from the vicinity of the 7.6-km Shoemaker crater [START_REF] Thomas | Seismic resurfacing by a single impact on the asteroid 433 Eros[END_REF], and that some mechanism-perhaps impact-induced seismic shaking-also erased craters with Dcrater < 0.1 km [START_REF] Richardson | Impact-Induced Seismic Activity on Asteroid 433 Eros: A Surface Modification Process[END_REF]). The issue of small crater erasure on asteroids is a fascinating one, but work on this topic is beyond the scope of this paper.

Empirical scaling law derived by fitting model and observed crater SFDs (Eros)

Given our limited information on the origin of Eros, we assume for now that Eros's original orbit was in the Flora asteroid family within the inner main belt. Accordingly, we assign it a starting orbit similar to Gaspra and give it a Pi value of 2.635 × 10 -18 km -2 yr -1 (Table 3). Our results indicate that we obtain our best fit using SFD #2, yielding 𝜒 6 BC = 3.72 (Table 4). All of the other SFDs also fit within 1s of this best fit case. Our best fit f value was 10.90 [-3.00, +3.00],

and our best fit crater retention age is Tast = 2.03 [-2.01, +0.86] Ga. If we were to instead adopt SFD #5 and #6, these mean values would change to f = 11.90 to 12.20 and Tast = 2.32 to 2.55 Ga.

Housen and Holsapple crater scaling fit (Eros)

Using the HH scaling law, a target density of 2.67 g cm -3 , and the same input parameters as for Gaspra and Ida, except for gravity, we obtain a best fit of 𝜒 6 B) = 3.85 for SFD #1. The SFDs #2-#7 are within 1s of the best fit case. Collectively, they yield mean ages between 2.1 < Tast < 4.5 Ga. The typical f values of this scaling law are close to f ~ 10 for the projectile sizes used here (Fig. 15). For SFDs #5 and #6, the mean ages are between 3.34 < 𝑇 :/3 < 3.81 Ga.

Ivanov scaling law fit (Eros)

The Ivanov scaling law provides the best fit of the three cases, with 𝜒 6 "X = 2.27 for SFD #8 (Table 4). Its higher f values yield an age of Tast = 1.4 [-0.19, +0.18] Ga (Table 5).

Interpretation of Eros's crater record

Our best fit mean crater retention ages from the main belt fit and the Marchi scaling law results are between 2.3 and 3.8 Ga. These values are substantially older than the estimated age of the Flora family (i.e., 1.35 ± 0.3 Ga; Sec. 3.6). One could dispute this, given that the best fit Ivanov scaling law's crater retention age is 𝑇 :/3 = 1.4 [-0.19, +0.18] Ga. The counterarguments are that Eros has notably higher crater spatial densities than Gaspra for similar-sized craters and Gaspra is a Flora family member (Figs. 13,15). Accordingly, we rule out Eros as a Flora family member on this basis.

It has also been postulated that Eros could come from the Maria family, located adjacent to the J3:1 resonance at high inclinations [START_REF] Zappalà | Maria's family: Physical structure and possible implications for the origin of giant NEAs[END_REF][START_REF] Zappalà | A search for the collisional parent bodies of large NEAs[END_REF]. As reported above, the odds of asteroids from the J3:1 reaching an Eros-like orbit are ~20%, less likely than the inner main belt but not unreasonably low. Several published dynamical ages for the Maria family have suggested that it is ~2 Ga (e.g., [START_REF] Spoto | Asteroid family ages[END_REF][START_REF] Aljbaae | The Maria Asteroid family[END_REF] or possibly 3 ± 1 Ga [START_REF] Brož | Constraining the cometary flux through the asteroid belt during the Late Heavy Bombardment[END_REF]. The collision probability of main belt asteroids with (170) Maria is Pi = 2.923 × 10 -18 km - 2 yr -1 , about 1.1 times the value used for the above age calculation. Multiplying our ages by this factor gives 2.2-4.1 Ga. These crater retention ages are in the same ballpark as the dynamical age estimates for Maria, given uncertainties.

The problem is that the [START_REF] Binzel | Compositional distributions and evolutionary processes for the near-Earth object population: Results from the MIT-Hawaii Near-Earth Object Spectroscopic Survey (MITHNEOS)[END_REF] spectral model predicts that (170) Maria has probabilities of 75%, 23%, and 2% of being an H, L, and LL chondrite, respectively. Given that Eros is likely a LL chondrite, it would appear that Maria can be ruled out as a candidate family on the basis of its spectral signature.

Accordingly, at this time, we have no independent constraints on the age of Eros. Given the comparable qualities of the fits in Fig. 15, we cannot use its crater SFD as a measure of which crater scaling law is preferred.

With that said, our crater studies above indicate that the empirical main belt and HH scaling law fits are preferred over the Ivanov scaling law fit. This suggests that the most likely scenario is that Eros was formed from the breakup of an asteroid in the inner main belt ~2.3 Ga to 3.8 Ga ago. The family has yet to be identified. If Eros came from a parent body that disrupted in the inner main belt, the reason that the family has not been found is plausibly because it disrupted in a highly diffusive region of this zone. Given the complicated network of resonances that exist in the inner main belt, we find it interesting but perhaps not surprising that an Eros precursor might avoid leaving behind clues to its existence after such a long time period. Less likely but still possible are that Eros came from the breakup of a body near the 3:1 resonance or from the Hungaria asteroid region. Additional work on this issue is warranted. First, these NEOs have experienced an erasure process that eliminated craters smaller than many tens of meters to possibly up to hundreds of meters. We do not examine erasure mechanisms in this paper, but it has been suggested that impact-induced seismic shaking or perhaps regolith mobility driven by thermal cycling may be responsible for this deficit (e.g., [START_REF] Richardson | Impact-Induced Seismic Activity on Asteroid 433 Eros: A Surface Modification Process[END_REF][START_REF] Marchi | Cratering on asteroids[END_REF]. Accordingly, only the largest craters on NEOs may stretch back to deep time.

To make progress in our work below, we avoid craters that may have been affected by crater erasure mechanisms. We define this as the crater diameter size range on the SFD plots discussed below where the observed and model crater SFDs diverge from one another. For Itokawa, Toutatis, Bennu, and Ryugu, this occurs for crater diameters smaller than ~100 m, ~250 m, ~150 m, and ~70-150 m, respectively.

Second, the dynamical histories of NEOs are uncertain. If we do not know where these asteroids came from in the main belt, it is difficult to predict how their collision probabilities and impact velocities varied with time. To make progress in our work below, we estimate the likely source regions and dynamical pathways followed by the bodies using our knowledge of main belt families, asteroid spectral signatures, and asteroid dynamics.

Third, as discussed in the introduction, the crater scaling laws that should be employed on asteroids that are a few hundred meters to several kilometers in diameter are uncertain. For the work below, we focus our attention on two possibilities:

• Option 1. The crater retention ages of small asteroids, based on their largest craters, are long.

For NEOs, they potentially provide a record of each world's traverse from their main belt starting orbit to the resonance that pushed them onto a planet-crossing orbit (e.g., [START_REF] Bottke | Debiased orbital and size distributions of the near-Earth objects[END_REF]Bottke et al., , 2006a)). In some cases, they may even reach back to the time of their parent body's disruption event.

• Option 2. The crater retention ages of small asteroids are short. For NEOs, the largest craters may only tell us about the last part of their journey in the main belt and perhaps the few millions to tens of millions of years it took to obtain their current orbits.

Impact experiments and numerical hydrocode modeling work have suggested that f values for small asteroids could be much higher than those derived from our Sec. 3 work, with values of 20, 40, or even 100 possible, depending on the physical properties of the target (e.g., Tatsumi and Sugita 2018).

As an example, consider the recent numerical hydrocode impact experiments of [START_REF] Davison | Collisional Histories of Small Planetesimals: Implications for Constraining the Age of the Surface of OSIRIS-REx Target, Bennu[END_REF]. They created a Bennu-like target with a basalt equation of state that was 20% porous and hit it with a 0.7 m diameter projectile at 7 km s -1 . When the strength of the target was set to 0.1, 1, 10, and 100 kPa, the simulations yielded crater diameters of 37, 22, 14, and 8.9 m, respectively, which translates into f values of 53, 32, 20, and 13, respectively. The unresolved issue is what strength value is appropriate to modeling the largest craters on Bennu.

If the low strength values are correct, Option 2 is preferred and our spacecraft-observed NEOs should have young crater retention ages. Option 2 might even be the expected outcome, given the estimated short timescales needed to spin small asteroids up to mass shedding via YORP.

On the other hand, our modeling results above indicate that f ~ 10 values can explain most main belt crater SFDs on Dast > 10 km asteroids (with the possible exception of Dcrater < 2-3 km craters on Ceres). This suggests that Option 1 may be viable.

It is interesting to consider that craters on Vesta (Dast ~ 530 km) and Gaspra (Dast ~ 12 km) can be fit with f ~ 10 for Dcrater > 0.1 to 0.2 km, even though they are very different in size (i.e., the ratio of the diameter of Vesta to that of Gaspra is 44). The diameter ratio between Gaspra and our spacecraft-observed NEOs, which are 0.3 < Dast < 2.5 km, ranges from 5 to 40. This raises the possibility that the largest NEOs in our Sec. 4 sample might follow the same trend.

With all of these issues in mind, we start our investigation with Itokawa.

Itokawa

(25143) Itokawa, an Apollo-type S (IV)-type NEO, was the target of JAXA's Hayabusa mission (Fig. 16). Itokawa appears to be a rubble-pile asteroid, with an elongated shape (dimensions of 0.535 km × 0.294 km × 0.209 km) and an estimated bulk density of 1.9 ± 0.13 g cm -3 [START_REF] Fujiwara | The rubble-pile asteroid Itokawa as observed by Hayabusa[END_REF]. It has two main components covered with rocks and boulders, and it seems likely that it was reassembled from debris produced by the disruption of the Itokawa parent body.

PLACE FIGURE 16 HERE

Hayabusa rendezvoused with Itokawa in 2005 and returned numerous grains from this body to Earth in 2010. An analysis of Itokawa samples showed this asteroid has a composition similar to LL ordinary chondrites, results that were consistent with its spectroscopic signature [START_REF] Nakamura | Itokawa dust particles: A direct link between S-type asteroids and ordinary chondrites[END_REF]. Three grains yield 40 Ar/ 39 Ar shock degassing ages of 1.3 ± 0.3 Ga [START_REF] Park | 40 Ar/ 39 Ar age of material returned from asteroid 25143 Itokawa[END_REF], whereas seven phosphate grains dated using the U-Pb system provide reset ages of 1.51 ± 0.85 Ga [START_REF] Terada | Thermal and impact histories of 25143 Itokawa recorded in Hayabusa particles[END_REF]. The interpretation is that a large, possibly catastrophic impact event affected the Itokawa precursor 1.3-1.4 Ga ago.

Dynamically, it can be shown that Itokawa, which currently resides on a fairly Earth-like orbit with (a, e, i) = (1.324 au, 0.280, 1.621°), most likely came from the innermost region of the main belt. The NEO models of [START_REF] Bottke | Debiased orbital and size distributions of the near-Earth objects[END_REF] and [START_REF] Granvik | Debiased orbit and absolute-magnitude distributions for near-Earth objects[END_REF] suggest that the odds of Itokawa having this provenance are ~86-100%.

The most prominent S-type family in the inner main belt is Flora, which formed from the catastrophic collision of a parent body larger than 150 km in diameter (e.g., [START_REF] Durda | Size-frequency distributions of fragments from SPH/ N-body simulations of asteroid impacts: Comparison with observed asteroid families[END_REF].

As discussed in Sec. 3.6, dynamical models suggest that the age of the Flora family is 1.35 ± 0.3 Ga [START_REF] Vokrouhlický | Forming the Flora family: Implications for the near-Earth asteroid population and large terrestrial planet impactors[END_REF]; see also [START_REF] Dykhuis | Defining the Flora family: Orbital properties, reflectance properties and age[END_REF]. This age is consistent with the sample reset ages of the Itokawa samples and the inferred crater retention age of Gaspra from Sec.

3.6. Finally, Itokawa is an excellent match with the LL-like spectral signature of Flora itself (e.g., [START_REF] Reddy | Chelyabinsk meteorite explains unusual spectral properties of Baptistina Asteroid Family[END_REF][START_REF] Binzel | Compositional distributions and evolutionary processes for the near-Earth object population: Results from the MIT-Hawaii Near-Earth Object Spectroscopic Survey (MITHNEOS)[END_REF].

Putting these clues together, we predict that Itokawa was produced by the disruption of the Flora parent body, and that it was once a member of the Flora asteroid family. A likely evolution scenario is that after formation, Itokawa resided near ~2.2 au until it drifted inward far enough by Yarkovsky thermal forces to escape from the main belt via the ν6 secular resonance (perhaps near 2.14 au).

We can use this concept to calculate Itokawa's mean collision probability and impact velocity with the rest of the main belt population. In Fig. 17, we show Pi values (and report mean

Vimp values) for Itokawa model asteroids encountering the 682 main belt asteroids with Dast ≥ 50 km discussed in Sec. 2.6. The model asteroids were assigned semimajor axes of 2.14-2.2 au and proper eccentricity and inclination values similar to Gaspra's, a Flora family member (Table 3).

They yield a mean value of Pi = 2.401 (± 0.11) × 10 -18 km -2 yr -1 and Vimp = 4.98 (± 0.04) km s -1 .

For reference, Gaspra's Pi value of 2.635 × 10 -18 km -2 yr -1 is 1.1 times that of Itokawa's. We do not model the portion of Itokawa's orbit where it was collisionally decoupled from the main asteroid belt; dynamical models suggest that this portion of its evolution was short compared to its journey within the main belt region (e.g., Bottke et al., 2015b).

PLACE FIGURE 17 HERE

Our prediction is that if Option 1 is correct, the largest craters in Itokawa's crater SFDs, which are plotted per square kilometer, should lie on top of our model fit to Gaspra's crater SFD in Fig. 13, provided we compensate for Itokawa's slightly smaller net collision probability value (i.e., the ratio of collision probabilities for Gaspra and Itokawa is 1.1, so we need to multiply 𝑁 RU/=M0:3<0 (> 𝐷 . ) for Itokawa by this value). Using craters derived from the work of Hirata et al. ( 2009) (see also [START_REF] Marchi | Cratering on asteroids[END_REF], we tested this idea in Fig. 18. The plot shows that Itokawa's craters with Dcrater ~ 0.1 km-those that presumably are least likely to have been affected by crater erasure-appear to be an extension of the crater SFD found on Gaspra. This result fulfills the predictions of Option 1, though it does not prove it. Given the available information at this time, it can only be considered an interesting coincidence.

PLACE FIGURE 18 HERE

For the moment, let us assume that Option 1 is true. Doing so knocks down several additional "logical dominos" and forces us to make several additional predictions that can help us interpret the other NEOs observed by spacecraft:

• The crater retention age derived from the largest craters on Itokawa is the same age as the Flora family-forming event that took place approximately 1.35 ± 0.3 Ga ago.

• Despite Itokawa's small size, a crater scaling law of f ~ 10 allows us to reasonably estimate its crater retention age.

• Given that our preferred scaling law appears to work reasonably well for asteroids larger than 10 km with a variety of sizes and taxonomic types, and now appears to work for the largest craters on the smallest NEO investigated in our sample (i.e., ~0.3 km; Itokawa), it seems plausible that it will work in similar ways on NEOs such as Toutatis, Bennu, and Ryugu.

• Despite the fact that YORP spin-up timescales are thought to be fast on Itokawa and other small worlds (see, e.g., the detected spin-up strength in [START_REF] Lowry | The internal structure of asteroid (25143) Itokawa as revealed by detection of YORP spin-up[END_REF][START_REF] Vokrouhlický | The Yarkovsky and YORP effects[END_REF] for other small near-Earth asteroids), the effects of YORP spin-up apparently did not lead to the erasure of craters of Dcrater ~ 0.1 km.

At the moment, though, all we have is an interesting coincidence, and Option 2 must still be considered viable. With these ideas in mind, we move to Toutatis, which is much closer in size to Gaspra than Itokawa.

Toutatis

(4179) Toutatis is an S-type NEO that is currently residing within the 3:1 mean motion resonance with Jupiter (Fig. 19). Its osculating orbit of (a, e, i) = (2.53 au, 0.63, 0.45°) places it on an Earth-crossing orbit. The Chinese mission Chang'e-2 flew by Toutatis in 2012 and reported dimensions of 4.354 km × 1.835 km × 2.216 km [START_REF] Bu | New insights of asteroid 4179 Toutatis using China Chang'e-2 close flyby optical measurements[END_REF]; see also [START_REF] Huang | The Ginger-shaped asteroid 4179 Toutatis: New observations from a successful flyby of Chang'e-2[END_REF].

These values largely confirmed estimates made from shape models derived using radar data [START_REF] Ostro | Radar images of asteroid 4179 Toutatis[END_REF][START_REF] Ostro | Asteroid 4179 Toutatis: 1996 radar observations[END_REF][START_REF] Hudson | High-resolution model of asteroid 4179 Toutatis[END_REF][START_REF] Takahashi | Spin state and moment of inertia characterization of 4179 Toutatis[END_REF]. The estimated bulk density of Toutatis is between 2.1 and 2.5 g cm -3 [START_REF] Ostro | Asteroid 4179 Toutatis: 1996 radar observations[END_REF][START_REF] Birlan | Dynamic and physical considerations on the asteroids density[END_REF]. This value is between Itokawa's bulk density (1.9 ± 0.13 g cm -3 ) and that of Eros (2.67 ± 0.03 g cm -3 ). This could mean that Toutatis has an internal structure somewhere between a classical rubble-pile (possibly like Itokawa) and a fractured or possibly shattered object with coherent fragments (i.e., [START_REF] Richardson | Gravitational aggregates: Evidence and evolution[END_REF].

PLACE FIGURE 19 HERE

The source of Toutatis is unknown. Given the available information, it is certainly possible that Toutatis is simply a background object that recently escaped the main belt. If we assume that Toutatis was once part of a prominent family, however, we can deduce a plausible parent from the available clues.

We start with spectra. [START_REF] Reddy | Composition of near-Earth asteroid (4179)[END_REF] found that the spectral signature of Toutatis is most similar to undifferentiated L-chondrites, though it lies close to various boundaries between H and L chondrites. The model described by [START_REF] Binzel | Compositional distributions and evolutionary processes for the near-Earth object population: Results from the MIT-Hawaii Near-Earth Object Spectroscopic Survey (MITHNEOS)[END_REF] suggests a similar result, with probabilities of 46%, 46%, and 8% of Toutatis being an H, L, and LL chondrite, respectively. From these data, we can probably rule out LL chondrite source families in favor of H and L chondrite families.

Next, we consider the results of dynamical models. The NEO model of [START_REF] Granvik | Super-catastrophic disruption of asteroids at small perihelion distances[END_REF]2018) indicates that Toutatis has chances of 22%, 53%, 15%, and 9% of coming from the u6 secular resonance near 2.2 au, the 3:1 mean motion resonance with Jupiter, the 5:2 mean motion resonance with Jupiter, and being a Jupiter-family comet, respectively. There are no known Jupiter-family comets that look like S-type asteroids, so we can rule out that possibility at this time. Using the NEO model of [START_REF] Bottke | Debiased orbital and size distributions of the near-Earth objects[END_REF], we find chances of 11%, 18%, 33%, and 39%

of coming from the ν6 secular resonance; the intermediate source Mars region, mostly in the inner main belt; the 3:1 resonance; and outer main belt sources (e.g., 5:2 resonance). These two models mostly agree with one another, though the outer main belt source is substantially higher in [START_REF] Bottke | Debiased orbital and size distributions of the near-Earth objects[END_REF]. Typically, NEOs in these models are favored to come from the inner main belt, which dominates the production of all NEOs, so the higher values found for the 3:1 and 5:2 resonances are intriguing.

We are now ready to consider the observed craters on Toutatis's surface, which range from 40 < Dcrater < 530 m [START_REF] Huang | The Ginger-shaped asteroid 4179 Toutatis: New observations from a successful flyby of Chang'e-2[END_REF][START_REF] Jiang | Boulders on asteroid Toutatis as observed by Chang'e-2[END_REF] (Fig. 20). Craters smaller than a few hundred meters appear to have been depleted by some kind of crater erasure mechanism, like those seen on other small asteroids observed by spacecraft (e.g., Eros, Steins, Itokawa, Ryugu, Bennu; [START_REF] Marchi | Cratering on asteroids[END_REF]Sugita et al., 2019;[START_REF] Walsh | Craters, boulders and regolith of (101955) Bennu indicative of an old and dynamic surface[END_REF]. If we only consider the largest craters, we find that their spatial densities are comparable to those on Ida or Eros (Figs. 12,15).

Accordingly, depending on the assumed collision probability Pi, and assuming Option 1 is valid, the crater retention age of Toutatis based on the largest craters is likely to be on the order of several billion years.

PLACE FIGURE 20 HERE

This result helps our deductive process because there are not many prominent S-type families in the main belt that are old enough to be plausible sources for Toutatis. If we only consider age as a discriminant, the candidate families include Flora, Maria, Eunomia, and Koronis.

Flora can be eliminated because it has an LL-type composition. The Maria family is a possibility;

as discussed above, it is adjacent to the 3:1 resonance, it may be > 2 Ga old, and the Binzel et al.

(2019) model suggests that (170) Maria is most likely to resemble an H chondrite. The Eunomia family appears to be in the LL camp, though L chondrites cannot be ruled out (Verzazza et al., 2014). It is also adjacent to the 3:1 resonance. Finally, there is the Koronis family located next to the 5:2 resonance. The [START_REF] Binzel | Compositional distributions and evolutionary processes for the near-Earth object population: Results from the MIT-Hawaii Near-Earth Object Spectroscopic Survey (MITHNEOS)[END_REF] model indicates that (158) Koronis, like Toutatis, has comparable odds of being an H or L chrondrite. This apparent similarity could indicate that the Koronis family is a good spectral match with Toutatis.

Additional clues to the origin of Toutatis may come from its inclination, which is very low (0.45°). In general, asteroids evolving into the 3:1 resonance from high-inclination sources such as the Maria or Eunomia families (proper inclinations > 12°) have difficulty reaching such low inclinations. Instead, it is much easier to reach a low inclination orbit by starting with a low inclination source such as the Koronis family (proper inclinations ~1-3°). However, Koronis family members evolving into the powerful 5:2 resonance at 2.8 au are less likely to obtain a Toutatis-like semimajor axis unless they are perturbed by the Earth during an encounter. This result is reflected in the dynamical results above. Still, the fact that Toutatis is currently on the Earthcrossing line may be a hint that such an event or events took place.

To further quantify our dynamical arguments, we used results from the Granvik et al.

(2018) numerical runs to determine how often test asteroids from the Maria/Eunomia and Koronis families reach the (a, e, i) orbit of Toutatis. For the former, we examined the evolution of 1759 test asteroids that reached the 3:1 resonance from starting orbits of a > 2.5 au, 0.05 < e < 0.18, and 11° < i < 16°. These objects reached resonance in the simulation via Yarkovsky drift forces. We found that 17 of these bodies passed within Δa = 0.03 au, Δe = 0.03, and Δi = 1° of the (a, e, i) orbit of Toutatis (i.e., 2.53 au, 0.63, 0.45°).

For the Koronis family, we tracked the evolution of 274 test asteroids that entered into the 5:2 resonance from starting orbits of a > 2.5 au, 0.0 < e < 0.10, and 0° < i < 3°. Here 9 test asteroids met our threshold. If we assume that the strength of the two sources above was equal, and we normalize the results by the number of test asteroids used, we find that the ratio favoring Koronis as a source for Toutatis over Maria/Eunomia is 3.4 (i.e., 1759 / 274 × 9 / 17).

Finally, we return to our cratering results in Fig. 20 and compare craters on Toutatis to Koronis family member Ida. We find it appealing that Toutatis crater counts are a good match to an extrapolation of the model crater SFD from Ida, though we cannot rule out a coincidence.

Putting the evidence together, we postulate that on the basis of spectral, dynamical, and cratering evidence, Toutatis is a lost member of the Koronis family. If Option 1 is correct, the crater retention age of Toutatis based on its largest craters should be the same as Ida, approximately 2-3 Ga.

Bennu and Ryugu

For our final test, we examine the crater histories of the NEOs Bennu and Ryugu, the targets of the OSIRIS-REx and Hayabusa2 sample return missions, respectively. They share a number of similarities, so we discuss them in tandem. The constraints discussed below mainly come from • They both have top-like shapes. Ryugu's dimensions are 1.04 × 1.02 × 0.88 km, while those of Bennu are 0.506 × 0.492 × 0.457 km. This makes Ryugu about twice and eight times as large as Bennu from a diameter and volume perspective, respectively.

• Both asteroids have a composition similar to primitive carbonaceous chondrites (e.g., CM or CI chondrites). The bodies are spectrally different, but in modest ways: Bennu is taxonomically classified as a B-type asteroid, whereas Ryugu is considered a Cb-type asteroid. Bennu's composition is akin to aqueously altered CM-type carbonaceous chondrites, whereas Ryugu's spectral signature is consistent with thermally and/or shock-metamorphosed CMs.

• They have the same bulk densities and geometric albedos: 1.19 g cm -3 and ~4.5%, respectively.

• Both Bennu and Ryugu appear to be rubble-pile asteroids, with the definition given by [START_REF] Richardson | Gravitational aggregates: Evidence and evolution[END_REF]. Their surfaces are a jumble of rocks and boulders that were likely produced in the aftermath of a family-forming event.

• Their orbital parameters are similar to those of Earth and each other: Bennu's (a, e, i) is (1.126 au, 0.204, 6.035°), and Ryugu's is (1.190 au, 0.1902, 5.884°).

• The most likely candidate families to produce these bodies are Eulalia and New Polana, located in the inner main belt at low inclinations near the 3:1 mean motion resonance with Jupiter (J3:1) at 2.5 au (e.g., [START_REF] Campins | The origin of asteroid 101955[END_REF][START_REF] Walsh | Introducing the Eulalia and new Polana asteroid families: Re-assessing primitive asteroid families in the inner Main Belt[END_REF]Bottke et al., 2015b). Bottke et al. (2015b) argues that the New Polana and Eulalia families having an approximate ~70% and ~30% probability of producing Bennu and Ryugu, respectively. The dynamical ages of the New Polana and Eulalia families are modestly different from one another; the former is 1400 [+150, -150] Ma, whereas the latter is 830 [+370, -100] Ma (Bottke et al., 2015b).

• Both Bennu and Ryugu have obliquities that are nearly 180°; Bennu's is 177.6°, and Ryugu's is 171.64°. This orientation, a probable outcome of YORP evolution (e.g., Bottke et al., 2006a;[START_REF] Vokrouhlický | The Yarkovsky and YORP effects[END_REF], indicates that both objects were migrating inward via the Yarkovsky effect when they escaped the main belt (e.g., Bottke et al., 2015b).

• Given the location of the New Polana and Eulalia families between ~2.4-2.49 au, the strongest likelihood is that both Bennu and Ryugu drifted inward across the inner main belt from these starting orbits. A likely departure zone from the main belt was through the ν6 secular resonance that defines the innermost boundary of the inner main belt. At low inclinations, the ν6 resonance escape zone is near 2.15-2.2 au [START_REF] Bottke | Debiased orbital and size distributions of the near-Earth objects[END_REF]. From there, they reached their current orbits via planetary encounters and interactions with resonances. Given the short lifetime of most NEOs (e.g., a few million to a few tens of millions of years; [START_REF] Bottke | Debiased orbital and size distributions of the near-Earth objects[END_REF][START_REF] Granvik | Debiased orbit and absolute-magnitude distributions for near-Earth objects[END_REF], Bottke et al. (2015b) predicted both Bennu and Ryugu escaped the main belt relatively recently compared to their long transit across the inner main belt.

To glean insights into the crater retention ages of Bennu and Ryugu, it is useful to compare their crater SFDs to an asteroid whose age and crater history are arguably well constrained. Here we choose Gaspra, a member of the Flora family, as our reference surface (Sec. 3.6;Figs. 11 and 13). The reasons why are as follows: Gaspra is located near the likely escape route of Bennu and Ryugu from the main belt, and Gaspra's crater retention age of 1.3 ± 0.3 Ga is both close to the estimated family ages of Eulalia and New Polana and is arguably well defined (i.e., it is consistent with sample ages from Itokawa and the likely age of the Flora family; see Secs. 3.6 and 4.1). Here we will superpose the crater SFDs of Bennu and Ryugu on Gaspra's and examine the similarities and differences.

For the comparison to be meaningful, we will assume that the crater scaling laws for Bennu and Ryugu are the approximately the same as that of Gaspra. This may be incorrect, with the worlds having different diameters (0.5 and 1.0 km vs. 12.2 km) and compositions (primitive carbonaceous chondrite vs. LL chondrite), but for the moment we accept this premise. We return to this issue in Sec. 5.

Second, we need to scale the Bennu-Ryugu crater SFDs for the different collision probabilities that they experienced compared to Gaspra over their orbital histories. As discussed above, Bennu and Ryugu probably came from the low-inclination Eulalia or New Polana families.

Accordingly, Bennu and Ryugu probably started with semimajor axes of 2.4 < a < 2.48 au and inclinations i ~ 2-3°. Next, they would have slowly migrated inward across the inner main belt until escaping out of the ν6 resonance. To account for this evolution, we ran collision probability simulations of test asteroids with (a, e, i) = (2.14 to 2.48 au, 0.1, 3°) against our asteroid population with Dast > 50 km using the methodology discussed in [START_REF] Bottke | Velocity distribution among colliding asteroids[END_REF] (Fig. 17). We ignore the small portion of time that Bennu and Ryugu were on planet-crossing orbits. If we assume that Bennu and Ryugu had a starting orbit in the Eulalia family, the mean Pi value for their evolution is (3.3 ± 0.46) × 10 -18 km -2 yr -1 . If we instead assume that they started in the New Polana family at 2.4 au, the mean Pi for their evolution is (3.1 ± 0.39) × 10 -18 km -2 yr -1 .

The collision probability of Gaspra, Pi = 2.635 × 10 -18 km -2 yr -1 , is 0.8 to 0.85 times that of Bennu/Ryugu. If we want to compare the crater SFD of Gaspra to those of Bennu and Ryugu, we need to multiply 𝑁 RU/=M0:3<0 (> 𝐷 . ) for Bennu and Ryugu by one of these values. This will allow us to visually determine whether Bennu and Ryugu are younger or older than Gaspra in a relative sense. If Bennu and Ryugu's crater SFD is lower on a cumulative plot than Gaspra's, its crater retention age is younger than 1.3 ± 0.3 Ga, whereas if it is higher, it is older. Our results are shown in Figs. 21 and22.

PLACE FIGURE 21 HERE PLACE FIGURE 22 HERE

We find that Ryugu's craters with Dcrater > 0.15 km (Sugita et al. (2019) appear to be a good fit with Gaspra's craters and the model production crater population (Fig. 21). Taken at face value, this would suggest that Ryugu's retention age for these craters is in the range of 1.3 ± 0.3 Ga. This result is interesting because the dynamical ages of Ryugu's postulated source families, New Polana and Eulalia, also match these values (i.e., 1400 [+150, -150] Ma, 830 [+370, -100] Ma, respectively; Bottke et al. 2015b).

Bennu's largest craters [START_REF] Walsh | Craters, boulders and regolith of (101955) Bennu indicative of an old and dynamic surface[END_REF]) also appear to produce a comparable fit, though there is room for interpretation based on the SFD of craters with Dcrater > 0.05 km (Fig. 22). Many of Bennu's proposed craters in this size range have subdued topography, which makes it difficult to know whether all of them are valid. In addition, the crater counts from [START_REF] Walsh | Craters, boulders and regolith of (101955) Bennu indicative of an old and dynamic surface[END_REF] are based on data acquired early in the Bennu encounter. Since then, the asteroid has been imaged and the topography measured via lidar to much higher resolution. Additional work on this issue is needed as the crater population is updated using these more recent data.

Summary

Our results suggest that Option 1 may be valid, and that the crater retention ages of NEOs, based on their largest craters, could be surprisingly old. The largest craters on Itokawa, a possible member of the Flora family whose age is ~ 1.3 Ga, appear to line up with the crater SFD found on Gaspra, a confirmed member of the Flora family. The same can be said for Toutatis, a possible member of the Koronis family (~2 to 3 Ga), and Ida, a confirmed member of the Koronis family.

The largest craters on Bennu and Ryugu also appear to be as ancient as those found on Gaspra, once we account for the different collision probabilities of the impacting population. For each one, this could suggest that some aspects of their surfaces go back as far as the family-forming event that made them. It may be possible to check Option 1 by analyzing the samples returned by the OSIRIS-REx and Hayabusa2 missions.

In Sec. 7, we will discuss the implications of our findings and whether they make sense given what we know about other modeling and observational data.

Discussion

In this paper, we use a new formulation of the main belt size distribution to examine the crater histories of asteroids observed by spacecraft. Some of the key takeaways from our work are as follows.

1. We have used a disruption scaling law that allows asteroids of Dast ~ 0.2 km to break up more easily than with the scaling law used by Bottke et al. (2005b) (or Benz and[START_REF] Benz | Catastrophic disruptions revisited[END_REF]. It is capable of producing a main belt size distribution that is more consistent with crater constraints than previous work. We find our best results for 𝑄 ! * functions that are higher in number, with #5 and #6 favored when Table 1 probabilities are also considered.

A question that emerges from our research is whether our new asteroid disruption law reflects reality. Unfortunately, it is difficult to assess this issue here without obtaining additional main belt and asteroid constraints. The mismatch between our scaling law and those derived using hydrocode simulations, such as in [START_REF] Benz | Catastrophic disruptions revisited[END_REF], may be a clue that our methodology is missing something, previous scaling laws are missing something, or that everyone is missing something.

In addition to possible issues with methodologies, there may be an issue of how the problem we are investigating is framed. For example, it has been argued that the spin-up of small asteroids by the YORP thermal torques produce frequent mass shedding events, and that this extra source of disruption acts to make the power law slope of the main belt size distribution more shallow than expected between 0.2 < Dast < 2 km [START_REF] Marzari | Combined effect of YORP and collisions on the rotation rate of small Main Belt asteroids[END_REF][START_REF] Jacobson | Effect of rotational disruption on the size frequency distribution of the main belt asteroid population[END_REF][START_REF] Penco | Yarkovsky depletion and asteroid collisional evolution[END_REF]; see discussion in Bottke et al., 2015a). If true, the new asteroid scaling law would essentially replace YORP disruption with impact disruption. From the physics perspective, this works against our collisional evolution model, which does not include this effect, but from the frame of trying to model crater size distributions on asteroids, it may not matter if the two methods yield the same main belt SFD.

It is also conceivable that both processes work together. Laboratory experiments indicate that asteroids spinning near their rotational breakup limit are much easier disrupt by impact than slow-spinning bodies (e.g., Holsapple 2007). This effect was recently studied by Ševeček et al. (2019), who performed a large number of numerical impact simulations with rotating targets using a smoothed particle hydrodynamics code coupled to an N-body code (e.g., see Durda et al. 2004).

They found that the critical energy needed to disrupt a target (i.e., the 𝑄 ! * function) changed rapidly when one approaches the critical spin rate of an asteroid. Unfortunately, their study was limited to bodies 10 to 100 km in diameter, substantially larger than the ~0.2-km bodies whose disruption threshold changes the most in Fig. 2. Future work on this issue is needed.

For the moment, let us assume that the YORP effect drives a modest fraction of rubble-pile asteroids with 0.2 < Dast < 2 km to spin near their disruption limit (e.g., [START_REF] Pravec | Spin rate distribution of small asteroids[END_REF][START_REF] Pravec | Formation of asteroid pairs by rotational fission[END_REF], and that this makes them easier to break up by impacts. From a one-dimensional collisional modeling perspective, when their short collisional lifetimes are combined with longer ones from slow rotators, the net effect is that these asteroids are easier to disrupt on average than before. In effect, this rationale can explain the shape of our new asteroid disruption law, even if our purely collisional model does not account for the physics producing it.

New work will be needed to see how our asteroid disruption law holds up when new constraints become available, more is known about how asteroids disrupt via YORP spin-up, or more in known about the impacts and YORP spin-up working in tandem.

2. Our favored main belt size distribution (#5 or #6) can be successfully fit to the observed crater size distributions found on Ceres (Kerwan basin), Vesta (Rheasilvia basin), Lutetia, Mathilde, Ida, Gaspra, and Eros. Our results indicate that the ratio of crater diameter to projectile diameter f for craters with Dcrater > 0.1 km on these worlds is ~10. The only exception found so far may be craters smaller than a few km on the 930-km-diameter asteroid Ceres.

Our fits between model and crater data indicate that a simple relationship exists between projectiles and crater sizes on asteroids of Dast > 10 km. Our test set includes a wide range of asteroid sizes and compositions. Our results also appear broadly consistent with the crater scaling laws of Holsapple and Housen (2007), provided we use certain input parameters (i.e., those suggested by Marchi et al. (2012b) for stony asteroids).

With this said, we add some cautionary notes on this interpretation. The size distribution of craters smaller than a few km on Ceres appears to differ from those on Vesta, with the Ceres SFD having a much steeper power law slope. If this is not a byproduct of secondary craters, the easiest explanation is that these craters on Ceres are a byproduct of an increasing f value and possibly different material properties for the terrains in question.

It should also be said that our f ~10 solution for Mathilde's crater SFD is (i) not unique, (ii) based on limited crater data, and (iii) fit to craters that are close to saturation. This situation is unfortunate because Mathilde is the only carbonaceous chondrite-like asteroid imaged by spacecraft to date that is larger than Ryugu (0.9 km) and smaller than Ceres (930 km). Still, we consider our solution reasonable because f >> 10 scenarios would require some kind of surface reset event within the past billion years that so far lack supporting evidence. The fact that Mathilde has no observed family [START_REF] Nesvorný | Identification and dynamical properties of asteroid families[END_REF] is an argument against the occurrence of such a reset event.

The reasons why the f ~10 scaling law works as well as it does for many different asteroid sizes and compositions will require additional study, but certain factors probably play a role:

• Most of the asteroids investigated to date probably have comparable material strengths, at least against impact events. They can be considered fractured or shattered versions of the stony meteorites in our collection.

• The collision velocities between asteroids in the main belt does not vary strongly from world to world, with typically mean velocities near ~5 km s -1 [START_REF] Bottke | Velocity distribution among colliding asteroids[END_REF].

•

The surface histories of the asteroids investigated here have been dominated by impacts from a main belt SFD whose shape has been in quasi-steady state for billions of years. Although different asteroids disrupt over time, the broad-scale properties of the main belt population have not changed substantially from a collisional evolution perspective over that interval.

The influence of impactors embedded in the early main belt population (e.g., comets implanted in the primordial asteroid belt; [START_REF] Levison | Contamination of the asteroid belt by primordial trans-Neptunian objects[END_REF][START_REF] Vokrouhlický | Capture of trans-Neptunian planetesimals in the main asteroid belt[END_REF]) has yet to be detected in the crater histories of main belt asteroids. The reason is probably because few if any of the surfaces investigated in this paper go back to the primordial days of the main belt. To explore earlier bombardment eras, we would need to examine the most ancient surfaces on Vesta or other as-of-yet unobserved large asteroids. This remains a fascinating topic for future work.

3. Our derived main belt size distribution, combined with a crater scaling law of f ~10, can be fit to the largest craters on Itokawa, Toutatis, Bennu, and Ryugu. All of these asteroids are smaller than Gaspra. The match yields a crater retention age on the order of ~1 Ga for Itokawa, Bennu, and Ryugu and ~2.5 Ga for Toutatis. These values are consistent with the computed formation ages of their source families (Itokawa from the Flora family, Bennu and Ryugu from the Eulalia or New Polana families, Toutatis from the Koronis family), though that does not prove they actually have these surface ages.

The possibility that the largest craters on these asteroids are ancient matches modeling work (e.g., [START_REF] Walsh | Craters, boulders and regolith of (101955) Bennu indicative of an old and dynamic surface[END_REF] but we still consider it something of a surprise. Therefore, we will discuss this topic further below. Our work shows four out of four examples where the crater retention ages of small asteroids match their predicted family ages. If this is merely a coincidence, it is a good one. With that said, there are other factors to consider here, and they may suggest that the crater retention age of Itokawa, Toutatis, Bennu, and Ryugu are younger than postulated. We present the arguments and their possible counters below.

Factor 1. Does Spin Up from YORP Reset the Surfaces of Most Small NEOs?

Over the past decade, it has become increasingly apparent that the YORP thermal torques modify the spin rates and obliquities of asteroids smaller than Gaspra (e.g., reviewed in Bottke et al. 2002Bottke et al. , 2006a;;[START_REF] Vokrouhlický | The Yarkovsky and YORP effects[END_REF]. They provide the easiest way to explain the spin rate distribution of NEOs and MBAs-which include numerous bodies spinning near fission speeds and other bodies that have almost no rotational angular momentum (e.g., [START_REF] Pravec | Spin rate distribution of small asteroids[END_REF]and their obliquity distribution, with small MBAs having preferentially extreme obliquity values (e.g., [START_REF] Hanuš | Asteroids' physical models from combined dense and sparse photometry and scaling of the YORP effect by the observed obliquity distribution[END_REF][START_REF] Ďurech | Asteroid models reconstructed from the Lowell Photometric Database and WISE data[END_REF] and many NEOs having values near ~180° (e.g., La [START_REF] Spina | Retrograde spins of near-Earth asteroids from the Yarkovsky effect[END_REF][START_REF] Farnocchia | Near Earth Asteroids with measurable Yarkovsky effect[END_REF]. They may also provide a ready explanation for the spinning top-like shapes of Bennu, Ryugu, and many other asteroids (e.g., [START_REF] Walsh | Formation and evolution of binary asteroids[END_REF].

Direct measurement of YORP accelerations indicate that they should frequently cause asteroids to undergo mass shedding events (e.g., [START_REF] Pravec | Formation of asteroid pairs by rotational fission[END_REF][START_REF] Jewitt | The active asteroids[END_REF]). As a second example, in situ studies of Bennu indicate that its rotation rate is currently accelerating at a rate of (3.63 ± 0.52) × 10 -6 deg day -2 , enough to double Bennu's rotation rate in 1.5 Myr (e.g., [START_REF] Hergenrother | The operational environment and rotational acceleration of asteroid (101955) Bennu from OSIRIS-REx observations[END_REF][START_REF] Scheeres | The dynamic geophysical environment of (101955) Bennu based on OSIRIS-REx measurements[END_REF], and a similar acceleration was detected for a number of other small NEAs (reviewed in [START_REF] Vokrouhlický | The Yarkovsky and YORP effects[END_REF].

The implication of these results is that small asteroids can be quite dynamic places. Many should undergo surface changes on timescales that are short compared to our estimated crater retention ages for small NEOs.

A caveat is so-called stochastic YORP, namely that asteroid shape changes driven by mass movement, craters, and mass shedding can cause an asteroid's spin acceleration to undergo a random walk [START_REF] Statler | Extreme sensitivity of the YORP effect to small-scale topography[END_REF]Bottke et al., 2015b). This effect may prevent some asteroids from reaching the kinds of rotation speeds that allow for frequent mass shedding events, or at least considerably delay it, whereas others may enter into mass shedding events again and again.

Another intriguing possibility is that YORP self-regulates itself into a long-lived equilibrium states for the surviving population of bodies, or at least a sub-population, weakening considerably the YORP effects (e.g., [START_REF] Golubov | Systematic structure and sinks in the YORP effect[END_REF].

Until this process is better understood, we must leave open the possibility that Itokawa, Toutatis, Bennu, and Ryugu fall into the former class of objects, and that their largest craters are indeed ancient.

There is also the possibility that some top shapes are formed in the family-forming event,

with small objects growing by the gravitational reaccumulation of debris [START_REF] Michel | Disruption and Reaccumulation as the Possible Origin of Ryugu and Bennu Top Shapes[END_REF]. In this scenario, ejected material grows by the gravitational accretion of nearby bodies, and this leads to more mass being accreted near the equator of the growing rubble-pile asteroid than near the pole. If true, the shapes of Bennu and Ryugu may indeed be ancient, which would lend credence to the idea that they have old retention ages for their largest craters.

Factor 2. Are Near-Earth Objects Preferentially Long-Lived Asteroids?

Our CoDDEM model results allow us to make predictions of the typical collisional lifetime of asteroids based on a choice of our asteroid disruption law (𝑄 ! * ). We find that the sizes of the smallest asteroids in our sample, Itokawa and Bennu, are relatively close to the minimum in the 𝑄 ! * function shown in Fig. 2. Accordingly, these bodies should be easier to disrupt on an energy per mass scale by impact events than asteroids that are considerably smaller or larger.

There are feedbacks here, however, with the shallow slope of the main belt SFD between 0.2 and 2 km leading to fewer projectiles for Dast < 0.2 km than estimated by Bottke et al. (2005b).

This leads to an average collisional lifetime for Itokawa-and Bennu-sized asteroids on the order of one hundred to a few hundred Myr, comparable to the values estimated in Bottke et al. (2005b).

Regardless, these intervals are much shorter than the estimated crater retention ages of Itokawa and Bennu. This mismatch is a potential argument that the crater retention ages of the largest craters are not the order of 1 Ga, but instead are much younger.

The counter to this argument is to consider how Itokawa-and Bennu-sized asteroids from a given family escape the main belt after a family-forming event. As a useful example, we refer the reader to the model results from Bottke et al. (2015b). They used numerical simulations to track how Bennu-sized asteroids from the Eulalia, New Polana, and Erigone families evolved in semimajor axis by the coupled Yarkovsky/YORP effects. Their work accounted for the likely collisional lifetime of their model asteroids; those that disrupted were removed from the simulation. Their goal was to reproduce the distribution of these families in semimajor axisabsolute magnitude (a, H) space, where the observed families make a quasi-"V"-shape.

The setup for the Bottke et al. (2015b) simulations were as follows. They assumed a large breakup event created a size distribution of asteroids representing the Eulalia, New Polana, and Erigone families. As shorthand, we classify 0.3-to 1-km bodies as "small", 1-to 4-km bodies as "modest-sized", and objects larger than 4 km as "big". All of these bodies begin to drift inward or outward towards resonances via the Yarkovsky and YORP effects. YORP torques cause the obliquities of the bodies to evolve toward 0° or 180°, where they obtain their maximum semimajor axis drift velocities from the Yarkovsky effect. The smaller fragments migrate more quickly than the modest-sized ones, which in turn move faster than the big ones. This creates a V-shape in (a, H) space, with smaller objects evacuated from the middle of the family (e.g., Fig. 18 of Bottke et al. 2015b). We refer to the two sides of the V-shape as "ears".

As small and modest-sized bodies disrupt, according to their assumed collisional lifetimes, attrition takes its toll on the leading edge of each ear. Only a fraction live long enough to make it to an "escape hatch" resonance that will take them out of the main belt. Some of these bodies also go through YORP cycles, where their obliquity values are reset by various processes. The small bodies have the shortest timescale to undergo YORP cycles, and it causes them to undergo a random walk in semimajor axis, slowing their progress toward an escape hatch resonance.

Eventually, though, this wave of surviving bodies reaches an escape hatch resonance. A specific example is shown in Fig. 19 of Bottke et al. (2015b), where Bennu-sized bodies from

Eulalia and New Polana reached various inner main belt resonances over timescales of 0.1 to > 2

Gyr after the family forming event. The expected flux of these escaping bodies at the estimated ages of the families was found to be consistent with the available constraints.

Taken together, these results indicate that the NEO population might be dominated by "lucky" asteroids that are long-lived survivors. An analogy might be World War I soldiers running across no man's land to reach the trenches of their enemy; most soldiers fall during the assault, but a few make it. If this scenario is accurate, Itokawa, Toutatis, Bennu, and Ryugu would be survivors of this gauntlet, which would make their putative ancient large craters less surprising and more of a selection effect. On the other hand, a small asteroid can survive and still not have an ancient crater retention age. The takeaway is that it may not be straightforward to interpret the surfaces of these small asteroids.

It is also possible that Itokawa, Toutatis, Bennu, and Ryugu are second-generation family members, but we predict that is unlikely. When big family members disrupt, they will create some modest-sized bodies and lots of small bodies, all which can now drift more rapidly. The starting location of these fragments, however, will often be closer to the center of the family than at great distances from the center. The initial second-generation population will also be in smaller in numbers than the initial first-generation population, and they experience the same attrition factors as first-generation bodies. All of this suggests that the second-generation bodies that escape the main belt are unlikely to outnumber the first-generation bodies for a considerable time after the family-forming event.

Factor 3. What Crater Scaling Laws are Applicable to Small Asteroids?

There is considerable debate about the projectile sizes needed to make small craters on asteroids, primarily because the strength of the surface materials is unknown. Estimates in the literature differ by many orders of magnitude, and this can change f values from 10 to 40 or more (e.g., [START_REF] O'brien | Craters on asteroids: Reconciling diverse impact records with a common impacting population[END_REF][START_REF] Tatsumi | Cratering efficiency on coarse-grain targets: Implications for the dynamical evolution of asteroid 25143 Itokawa[END_REF]).

An innovative way to glean insights into this issue was through the Small Carry-on Impactor (SCI) experiment carried by JAXA's Hayabusa2 mission to Ryugu. The SCI consisted of a 30-cm disk impactor made of copper that was accelerated by an explosion to an impact speed with Ryugu of 2 km s -1 . The impact produced a cone-shaped debris curtain and a crater-like feature that was approximately 17.6 ± 0.7 m in diameter rim to rim [START_REF] Arakawa | An artificial impact on the asteroid 162173 Ryugu formed a crater in the gravity-dominated regime[END_REF]). This value is so large that the Hayabusa2 team suggested that Ryugu's surface acts like it has the same strength as cohesionless sand upon impact, corresponding to f >> 10. If so, and if we can assume this result is applicable to all Ryugu craters, it would imply that Ryugu's crater retention age is considerably younger than 1 Ga.

The caveat that we can provide at this time is that the creation of craters with Dcrater > 100 m on Ryugu may be different than the formation of much smaller craters a few tens of meters across. If true, different projectile sizes are sensitive to how the terrains can change as they become larger (e.g., boulder concentrations per unit area, etc.).

Factor 4. Do Scaling Laws Change the Size Distributions of Small Craters?

Assuming that the crater retention ages for Itokawa, Toutatis, Bennu, and Ryugu are young carries its own implications. For example, it means that crater scaling laws must change between Eros-and Gaspra-sized and larger objects (Dast > 10 km), which have f ~ 10, and smaller asteroids, where presumably f >> 10. From our work above, it is not yet clear how one would tell the difference between the two.

As an example, consider two asteroids with the same age, a large one where f ~ 10 and a small one where f ~ 50. Our model main belt SFD in Fig. 1 shows that the slope of < 0.1-km asteroids follows a Dohnanyi-like power law slope of -2.6 (e.g., O'Brien and [START_REF] O'brien | Steady-state size distributions for collisional populations: Analytical solution with size dependent strength[END_REF]Bottke et al., 2015a). If we assume that projectiles striking from this portion of the SFD will create a crater SFD on these two bodies with the same slope, the only outward difference between the two would be in crater spatial densities.

The hope would be to find a transitional target body where large craters follow f ~ 10 and smaller craters follow f >> 10. This would produce a change in the slope of the crater SFD that would be observable. The issue is that crater erasure on small bodies makes such signatures rare or hard to interpret. We speculate that crater erasure may even be a byproduct of this change in f; models indicate that increasing f for small craters might make them more effective at crater erasure.

This scenario warrants further work.

The prevalence of a Dohnanyi-like power law slope of -2.6 at small asteroid sizes may also help explain the coincidence of why the largest craters on Itokawa, Toutatis, Bennu, and Ryugu have the same crater retention age as their putative source families. For every target body, there must be a largest crater than can form; asteroids striking a surface that are larger than a critical threshold will produce an asteroid-wide surface resetting event and/or disruption. Our work

indicates that f values near 10 or >> 10 will be drawn from the same -2.6 slope in the main belt SFD. If f >> 10, we can expect a young crater retention age, with the largest possible crater forming relatively quickly compared to the formation age of the asteroid.

If this logic holds, it could be rare for the surface of a small asteroid to be caught between a reset event and the formation of its largest crater. This would make the true crater retention ages of Itokawa, Toutatis, Bennu, and Ryugu much younger than suggested.

One way to test this hypothesis is to determine the shock resetting (e.g., Ar-Ar) and surface (cosmic ray exposure) ages of samples returned from Bennu and Ryugu. We might expect to find the formation age of the family in the samples, as was attempted for Itokawa (e.g., [START_REF] Park | 40 Ar/ 39 Ar age of material returned from asteroid 25143 Itokawa[END_REF], but there may also be evidence of other impact events. If they skew toward younger ages, this may be evidence that impacts only affected those rocks that had been brought to the surface over relatively recent times. On the other hand, finding numerous impact ages could tell us that (i) the surface is indeed ancient or (ii) there have been a multitude of resetting events, with the samples recording impacts back to the formation of the target asteroid.

Future Work

Future work promises to increase the fidelity of the four components of our crater production models presented in Sec 2.1. For example, one way to improve our calculation of Component 1 (assessment of main belt SFD) is for new surveys of main belt asteroids to come online that can more easily detect sub-km asteroids. A present-day example of this comes from [START_REF] Heinze | The Flux Distribution and Sky Density of 25th Magnitude Main Belt Asteroids[END_REF], who discussed how Dark Energy Camera observations have probed the small body size distribution of the main belt. It is also expected that future wide field surveys such as the Vera C. Rubin Observatory, previously referred to as the Large Synoptic Survey Telescope (LSST), and NASA's Near-Earth Object Surveillance Mission will be able to detect numerous main belt bodies in the sub-km range across the main belt.

These surveys will also give us increased knowledge of the nature of the km-and sub-km main belt populations. This information is needed to calculate improved collision probabilities and impact velocities between our target asteroids and a representative population of main belt asteroids (Component 3).

In terms of crater scaling laws for small asteroids (Component 2), new data will be provided by NASA's Double Asteroid Redirection Test (DART) mission [START_REF] Cheng | AIDA DART asteroid deflection test: Planetary defense and science objectives Planet[END_REF]) and ESA's Hera mission [START_REF] Michel | The Hera Mission: European Component of the Asteroid Impact and Deflection Assessment (AIDA) Mission to a Binary Asteroid[END_REF]. DART is a kinetic impactor that will hit the 160 m moon of Didymos, a 780-m-diameter near-Earth asteroid. This event will make a crater and change the angular momentum of the system. Hera will follow-up the DART mission with a detailed post-impact survey of the Didymos system. From Hera, we will gain new insights into the projectile sizes needed to make small craters on asteroids, and this will translate into new constraints for crater scaling laws.

Finally, we predict that new dynamical modeling work and meteorite/sample analysis may lead to additional constraints for Component 4, the crater retention age of target asteroids. For example, future work may help us identify the likely source family of Eros. If that happens, we can use dynamical models to determine the age of the family-forming event and compare that age to the predicted crater retention age of Eros (i.e., as discussed in Sec. 5.7, we modestly favor a crater retention age of ~2.3 to 3.8 Ga, but that is based on an assumed starting location and evolutionary history for Eros). In addition, new work on the shock degassing ages of meteorites may allow us to pin down the starting location of certain asteroids and refine our predicted crater retention ages. We look forward to an analysis of the samples return from Bennu by the OSIRIS-REx spacecraft and from Ryugu by the Hayabusa2 spacecraft. They may retain a record of the family-forming event that created them.

Conclusions

Here we summarize the main conclusions of our paper.

We have derived a new formulation of the main belt SFD that is modestly different from the one proposed by Bottke et al. (2005b) (Sec. 3). We have successfully fit it to numerous asteroid crater SFDs observed by spacecraft. Our results for Ceres (multi-km craters only), Vesta, Lutetia, Mathilde, Ida, Eros, and Gaspra yield a crater scaling law where the ratio of crater sizes to projectile sizes is a factor f ~ 10 (Sec. 5).

From a probability standpoint, our results favor the main belt SFDs #5 and #6 (Fig. 1; see also Fig. 2). They are likely to reproduce main belt SFD and family constraints as well as constraints from asteroid craters (Sec. 5).

Our derived empirical scaling law largely match the results produced by the crater scaling law of Holsapple and Housen (2007), provided certain parameters are used (i.e., we assume that the surface material acts like cohesive soils, and that the yield strength of most asteroids is on the order of Y ~ 2 ×10 7 dynes cm -2 . Typical values from our input parameters yield f ~ 10 (Sec. 3).

Conversely, the scaling law of Ivanov (2004) yields poorer fits than the other tested scaling laws for large asteroids, mainly because its parameter choices yield f >> 10 for target asteroids that are Dast > 10 km. The f >> 10 values also produce crater retention ages on asteroids that are inconsistent with independent ages derived using dynamical methods and/or sample evidence (Sec.

5).

For the spacecraft-observed asteroids and surfaces tested in this paper, we summarize our findings below and in Table 6.

PLACE TABLE 6 HERE

Vesta (Rheasilvia Basin). We find that the crater retention age of the ~500-km Rheasilvia basin on Vesta is probably younger than 1.3 Ga, with a plausible value near ~1 Ga. These ages are consistent with the 40 Ar/ 39 Ar ages found in feldspar grains taken from eucrite meteorites, which range between 0.6 and 1.7 Ga (Lindsey et al., 2015) (Sec. 5.1).

Ceres (Kerwan Basin). We predict that the crater retention age of the ~280-km Kerwan basin on Ceres, based on fits to the multi-km and larger craters, is approximately 0.8 to 0.9 Ga.

Our preferred model of f ~ 10, however, may not fit the crater SFD for Dcrater < 2 km. The inconsistency in the power law slopes of sub-km craters on Ceres and Vesta could suggest that (i) secondary craters dominate the crater SFD of small craters on Ceres or (ii) the physical properties of Ceres's surface are so different from that of Vesta that they allow f >> 10 for small craters (Sec.

5.2).

Lutetia (Achaia Region). The crater retention age of the oldest surface on M-type asteroid Lutetia (Achaia) is the order of 2.5 to 3.5 Ga, though error bars indicate it could also be as old as the formation of Lutetia itself (Sec. 5.3).

Mathilde. The crater retention age of Mathilde, a carbonaceous chondrite asteroid that is larger than Bennu and Ryugu and smaller than Ceres, appears to be somewhere between 2.2 Ga and the formation age of the body itself, which could go back to the planetesimal formation era.

Our interpretation is that f ~ 10 for this body makes the most sense; f >> 10 would yield crater retention ages so young that we would expect to see a family associated with Mathilde (Sec. 5.4).

Ida. The crater retention age of S-type asteroid Ida, a member of the Koronis family, appears to be between 2-3 Ga, though it could be older depending on the cratered surface analyzed.

This age is consistent with the expected dynamical age of the Koronis family (Sec 5.5).

Gaspra. The crater retention age of S-type asteroid Gaspra, a member of the Flora family, appears to be ~ 1.3 Ga. This age is consistent with the expected dynamical age of the Flora family.

It also matches the 40 Ar/ 39 Ar ages of LL chondrite grains returned from Itokawa by the Hayabusa spacecraft: 1.3 ± 0.3 Ga [START_REF] Park | 40 Ar/ 39 Ar age of material returned from asteroid 25143 Itokawa[END_REF]. We predict that Itokawa is very likely to have been a

Flora family member prior to becoming an NEO (Sec. 5.6).

Eros. The crater retention age of S-type asteroid Eros is between ~2.3 to 3.8 Ga. Although

Eros is an NEO with a likely LL chondrite composition, it does not appear to be a lost member of the Flora family. The origin of Eros is currently an unsolved problem (Sec. 5.7).

Itokawa. Itokawa is a small S-type NEO that was likely to have once part of the Flora family (see Gaspra above). If we limit our analysis to its largest craters, we find that Itokawa's crater spatial densities match those found on Gaspra, another Flora family member. Moreover, if these large craters formed with f ~ 10, they suggest that components of Itokawa's crater history tell the story of the billion-year interval between the formation of the Flora family and the present day (Sec. 6.1). On the other hand, the crater retention age based on the largest craters could be young if f >> 10.

Toutatis. Based on the results of dynamical models and spectroscopic interpretation, we predict that the 2.5-km S-type NEO Toutatis was once part of the Koronis family. We find that the crater spatial densities of its largest craters match trends in the SFD of Ida, a Koronis family member. If these large craters formed with f ~ 10, components of the observed surface on Toutatis could be as old as the Koronis family-forming event (Sec. 6.2).

Bennu and Ryugu. Bennu and Ryugu are 0.5-and 0.9-km carbonaceous chondrite asteroids, respectively. The crater spatial densities of their largest craters were found to be comparable to the crater SFD found on Gaspra, whose surface is the order of 1.3 Ga old. If these large craters formed with f ~ 10, both worlds could have observed surface features as old as the predicted source families for these worlds, Eulalia and New Polana, which are nearly 1 Ga old (Sec. 6.3).

The largest craters on Itokawa, Toutatis, Bennu, and Ryugu yield ages of ~1, 2.5, 1, and 1

Ga, respectively, provided f ~ 10. There are some reasons to think it could be illusionary, with the reality being that the largest craters on these worlds formed with f >> 10, but we are unable to dismiss these ancient ages out of hand.

If the largest craters on these worlds are indeed old, all four of these NEOs, and in fact most NEOs, are probably likely the fortunate survivors of collisional evolution within the main belt region. Tracking the collisional and dynamical evolution of individual main belt bodies across the main belt, and accounting for the possibility of a collisional cascade that allows breakup event to produce daughter fragments, can potentially be simulated by the next generation of numerical models (Sec. 7).

On the other hand, if the largest craters on these four NEOs are instead relatively young, f either had to go from a value of 10 for Eros/Gaspra/larger bodies to f >> 10 for only slightly smaller bodies or it had to change on the surfaces of the four NEOs themselves. In the latter case, the evidence for this change would be found by modeling the SFD of the smaller craters that have been strongly affected by crater erasure. Finding evidence of how and why this happens could lead us to a better understanding of the physical nature of asteroids and how they are affected by collisions (Sec. 7). 1). The results are given for the three different crater scaling laws discussed in Sec. 2.5. The second column is M, the number of crater data points on the asteroid/region in question (see Eq. 14). * #1, and so on). SFD #1 was designed to match the one used in Bottke et al. (2005b). The blue dots represent the debiased main belt SFD as discussed in Bottke et al. (2005a, b). For reference, the red line shows the NEO SFD as defined by [START_REF] Harris | The population of near-Earth asteroids[END_REF]. As shown in Fig. 2, the higher index numbers correspond to lower minimum values for 𝑄 ! * . This allows more asteroids to disrupt between 0.1 < Dast < 1 km, which in turn means the power law slope of the SFD becomes shallower in that range. The upturn in slope occurs at larger sizes for low index numbers (~0.5 km for #1) and smaller sizes for high index numbers (~0.2 km for #8). 3). The predicted ages for Kerwan basin are ~0.8-0.9 Ga from the empirical main belt fit and the HH scaling law. The Ivanov scaling law predicts an age of ~0.5 Ga, smaller than the other two scaling laws because its f values in the inset figure are much higher. Plot components are as in Fig. 4. The best fit model crater SFD is the empirical main belt fit and the HH scaling law (Table 3), both of which have f values near 9 or 10 for the majority of observed craters. The predicted mean age of the Achaia region from the empirical main belt fit and HH scaling law model is ~2.5-3.5 Ga. 3). They both yield f values near 10-12 for the majority of observed craters and crater retention ages of ~2 to > 4 Ga. Many other SFDs fit the data within 1σ of the best fit cases (Table 3). If we apply higher-number SFDs to make our model crater SFD, the crater retention ages derived from the empirical main belt fit and the HH scaling law approach the age of the Solar System.

Figure 11. Ida and Gaspra to the same scale. Gaspra (right) was imaged at a range of 5,300 km, and Ida (left) was imaged from 3,000 to 3,800 km, both by the Galileo spacecraft. Gaspra's dimensions are 18.2 km × 10.5 km × 8.9 km, and Ida's are 59.8 km × 25.4 km × 18.6 km [START_REF] Belton | Galileo encounter with 951 Gaspra: First pictures of an asteroid[END_REF]1996). Both bodies are S-type asteroids that are members of asteroid families: Gaspra is part of the Flora family, which may be 1.3 Ga old, and Ida is part of the Koronis family, which is 2-3 Ga old [START_REF] Bottke | Dynamical spreading of asteroid families via the Yarkovsky effect: The Koronis family and beyond[END_REF][START_REF] Vokrouhlický | The vector alignments of asteroid spins by thermal torques[END_REF]2017;[START_REF] Brož | Constraining the cometary flux through the asteroid belt during the Late Heavy Bombardment[END_REF][START_REF] Nesvorný | Identification and dynamical properties of asteroid families[END_REF][START_REF] Spoto | Asteroid family ages[END_REF]. Courtesy NASA/JPL/USGS.

Figure 12. A comparison between the observed crater SFDs found on S-type asteroid Ida, a member of the Koronis family, and various crater models. The observed crater counts are from Fig. 5 of [START_REF] Chapman | Cratering on Gaspra[END_REF]. Plot components are as in Fig. 4. The best fit model crater SFD is the empirical main belt fit, but both the HH scaling law and the Ivanov scaling law yield comparable fits (Table 3). The empirical main belt fit scaling law and the HH scaling law are our preferred solutions, with the inset showing that both have f ~ 9-11.

Figure 13. A comparison between the observed crater SFDs found on S-type asteroid Gaspra, a member of the Flora family, and various crater models. The observed crater counts are from [START_REF] Runyon | Gaspra's craters: Implications for production functions and surface processes[END_REF]. Plot components are as in Fig. 4. The HH crater scaling law fit in red overlaps the main belt empirical fit in blue. All of our best fit models favor SFD #1 (Table 3), but all yield crater retention ages that are substantially lower than the age constraints from Itokawa samples and our estimated dynamical age of the Flora family. We prefer SFD #5 or #6, which yield results within 1σ of the best fit case, yet yield mean crater retention ages between ~1.1 and ~1.4 Ga, values that are close to Gaspra's estimated age of 1.35 ± 0.3 Ga from additional constraints (Fig. 18). For these latter runs, our preferred value of f is ~10. [START_REF] Robinson | The geology of 433 Eros[END_REF]. Craters that are in probable saturation (i.e., those with Dcrater < 0.6 km) are not shown. Plot components are as in Fig. 4. The HH crater scaling law fit in red largely overlaps the main belt empirical fit in blue. The best fit model crater SFD here is the Ivanov scaling law fit, though all of the fits are fairly comparable (Table 3). See text for discussion of the crater retention age of Eros and possible source families.

Figure 16. Three near-Earth asteroids observed by spacecraft that are smaller than 1 km. Ryugu (left) is a Cb-type asteroid with a mean diameter of 896 m and a bulk density of 1.19 ± 0.02 g cm -3 [START_REF] Watanabe | Hayabusa2 arrives at the carbonaceous asteroid 162173 Ryugu --A spinning top-shaped rubble pile[END_REF]. Its spectral signatures are consistent with thermally and/or shock-metamorphosed carbonaceous chondrite meteorites (Kitazato et al., 2019). Bennu (center) is a B-type asteroid with a mean diameter of 492 m and a bulk density of 1.190 ± 0.013 g cm -3 [START_REF] Lauretta | The unexpected surface of asteroid (101955) Bennu[END_REF]. Its composition is similar to aqueously altered CMtype carbonaceous chondrites. Itokawa (right) is an S-type asteroid with dimensions of 0.535 km × 0.294 km × 0.209 km and an estimated bulk density of 1.9 ± 0.13 g cm -3 [START_REF] Fujiwara | The rubble-pile asteroid Itokawa as observed by Hayabusa[END_REF]. Grains from Itokawa indicate that it has the composition of an LL-type ordinary chondrite. Images courtesy of NASA/JAXA.

Figure 17. The model collision probabilities of Itokawa, Bennu, and Ryugu as they evolved inward toward the Sun across the main belt by the Yarkovsky effect. For Itokawa (red), we assumed that it started at 2.2 au, the center of the Flora family, and then migrated to 2.14 au where it escaped the main belt via the υ6 secular resonance. Our Itokawa model asteroids were assigned proper eccentricity and inclination values similar to Gaspra (Table 2). The collision probabilities Pi were calculated with 682 main belt asteroids with Dast > 50 km (Sec. 2.6). The mean values of all red points are Pi = 2.401 (± 0.11) × 10 -18 km -2 yr -1 and Vimp = 4.98 (± 0.04) km s -1 . For Bennu and Ryugu (blue), we assumed that they started at 2.487 au and 2.4 au, the centers of the Eulalia and New Polana families, respectively. The model asteroids were given eccentricities and inclinations of 0.1 and 3°, respectively. The rest of the method was same as with Itokawa. The mean Pi value for Bennu and Ryugu starting in the either Eulalia and New Polana families was (3.3 ± 0.46) × 10 -18 km -2 yr -1 and (3.1 ± 0.39) × 10 -18 km -2 yr -1 , respectively, and their mean impact velocities Vimp were 4.72 (± 0.10) km s -1 and 4.67 (± 0.12) km s -1 , respectively. [START_REF] Bu | New insights of asteroid 4179 Toutatis using China Chang'e-2 close flyby optical measurements[END_REF]; see also [START_REF] Huang | The Ginger-shaped asteroid 4179 Toutatis: New observations from a successful flyby of Chang'e-2[END_REF]. It was imaged by the Chang'e-2 spacecraft during a flyby that had a closest approach distance of 770 ± 120 meters [START_REF] Huang | The Ginger-shaped asteroid 4179 Toutatis: New observations from a successful flyby of Chang'e-2[END_REF]. Courtesy CNSA.

Figure 20.

A comparison between craters found on Toutatis and Ida and various crater models. Plot components are as in Fig. 4. Ida is a member of the Koronis family, and its crater counts are the black open circles (Fig. 12). Its estimated crater retention age is ~2.5 Ga.

Toutatis is currently a NEO, but its orbit and spectral signature can be plausibly linked to the Koronis asteroid family. The largest craters on Toutatis [START_REF] Huang | The Ginger-shaped asteroid 4179 Toutatis: New observations from a successful flyby of Chang'e-2[END_REF][START_REF] Jiang | Boulders on asteroid Toutatis as observed by Chang'e-2[END_REF] green circles) are aligned with those of Ida's crater SFD. If the largest craters on Toutatis were made with the same crater scaling law as Ida, they could represent the same crater retention age. The smaller craters on Itokawa have been strongly affected by a crater erasure mechanism (e.g., [START_REF] Marchi | Cratering on asteroids[END_REF]. The empirical main belt fit curve (blue line) represents a mean crater retention age of ~2.4 Ga, whereas the HH scaling fit is ~2.9 Ga (red line).

6

  Near-Earth Asteroids Smaller than 10 Kilometers Armed with the insights gleaned from Sec. 3, we now consider the crater retention ages of smaller NEOs observed by spacecraft: (25143) Itokawa, (4179) Toutatis, (101955) Bennu, and (162173) Ryugu. It is challenging to model the cratering history of any of these bodies for several reasons.

  Lauretta et al. (2019), Hamilton et al. (2019), Sugita et al. (2019), and Watanabe et al. (2019):

Figure 1 .

 1 Figure 1. Collisional evolution model results for the main asteroid belt, based on the assumptions and model results of Bottke et al. (2005b). The model SFDs #1-8 are assigned an index number corresponding to the asteroid disruption laws 𝑄 ! * #1-8 shown in Fig. 2 (i.e., SFD #1 was produced by 𝑄 ! * #1, and so on). SFD #1 was designed to match the one used in Bottke et al. (2005b). The blue dots represent the debiased main belt SFD as discussed in Bottke et al.(2005a, b). For reference, the red line shows the NEO SFD as defined by[START_REF] Harris | The population of near-Earth asteroids[END_REF]. As shown in Fig.2, the higher index numbers correspond to lower minimum values for 𝑄 ! * . This allows more asteroids to disrupt between 0.1 < Dast < 1 km, which in turn means the power law slope of the SFD becomes shallower in that range. The upturn in slope occurs at larger sizes for low index numbers (~0.5 km for #1) and smaller sizes for high index numbers (~0.2 km for #8).

Figure 2 .Figure 3 .

 23 Figure 2. The asteroid disruption laws used in our collisional evolution model runs. Each disruption law 𝑄 ! * is assigned an index number #1-8, and they produce the model main belt SFDs shown in Fig. 1 (i.e., SFD #1 was produced by 𝑄 !* #1, and so on). The parameters needed to generate the curves can be found in Table1. Disruption law 𝑄 ! * #1 matches the one used inBottke et al. (2005b). The green dot is a normalization point determined from laboratory impact experiments. It is defined as (𝑄 ! * D)C , 𝐷 D)C ) = (1.5 ×10 7 erg g -1 , 8 cm) (e.g.,[START_REF] Durda | Collisional models and scaling laws: A new interpretation of the shape of the main-belt asteroid size distribution[END_REF] The minimum 𝑄 ! * value (𝑄 !

Figure 4 .

 4 Figure 4. A comparison between the observed crater SFDs found on the floor and ejecta blanket of Vesta's Rheasilvia basin and various crater models. The observed crater counts are from Marchi et al. (2015). The best fit model crater SFDs for three different crater scaling laws are shown with the colored lines: the main belt empirical fit in blue (Sec. 3.1.1), the HH crater scaling law fit in red (Sec. 3.1.2), and the Ivanov crater scaling law in green (Sec. 3.1.3). The numbers in parentheses correspond to the index number of the model main belt SFD applied to produce the model crater SFD (Fig. 1). See Sec. 3.1.1 for discussion of error bars. The inset figure shows the ratio f = Dcrater / Dast for the different crater scaling laws as a function of the impacting asteroid's diameter. The best fit model crater SFDs is the HH scaling fit, though the main belt empirical fit matches everything but the very largest craters. The Ivanov scaling fit produces a model crater SFD that is lower than the observed data for Dcrater > 1.5 km.

Figure 5 .

 5 Figure 5. The largest crater on C-type asteroid Ceres is Kerwan basin. It is ∼284 km in diameter(Williams et al., 2017) and has a relaxed polygonal shape. The center of the basin is located at 10.8° south latitude and 123.9° east longitude. The image was taken during Dawn's Survey phase from an altitude of 4,400 km. Courtesy NASA/JPL-Caltech.

Figure 6 .

 6 Figure 6. A comparison between the observed crater SFDs found on or near the Kerwan basin on Ceres and various crater models. The observed crater counts are from Hiesinger et al. (2016) with updates from this study. Plot components are as in Fig.4. The best fit model crater SFD is the empirical main belt fit (Table3). The predicted ages for Kerwan basin are ~0.8-0.9 Ga from the empirical main belt fit and the HH scaling law. The Ivanov scaling law predicts an age of ~0.5 Ga, smaller than the other two scaling laws because its f values in the inset figure are much higher.

Figure 7 .

 7 Figure 7. ESA's Rosetta mission observed the M-type asteroid Lutetia during a flyby. The colors represent different regions of Lutetia as defined by geologic mapping. The oldest part of the asteroid, and the one investigated in the paper, is the heavily cratered Achaia region. This image was published in Thomas et al. (2012), who adapted it from Massironi et al., (2012). Copyright Elesvier.

Figure 8 .

 8 Figure 8. A comparison between the observed crater SFDs found on the Achaia region of Lutetia and various crater models. The observed crater counts are from Marchi et al. (2012b).Plot components are as in Fig.4. The best fit model crater SFD is the empirical main belt fit and the HH scaling law (Table3), both of which have f values near 9 or 10 for the majority of observed craters. The predicted mean age of the Achaia region from the empirical main belt fit and HH scaling law model is ~2.5-3.5 Ga.

Figure 9 .

 9 Figure 9. Mathilde is a large C-type main belt asteroid. The part of the asteroid shown is 59 by 47 km across. It was imaged by the NEAR spacecraft from a distance of 2,400 km. The surface exhibits many large craters, some which are partially shadowed. Courtesy NASA/JPL/JHUAPL.

Figure 10 .

 10 Figure 10. A comparison between the observed crater SFDs found on C-type asteroid Mathilde and various crater models. The observed crater counts are from O'Brien et al. (2005), who reformulated them from[START_REF] Chapman | Cratering on Gaspra[END_REF]. Plot components are as in Fig.4. The HH crater scaling law fit in red largerly overlaps the main belt empirical fit in blue. The best fit model crater SFDs are the empirical main belt fit and the HH scaling law (Table3). They both yield f values near 10-12 for the majority of observed craters and crater retention ages of ~2 to > 4 Ga. Many other SFDs fit the data within 1σ of the best fit cases (Table3). If we apply higher-number SFDs to make our model crater SFD, the crater retention ages derived from the empirical main belt fit and the HH scaling law approach the age of the Solar System.

Figure 14 .

 14 Figure 14. Eros is an S-type near-Earth asteroid. It has dimensions of 34.4 km × 11.2 km × 11.2 km (Thomas et al. 2002). This view of Eros's northern hemisphere is a mosaic of six NEAR spacecraft images taken from an orbital altitude of about 200 km. Psyche crater (5.3 km across) is located at the 12 o'clock position in the middle of the saddle-shaped region, and Himeros crater (11 km) can be seen on the opposite side of Eros at the 5 o'clock position. Courtesy NASA/JPL/JHUAPL.

Figure 15 .

 15 Figure15. A comparison between the observed crater SFDs found on S-type asteroid Eros, a near-Earth asteroid, and in various crater models. The observed crater counts are from the database of P. Thomas and M. Robinson[START_REF] Robinson | The geology of 433 Eros[END_REF]. Craters that are in probable saturation (i.e., those with Dcrater < 0.6 km) are not shown. Plot components are as in Fig.4. The HH crater scaling law fit in red largely overlaps the main belt empirical fit in blue. The best fit model crater SFD here is the Ivanov scaling law fit, though all of the fits are fairly comparable (Table3). See text for discussion of the crater retention age of Eros and possible source families.

Figure 18 .

 18 Figure18. A comparison between craters found on Itokawa and Gaspra and various crater models. Plot components are as in Fig.4. Gaspra is a member of the Flora family, and its crater counts are the black open circles (Fig.13). Its estimated crater retention age is ~1.3 Ga. Itokawa is currently an NEO, but it was probably a member of the Flora family in the past. Itokawa's crater SFD(Hirata et al., 2009; green dots) has been multiplied by the ratio of the collisional probabilities between Gaspra and Itokawa (Fig.17and Sec. 4.1). Within errors, the crater spatial densities of Gaspra and the largest craters on Itokawa appear the same. If the largest craters on Itokawa were made with the same crater scaling law as Gaspra, they could represent the same crater retention age. The smaller craters on Itokawa have been strongly affected by a crater erasure mechanism (e.g.,[START_REF] Richardson | Impact-Induced Seismic Activity on Asteroid 433 Eros: A Surface Modification Process[END_REF]). The model crater SFDs are discussed in Fig.13. Here the empirical main belt fit curve (blue line) represents a crater retention age of ~1.38 Ga, whereas the HH scaling fit is ~1.32 Ga (red line).

Figure 19 .

 19 Figure19. Toutatis is an S-type near-Earth asteroid with dimensions of 4.354 km × 1.835 km × 2.216 km[START_REF] Bu | New insights of asteroid 4179 Toutatis using China Chang'e-2 close flyby optical measurements[END_REF]; see also[START_REF] Huang | The Ginger-shaped asteroid 4179 Toutatis: New observations from a successful flyby of Chang'e-2[END_REF]. It was imaged by the Chang'e-2 spacecraft during a flyby that had a closest approach distance of 770 ± 120 meters[START_REF] Huang | The Ginger-shaped asteroid 4179 Toutatis: New observations from a successful flyby of Chang'e-2[END_REF]). Courtesy CNSA.

Figure 21 .

 21 Figure21. A comparison between craters found on Ryugu and Gaspra and various crater models. Plot components are as in Fig.4. Gaspra is a member of the Flora family, and its crater counts are the black open circles (Fig.13). Its estimated crater retention age is ~1.3 Ga. Ryugu is currently an NEO, but it is thought to have come from either the New Polana or Eulalia families, with an estimated age of 1.4 [+0.15, -0.15] Ga or 0.83 [+0.37, --0.10] Ga, respectively(Bottke et al. 2015b). Ryugu's crater SFD (Sugita et al. 2019) has been multiplied by the ratio of the collisional probabilities between Gaspra and Ryugu (Fig.17and Sec. 4.3). Within errors, the crater spatial densities of Gaspra and the largest craters on Ryugu are similar. If the largest craters on Ryugu were made with the same crater scaling law as Gaspra, they could represent a crater retention age near ~1.3 Ga. The smaller craters on Ryugu have been strongly affected by a crater erasure mechanism (e.g.,[START_REF] Marchi | Cratering on asteroids[END_REF]. The model crater SFDs are discussed in Fig.13and 18.

Figure 22 .

 22 Figure22. A comparison between craters found on Bennu and Gaspra and various crater models. Plot components are as in Fig.4. Gaspra is a member of the Flora family, and its crater counts are the black open circles (Fig.13). Its estimated crater retention age is ~1.3 Ga. Bennu is currently an NEO, but it is thought to have come from the New Polana and Eulalia families, with estimated ages of 1.4 [+0.15, --0.15] Ga and 0.83 [+0.37, -0.10] Ga, respectively(Bottke et al. 2015b). Bennu's crater SFD (Sugita et al. 2019) has been multiplied by the ratio of the collisional probabilities between Gaspra and Bennu (Fig.17and Sec. 4.3). Within errors, the crater spatial densities of Gaspra and the largest craters on Bennu are similar. If the largest craters on Bennu were made with the same crater scaling law as Gaspra, they could represent a crater retention age near ~1.3 Ga. The smaller craters on Bennu have been strongly affected by a crater erasure mechanism (e.g.,[START_REF] Marchi | Cratering on asteroids[END_REF]. The model crater SFDs are discussed in Fig.13and 18.

Table 3

 3 .90, with 1s of -1.40 and +1.00, whereas our best fit crater retention age for Rheasilvia was Tast = 0.85 Ga, with 1s errors of -0.23 and +0.24 Ga (Table5). The visual fit to the crater data in Fig.4is good except for Dcrater > 20 km craters. It is possible that the mismatch in this range stems from small number statistics.

	PLACE TABLE 5 HERE

, our best fit case was found for SFD #5. It yielded 𝜒 6 BC = 22.32 (Table

4

). The SFDs #4-8 yielded values within 1s of this best fit case. Our best fit f value was 9

2 Housen and Holsapple crater scaling fit (Rheasilvia)

  (Eqs. 10 and 11). Our input parameters for this equation were given in Sec. 3.1.2. The best fit is from SFD #8, but 𝜒 6 "X in this case ended up as 70.51, a value indicative of a poor fit. The reason is that this scaling law produces larger f values than the others tested for impactors of Dast < 0.1 to 0.2 km. To fit the smallest craters on Rheasilvia, the production function must substantially undershoot the craters with Dcrater > 1 km, as shown in Fig 4.The best fit crater retention age is Tast = 0.37 [-0.02, + 0.01] Ga, a value that is considerably younger than the previous two test cases. It falls outside the 40 Ar/ 39 Ar age range of the Kapoeta feldspar grains (0.6-1.7 Ga;Lindsey et al., 2015). It also does not match constraints on Vesta family's age from dynamics

  . It is located near the outer edge of the central

	main belt, with proper orbital elements of (a, e, i) = (2.77 au, 0.12, 9.6°). It is classified as a C-
	type asteroid, and observations from the DAWN spacecraft indicate that it is a volatile-rich rocky
	body. Studies based on DAWN spacecraft data have provided us with many critical parameters for
	Ceres, including its bulk density of 2.160 ± 0.009 g cm -3 and a surface gravity of 0.28 m s -2
	(Russell et al., 2016). The mineralogy and geochemistry of Ceres, as constrained by Dawn
	observations, appear consistent with the bulk composition of CM/CI carbonaceous chondrites
	(McSween et al., 2018).

Table 4

 4 The latter value is modestly higher than those in[START_REF] Hiesinger | Cratering on Ceres: Implications for its crust and evolution[END_REF], even though our collision probability Pi value is higher than their value of 2.84 × 10 -18 km -2 yr -1 . The reason is because SFD #8 has fewer small bodies than the SFD inBottke et al. (2005b) (SFD #1), which was used in their "asteroid-derived" model. Visually, the empirical main belt crater model reproduces the crater data reasonably well.

	). Both SFD

If we were to choose a more probable main belt SFD according to Table

1

, such as SFD #6, our best fit f value is 8.8 [-1.70, +1.20] and our best fit crater retention age 𝑇 :/3 is 0.86 [-0.18, + 0.18] Ga.

  >> 10, it would lead to a crater retention age for Mathilde that would be younger than ~1 Ga. Given the lack of evidence on Mathilde for any surface reset event,

SFD. In turn, those results suggest that Mathilde's crater retention ages go back to the earliest days of Solar System history. Arguments in favor of this interpretation come from Mathilde's lack of an asteroid family. Some might argue that the porous nature of C-type asteroids makes them less likely to produce families, but numerous C-type families have been identified across the main belt (e.g.,

[START_REF] Nesvorný | Identification and dynamical properties of asteroid families[END_REF][START_REF] Masiero | Asteroid family physical properties[END_REF]

, including those likely to have produced Bennu and Ryugu

(Bottke et al. 2015b

).

If our interpretation is valid, it obviates the need for crater erasure mechanisms that affect the observed craters smaller than several km (e.g.,

[START_REF] O'brien | Craters on asteroids: Reconciling diverse impact records with a common impacting population[END_REF]

. It would also argue that projectiles hitting carbonaceous chondrite-like materials do not necessarily lead to substantially larger craters than those made on other asteroids (e.g., Ceres craters; see Sec. 3.2). If Mathilde's crater scaling laws produced f

  Table1), the best fit f value moves to 10.00[-3.40, +3.00], and our best fit crater retention age is Tast = 2.36[-2.36, +2.37] Ga. The lower error bar is a formal number and

best fit f value is 10.90

[-2.70, +2.90]

, and our best fit crater retention age is Tast = 2.52

[-1.98, +0

.94] Ga. The f value is similar to our results in previous runs, not a surprise given the location of the knee in Ida's crater SFD at Dcrater ~ 2 km. If we apply a more probable main belt SFD, such as SFD #6 (

Table 1 .

 1 The parameters used to define the eight 𝑄 ! * asteroid disruption functions tested in this paper, as well as how they fared against constraints. The parameters in columns 2 to 4 are Dmin, defined as the location of the minimum 𝑄 ! * value (𝑄 ! * 5.E ), and the two variables 𝛼 and 𝛽, which are applied in Eq. 5, 6, and 7 to derive 𝑄 ! * . Columns 5 and 6 describe the number of trials out of 100 test runs that match both of our main belt constraints. The metric 𝜓 #(! 6 , defined by Eq. 4, describes how well the model main belt SFD compares to the observed main belt SFD. The metric

	Test #	Dmin (km)	𝛼		𝛽	% trials with 𝜓 #(! 6 and 𝜒 ()B 6 > 2σ < 20	6 % trials with 𝜓 #(! and 𝜒 ()B 6 > 2σ < 10
	1	0.2	-0.35	1.33		37	21
	2	0.2	-0.385	1.35		36	30
	3	0.2	-0.42	1.37		41	28
	4	0.2	-0.455 1.405		47	26
	5	0.21	-0.49	1.44		34	23
	6	0.21	-0.535 1.465		33	22
	7	0.2	-0.58	1.49		20	2
	8	0.2	-0.625	1.53		12	0
	𝜒 ()B
	Asteroid/Region	Effective surface	Strength to gravity	Reference
		Name			gravity		transition diameter
					(cm s -2 )		DSG (km)
		Ceres			28		1.75	Hiesinger et al. (2016)
	(Kerwan Basin)
		Vesta			25		1.94	Russell et al. (2012)
	(Rheasilvia Basin)
		Lutetia			4.7		49	Patzold et al. (2011)
	(Achaia Region)
		Mathilde			0.96		51	Thomas et al. (1999)
		Ida			0.7		69	Thomas et al. (1996)
		Gaspra			0.5		97	Thomas et al. (1994)
		Eros			0.4		120	Thomas et al. (2002)

6

is a 𝜒 6 test where the fit between the model and observed families (for parent bodies Dast > 100 km) is better than 2σ (i.e., probability >5%).

Table 2 .

 2 Compilation of asteroid surface gravities and strength to gravity transition diameter value used for the crater scaling laws in this paper.

	Ast. #	Asteroid/Region	Proper a	Proper e	Proper i	Ncross	Pi	Vimp
		Name	(AU)		(deg)		(10 -18 km -2 yr -1 )	(km s -1 )
	1	Ceres	2.7670963	0.1161977 9.6474113	642	3.455	4.860
		(Kerwan Basin)				(out of 681)
	4	Vesta	2.3615127	0.0987580 6.3923416	372	2.878	4.710
		(Rheasilvia Basin)				(out of 681)
	21	Lutetia	2.4352603	0.1292457 2.1461887	491	3.763	4.379
		(Achaia Region)				(out of 681)
	243	Ida	2.8616140	0.0456271 2.0883834	582	4.037	3.720
						(out of 682)
	253	Mathilde	2.6477821	0.2189155 6.5350556	666	3.723	5.237
						(out of 681)
	951	Gaspra	2.2097211	0.1475680 5.0786877	327	2.635	4.924
						(out of 682)
	2867	Steins	2.3635361	0.1082622 9.3526096	392	2.785	5.154
						(out of 682)

Table 3 .

 3 The intrinsic collision probabilities (Pi) and impact velocities (Vimp) for main belt asteroids observed by spacecraft. The first column is the asteroid number. The second column is the name of the asteroid, with the name of the region examined in parentheses where applicable.

	Asteroid	M	Empirical	𝜒 !" #	SFDs That	HH	𝜒 !$ #	SFDs That	Ivanov	𝜒 %& #	SFDs That
	Name	(# Crater	Fit:	Best Fit	Fit	Scaling:	Best Fit	Fit	Scaling:	Best	Fit
		Data	Best Fit		Within 1 s	Best Fit		Within 1 s	Best Fit	Fit	Within 1 s
		Points)	SFD			SFD			SFD
	Ceres-KB	22	8	2.84	6-7	8	4.61	7	8	9.44	-
	Vesta-RB	48	5	22.32	4,6-8	7	19.34	6	8	70.51	-
	Lutetia-AR	17	7	1.97	3-6, 8	7	2.24	5-6, 8	8	4.00	7
	Mathilde	12	2	2.67	1,3-6	4	2.60	1-3,5	8	4.98	3-7
	Ida	16	8	1.17	1-7	6	1.39	1-5, 7-8	8	1.93	1-7
	Gaspra	14	1	1.94	2-8	1	1.93	2-8	1	1.85	2-8
	Eros	12	2	3.72	1,3-8	1	3.85	2-7	8	2.27	1-8

The proper semimajor axis a, eccentricity e, and inclination i values were taken from the Asteroids Dynamic Site AstDyS (https://newton.spacedys.com/astdys/). The comparison population of 682 asteroids with Dast ≥ 50 km was taken from

[START_REF] Farinella | Collision rates and impact velocities in the main asteroid belt[END_REF]

, and Ncross describes the number of these bodies on crossing orbits with the target asteroid.

Table 4 .

 4 Compilation of results where model and observed crater SFDs were compared to one another. The index numbers 1-8 corresponds to the eight 𝑄 ! * asteroid disruption functions and eight model main belt SFDs (Figs. 1-2; Table

Acknowledgements

We thank Kevin Housen for his comments that helped us to better understand the nature of how asteroids might disrupt as they are spun up to high speeds, as well as the constructive comments from our referees. We also thank Peter Thomas and Mark Robinson for the use of their Eros crater database. The data generated in this work are reported in the text and figures. The work in this paper was supported by NASA's OSIRIS-REx mission through NASA Contract NNM10AA11C issued through the New Frontiers Program and grant no. 80NSSC18K0226 as part of the OSIRIS-REx Participating Scientist Program. The work of David Vokrouhlický was partially funded by a grant 18-06083S of the Czech Science Foundation.

Crater retention age of largest craters ~ 1 Ga, provided f ~ 10.

6.3

Table 6. Summary of results for asteroids investigated in this paper. For each asteroid, we list the taxonomic type, its dimensions, its association with a given asteroid family, comments on its probable crater retention age and the age match with auxiliary constraints, and the section of the paper where the asteroid is discussed.