Numerical simulations suggest asteroids (101955) Bennu and (162173) Ryugu are likely second or later generation rubble piles

Domain: earth_science
Published: 2024-07-05
Editorial outcome: accept

Abstract
Rubble pile asteroids are widely understood to be composed of reaccumulated debris following a catastrophic collision between asteroids in the main asteroid belt, where each disruption can make a family of new asteroids. Near-Earth asteroids Ryugu and Bennu have been linked to collisional families in the main asteroid belt, but surface age analyses of each asteroid suggest these bodies are substantially younger than their putative families. Here we show, through a coupled collisional and dynamical evolution of members of these families, that neither asteroid was likely to have been created at the same time as the original family breakups, but rather are likely remnants of later disruptions of original family members, making them second, or later, generation remnants. Our model finds about 80% and 60% of asteroids currently being delivered to near-Earth orbits from the respective families of New Polana and Eulalia are second or later generation. These asteroids delivered today in the 0.5-1 km size range have median ages since their last disruption that are substantially younger than the family age, reconciling their measured crater retention ages with membership in these families.

Introduction
Near-Earth asteroids (NEAs) are dynamically transient bodies that primarily originate in the main asteroid belt and survive for only an average of 10 Myr in near-Earth orbits before either impacting the Sun or a planet or being ejected from the solar system1. They transition from the Main Belt to near-Earth orbits through a handful of dynamical pathways that are typically orbital resonances with giant planets or Mars2. While any given NEA orbit is chaotic and impossible to directly trace back to a specific orbit or location in the Main Belt, it is possible, in some cases, to determine which dynamical pathway an NEA most likely followed2,3.

NEAs larger than a few hundred meters and smaller than 5–10 km are expected to be rubble piles, which are the accumulated debris from the collisional fragmentation of larger asteroids4–6. Models of the collisional history of the asteroid belt have shown that bodies smaller than diameter D ~ 10–20 km have collisional lifetimes shorter than the age of the solar system7. They are unlikely to have survived in their current state and are vastly more likely to be remnants of a more recent collisional event.

Bennu and Ryugu, rubble pile NEAs that have both been visited by spacecraft8,9, are highly likely to reach their specific orbits in near-Earth space via the ν6 secular resonance at the inner edge of the Asteroid Belt or via the Intermediate Mars Crossing population (IMC). The IMC population is filled largely by minor resonances, such as the overlapping resonances with Jupiter and Mars at 2.25 au (hereafter J7:2M5:9). Models of NEA delivery find that the probability of a body on Bennu’s current orbit has 82% and 18% probability to originate from the ν6 or IMC, respectively10. Ryugu has similar origins with 80% and 20% likelihood from the same pathways11. Updated NEO models, using a larger number of dynamical sources find similarly high probabilities of delivery from ν6 resonance where Bennu is 86% likely and Ryugu 88% (refs. 12,13).

Evidence of on-going asteroidal collisional evolution can be found by way of asteroid families, which litter the Main Belt (see refs. 14,15 for reviews). Here, the accumulated remnants of disruption events are linked by their similar orbital properties and sometimes through their physical properties such as spectral colors and albedo. Therefore, nearly all small (D < 10 km) asteroids in the Main Belt should be remnants of these big collisions, share a genetic link with a large parent body, and have a potentially huge number of relatives (family members) formed at the same time in the same event. Once a family is formed all the new members are subject to the same collisional environment as the parent body. The clock on the collisional lifetime of each new member starts ticking and eventually many will experience big impacts, resulting in the formation of new sub-families within the larger original family. It is common to find numerous smaller sub-families within larger, older, families (e.g refs. 16,17), and there are certainly large numbers of families and sub-families that are yet to be detected.

Similarly, all the new family members made in a collision are subject to size-dependent orbital drift due to the Yarkovsky effect18,19. The rate of semi-major axis change depends on an asteroid’s obliquity, thermal properties, and distance from the Sun, but importantly has a 1/D dependence, where D is the diameter of the asteroid, that drives the smallest members of a family to drift farther and faster than larger members. Constraints on drift rates allow for the spread of a family to act as a chronometer and allow estimates of families’ ages, or time of formation (see refs. 19–24) This semi-major axis drift is also the mechanism by which Main Belt asteroids transition to near-Earth orbits, as their drift often passes through orbital resonances with the giant planets or Mars. This can cause their orbits to be dynamically excited, i.e. their orbital eccentricity is increased, and some of them may subsequently encounter a terrestrial planet and be scattered into near-Earth space.

The knowledge of the dynamical pathway, combined with other spectroscopic and radiometric constraints, led to links between Bennu and Ryugu with the New Polana or Eulalia families10,11,25,26. These low-albedo families are found in the outer edge of the inner Main Belt (2.1–2.5 au), overlap with each other in their distribution of orbital elements and appear nearly indistinguishable in visible and near-infrared spectroscopic studies27–30. It is estimated that New Polana formed 1400 ± 150 My ago and Eulalia formed 830 +370 -100 Myr ago26, and either is currently capable of delivering 0.5–1 km objects like Bennu and Ryugu to the near-Earth population through the ν6 or IMC pathways.

An asteroid that formed after the initial collision that led to the formation of these families must endure throughout the entire lifespan of the family without being shattered by subsequent collisions or dynamically expelled before reaching its current near-Earth orbit. In the collisional environment of the Main Belt, asteroids the size of Bennu (D ~ 500 m) or Ryugu (D ~ 1 km) are estimated to have lifetimes against collisional disruption of only ~200 Myr and 440 Myr, respectively7,31. To be 1st-generation members of either family they must survive several collisional lifetimes. The alternative is that a larger 1st-generation remnant that was more resilient to collisional disruption due to its larger size, but that drifted more slowly for many millions of years, was itself eventually disrupted and produced Bennu or Ryugu as part of a second, or later, generation sub-family.

This difference has implications for interpreting the observed geology on both bodies and interpreting the nature of the returned samples. The crater retention age for Bennu and Ryugu are 10-65 Myr (ref. 32) and <30 Myr, respectively8,33. If it is probable that they are second generation or older remnants of a family then their young surfaces may be reconciled with their proposed membership in these families. If not, then membership in these families is challenged or geologic processes such as global resurfacing event(s) would need to be responsible for resetting their crater retention ages. Meanwhile, indications of a few or numerous large collisional events in the returned samples could illuminate a likely history or reveal gaps in modeling assumptions.

In this work we present results from a combined dynamical and collisional model of asteroid family evolution. It shows that it is more likely that asteroids currently being delivered from the Eulalia and New Polana families are second, or later, generation remnants and not 1st generation remnants.

Results
Modeling family generation and evolution
Previous efforts have carried out extensive N-body simulations of members of the New Polana and Eulalia families drifting across the main asteroid belt26. Those N-body simulations included the planets and thus incorporated the effects of the orbital resonances directly, but did not account for collisional removal or the generation of new simulated bodies in collisions. Other models have developed sophisticated routines to capture the effects of collisional evolution, and its large-scale effects on the size frequency distribution (SFD) of the asteroids in the main asteroid belt7,34. These large models cover such large timescales and large swaths of the Main Belt that they don’t take into account individual asteroid families and do not incorporate the ability to track individual asteroids and their changing orbits. The modeling approach here, like the former, tracks individual asteroids’ orbital drift and their removal by resonances, and like the latter, also considers their collisional disruption and the re-distribution of mass to smaller bodies.

Asteroid families are modeled here through Yarkovsky effect-driven orbit drift and a Monte-Carlo approach to collisional evolution. In this model, the initial asteroid family, and those created through later disruption, have half of the disrupted target’s mass in its largest remnant7. The next largest body, the largest fragment is modeled to always be a constant mass ratio of the largest remnant7. It is the largest fragment that is the anchor for the power law size-distribution of asteroids for the rest of the family (see Fig. 1). Previous models of the Main Belt SFD evolution used the largest size ratio of 0.8 (where mass of the largest fragment is 1/2 that of the largest remnant) and re-created many aspects of the Main Belt size distribution7, while others left it as a free parameter that varied as a function of the impact properties34. Here, the largest fragment relative to the largest remnant is left as a free parameter, varying between 0.2, 0.4, 0.6, and 0.8 the size of the largest remnant (0.008, 0.064, 0.022 and 0.5, respectively, the mass). This range includes the previously used values and brackets the observed values for a handful of well-studied inner Main Belt C-type families and is in line with numerical simulations of disruption and reaccumulation4,35,36 (Fig. 1). The models here assume that the slope of the cumulative SFD is constant, similar to ref. 34, and unlike ref. 7 who used a piecewise SFD. Values for that slope found in simulations36 and in observed asteroid families span a wider range14. Again, the power law exponent values for family SFD slope were treated as a free variable. Slope values of −3, −4, and −5 were tested to bracket the most likely typical initial family SFDs.

To initiate a simulation the diameter of the original parent body is selected and the fragments sizes are then determined using the formula escape velocity of the target bodies and this is included for each initial family and subsequent sub-family following the formulations of ref. 20. The objects are given random obliquities on a unit sphere that determine their Yarkovsky drift rate and direction when combined with their sizes. The evolution of the objects due to thermal torques is complex and poorly constrained at small sizes where spin rate and obliquity changes are expected due to the YORP effect18,19,26. Changes in obliquity, in turn, affect orbital drift rates. Thus, the asteroids were modeled in two ways. First, the initial obliquities did not change throughout the simulation and drift rates for each simulated asteroid were calculated from Bennu’s measured rate of −284 m/year that is scaled to its diameter, semi-major axis, and obliquity39,40. Second, they were modeled with the YORP cycles and stochastic YORP functions and parameters described in ref. 26. These formulations require tracking spin rates, which have direct interplay with changing obliquity and reorientation due to impacts. The net effect is for obliquities to evolve towards 0 or 180°, with some reorientations due to impacts on very slowly spinning bodies26. Unless otherwise noted results are shown for the first modeling method with unchanging obliquities.

An asteroid in the main asteroid belt will experience collisions with a size distribution of impactors, and eventually could be the impactor itself when hitting a larger body. For this effort we do not attempt to accumulate every impact, rather we follow previous works that modeled main belt collisional evolution26 and determine when a disruptive impact is likely to happen based on the collisional probabilities for the main asteroid belt. This is calculated via the size-dependent specific impact energy equal to Q*D, which is the ratio of the projectile’s kinetic energy to the target’s mass leading to half of the target’s mass being in the largest remnant. The timescale for a disruption to occur is treated with a disruption law based on measured physical properties of boulders on Bennu31 and is similar in these size ranges to that used in previous models of Main Belt collisional evolution7,41.

When a simulated asteroid is disrupted it is replaced with a size distribution of smaller remnants. These small remnants follow the same rules for the largest remnant, largest fragment and SFD slope as for largest remnant, largest fragment and the power law slope for the remaining distribution of family members. Nominal settings utilize a largest remnant to be 0.8 times the size of the target (1/2 the mass) and the largest fragment to be 0.4 times the size of the largest remnant and a power-law SFD slope of −4.0. The size of the parent body is not precisely known for either family but is constrained by the size of the observed largest family members, (495) Eulalia and (142) Polana (with diameters of about 40 km and 55 km respectively), and by the nature of their family SFDs showing that neither family appeared to be super-catastrophic in nature25. Integrating the total mass of known and possible family members and accounting for various loss mechanisms, the estimated parent body for the Eulalia family had a diameter Dpb = 100-160 km (ref. 25). The New Polana family is more challenging to estimate25. Given the asteroids similar size and their families’ similar SFDs, here we consider simulating Dpb = 80 km, 100 km and 160 km for both families to attempt to bracket all possibilities.

The collisional cascade in the Main Belt creates debris down to very small sizes, but Bennu has a diameter about 0.5 km and Ryugu about 1 km (refs. 37,38). Hence for the focus of this work it is only necessary to track the delivery of simulated asteroids down to those sizes, and not below. Collisional timescales for simulated asteroids are a function of the Main Belt background asteroid population as previously calculated2, rather than the number, size, or orbits of the asteroids simulated here. Thus, ignoring very small objects does not change the collision calculation results.

The initial semi-major axes of the objects in a family are set to the same value as the parent bodies, where Eulalia has semi-major axis a = 2.487 au and Polana a = 2.41 au. There is some initial size-dependent dispersal of orbital semi-major axis that scales with the for the initial family and explore the same range of values for each variable. The total integrated mass is limited so that no mass beyond that of the disrupted body is added to the simulation. Similarly, no bodies smaller than 0.5 km are added. Newly formed bodies are started near the location of their recently disrupted parent, whose size sets the scaling for the the size-dependent spread in the semi-major axis relative to the parent’s. They are given random obliquities to determine their drift rate and direction.

Simulated asteroids can be removed via collisional evolution or by orbital dynamics—specifically if their orbits drift across resonances with planets. As a simulated asteroid drifts across the main belt, it can encounter mean motion or secular resonances with planets, and face certain, or probable removal from the asteroid belt by eccentricity increase. Drifting inward from the inner main belt, the ν6 secular resonance will be 100% efficient at removing a body if it reaches a of ~2.15 au. The location of ν6 is inclination-dependent, and is, in places, close to the location where a body could start to cross the orbit of Mars, which is eccentricity-dependent. The Eulalia and New Polana families have relatively low eccentricity and inclination, allowing for a simple calculation of semi-major axis overlap with the location of ν6. The outer edge of the inner Main Belt is the 3:1 resonance with Jupiter, which, like ν6, is nearly 100% efficient at removing bodies from the belt. In the middle of the inner belt lies the overlapping Jupiter 7:2 Mars 5:9 mean motion resonances at a = 2.255 au (hereafter J7:2M5:9), which could be encountered from either direction (inward or outward drifting—depending on the simulated asteroid’s starting location) and is about 33% efficient at removal of km-sized asteroids26,42. Note that the removal efficiency is drift rate-dependent, which itself is size dependent26,43 and these removal calculations are valid only for the km-scale objects.

Furthermore, while the ν6 and 3:1 are effectively 100% efficient at removing asteroids from the asteroid belt, neither are as efficient at placing them on near-Earth orbits, as some are ejected from the Solar System or achieve some other fate. Here, since previous modeling has found Bennu and Ryugu’s likely dynamical pathway to their current orbits, there is no need to blend the relative probability of NEA delivery between these two removal sources. Instead, we are exploring the collisional history of objects that follow this same dynamical path. Since this effort will rely on a Monte Carlo approach the total number of initial simulated asteroids in a family can be as large as computationally feasible as long as it adheres to the family size distribution. The numbers and sizes of fragments made in each subsequent collision and disruption must follow the same rules and make proportionally similar numbers, while still adhering to the family size distribution. Thus, nominal simulations below begin with 300,000 simulated asteroids, although the number grows throughout each simulation depending on the number and nature of the collisions in each.

1st- or Nth-generation remnants
The delivery of different generations of asteroidal fragments varies over time due to the different family locations relative to the important escape routes. For example, families farther from the J7:2M5:9 and ν6 take longer to deliver their first fragments, allowing more time for collisional evolution. The evolution for the New Polana family illustrates the time-varying delivery from the different pathways, where the family location at 2.41 au requires asteroids to cross the J7:2M5:9 prior to reaching the ν6 resonance. This results in an initial spike at about 180 Myr that precedes a spike from ν6 coming nearly 100 Myr later (Fig. 2a). While the J7:2M5:9 is only 1/3 efficient at removing asteroids, the extra roughly 100 Myr that it takes them to reach the ν6 allows for many to be removed via further collisional evolution. As the time since family formation increases (family formation is at 0 Myr in the figures), the effects of collisional evolution increase the number of 2nd- and later-generation of fragments reaching the resonances (Fig. 2b, c).

A simple property extracted from each set of simulation parameters is the crossover time when the number of Nth generation asteroids begin to dominate in number over 1st generation bodies. For the set of nominal parameters where the largest fragment is 0.4 the size (and 6% the mass) of the largest remnant and the family cumulative SFD slope is −4.0, this occurs at 450 Myr for the New Polana family with an assumed parent body size of 100 km (Fig. 2b, c). For the same parameters, the Eulalia family reaches the crossover time at 510 Myr, due to its center location at 2.448 au, which is farther from either of the important resonances (Fig. 3b, c).

Delivery from the Eulalia and New Polana families have similar dependencies on the two parameters: where crossover times increase for 1) steeper family SFD slopes and 2) smaller largest fragment sizes (Fig. 4). The steeper family SFD slopes decrease the number of smaller bodies generated during each disruption event making the re-filling of the small end of the size distributions slower, which in turn decreases the number of fast-drifting bodies that can reach escape resonances. Similarly, the size of the largest fragment determines the mass of disrupted objects that are preferentially retained, but also can result in a large number of fragments being created below the resolution of the model at 0.5 km, where they are no longer tracked. The smaller largest fragments result in even more mass being put into asteroid with sizes below that limit and decrease the number of new, km-sized, fragments to be tracked, and subsequently causing later crossover times.

The New Polana family is estimated to be 1400 ± 150 Myr (ref. 26) and there are many millions of years that 1st-generation fragments dominate delivery. The relevant delivery time is 1400 Myr after family formation and corresponds to asteroids escaping the Main Belt today. For the nominal modeling scenario, the delivery at 1400 Myr after family formation is 21% from 1st generation, compared to 41% 2nd generation and 22% 3rd generation (Table 1). The younger Eulalia family had a delivery composition of 38% 1st generation, 39% 2nd generation and 18% 3rd generation after 830 Myr (Table 1).

Across all tested ranges in the largest fragment size, family SFD slope and parent body size there was variation on order of a few 100 Myr for the crossover times (Fig. 4). The New Polana family has crossover times that range from 340 Myr to 520 Myr while Eulalia spans crossover times from 410 Myr to 610 Myr for the same suite of parameters. The inclusion of YORP cycles and stochastic YORP shortens the range of crossover times to 330 Myr to 430 Myr for New Polana and 400 to 540 Myr for Eulalia. YORP cycles lead to obliquities of 0 or 180°, maximizing drift inward or outward which result in shorter crossover times. The range of crossover times for both families, for all tested parameters with or without YORP cycles and stochastic YORP indicates that Nth-generation objects dominate delivery to the inner solar system at timescales shorter than the expected family age.

Formation ages of delivered asteroids
In the models presented here, the formation age of each simulated asteroid is the time when it enters the simulation from either the family-forming event or the subsequent disruption of a larger fragment. For the D < 1 km objects that reached the inward resonances in the nominal New Polana simulation 33% and 21% had asteroid formation ages younger than 60 Myr and 10 Myr, respectively, with a median age of 190 Myr (Fig. 5). For the Eulalia simulation 19% and 15% were younger than 60 Myr and 10 Myr respectively with a median asteroid age of 330 Myr. The differences primarily reflect the different family ages, where the New Polana family age of 1400 Myr is nearly three collisional timescales for a 1 km object which greatly reduces the number of 1st generation, and very old objects, that could survive. The greater fraction of Nth generation fragments from New Polana (80%) decreases the median asteroid formation age relative to Eulalia, which has a larger (38%) population of very old 1st-generation bodies that are delivered.

The crater retention ages for Bennu and Ryugu are 10-65 Myr and <30 Myr respectively8,32,33. From the suite of model parameters considered here, the case with highest fraction of ages below 60 Myr and the lowest median age assumed a family SFD with a slope of −3 and a largest fragment size of 0.8 times the largest remnant for the New Polana family. Here the median age was 90 Myr, with 46% and 34% with 60 Myr and 10 Myr, respectively.

Simulated asteroids from the 2nd- or later generations were, by definition, not directly created in the original family breakup and thus have a direct parent smaller than the family parent body. Most Nth-generation objects (97%) are the largest remnant of a disruption of only a slightly larger body—they are the largest remnants of their direct parent. For the set of nominal parameters, the smallest parent had a diameter of 0.63 km, where the relationship between smallest allowable direct parent depended on the simulation settings (Fig. 6). Increasing the relative size of the largest fragment decreases this fraction as more larger bodies are produced at each disruption and intermediate sizes were subsequently more populated in the distribution. Steeper SFD slopes increase the fraction of objects that are the largest remnants.

Methods
Collisional disruption calculations
The code steps forward in 10 Myr time increments. Within each timestep each simulated asteroid is tested, based on its size, to determine if it suffers a disruptive collision. This was calculated with the Q*D formulation and coefficients from ref. 31, shown in to be similar to the results in ref. 7 in the size range of interest. Notably ref. 31 disruption law was determined from boulder strengths on Bennu and Ryugu. Changing the simulation timestep to 5 Myr did not change the crossover point for Eulalia with these settings and New Polana changed the crossover time by <2%.

Yarkovsky drift rates for simulated asteroids
The semi-major axis drift rates are anchored on the measured drift rate of 1.9e−3 au/Myr measured for Bennu39,40, and here they are scaled by their semi-major axis, their size ratio and their obliquity (where Bennu’s obliquity is nearly 180°).

da/dt = (da/dt)bennu * sqrt(a/a_bennu) * (D_bennu/D) * cos(obliquity) (1)

where a and D are the target bodies semi-major axis and Diameter, a_bennu is Bennu’s semi-major axis, D_bennu is Bennu’s diameter and obliquity is the target bodies orbital obliquity. The measured rate for Bennu is very accurate owing to the long-duration radio tracking provided by the OSIRIS-REx mission40, but it is possible that Bennu is not representative of the modeled families or that there is significant variation among members of any given family. Changing the drift rate by a factor of 3 in either direction (faster or slower) does change the outcomes. The faster drifting case had the same crossover time of 470 Myr, but with a strong shift to second generation remnants and a decreased median formation age for delivered objects. Notably there are very few remnants being delivered at 1400 Myr relative to the earlier peak delivery times as the total numbers in the family are heavily depleted due to the fast drift to resonances. The slow drifting case had a moderate shift to 3rd and later generations but a longer crossover time of 720 Myr, after which Nth-generation fragments dominate.

Flora and Didymos from Flora and Baptistina families
Itokawa, an LL chondrite72, is likely to have originated from the Flora family at 2.202 au at the innermost edge of the main belt41, and Didymos is a good match for the Baptistina family in a similar region of the inner Main Belt63. Using those families as possible origins, and the same drift rate as Bennu and the rest of the simulation nominal parameters the outcomes strongly support Itokawa as Nth-generation (95%) and Didymos 1st generation (77%).

YORP cycles and stochastic YORP
This implementation of the model requires tracking spin rates for each simulated asteroid as change in obliquity is spin-rate dependent. Both spin rate and obliquity change on timescales as short as ~1 Myr, and thus the simulation timestep was decreased to 0.01 Myr. The functions for obliquity change and spin rate change are those described in refs. 20,26, and collisional reorientation, described in ref. 20 and derived from ref. 73. These formulations require several additional parameters, briefly described here, but more detailed explanation in ref. 26, that provides numerous fits to the evolution of the Eulalia and New Polana families.

The strength of the YORP effect is controlled by c_YORP, and was set at 0.6 throughout. δ_YORP controls the bias for spin-up and spin-down following a reorientation, and is set at 0.5. τ_YORP is the timescale for the stochastic YORP process of selecting a new YORP condition for spin rate change but not altering the obliquity. In Bottke et al. 2015 this was done in a size-dependent manner and that is re-created here where for diameter larger than 3 km τ_YORP = 1 Myr, for diameters larger than 1.9 km and smaller than 3 km τ_YORP = 0.5 Myr and for diameters smaller than 0.9 km τ_YORP = 0.25 Myr.

Peer review rounds
[{"round": 1, "reviewer_comments": {"reviewer_1": "Authors: Walsh et al.,\nTitle: Bennu and Ryugu likely 2nd or later generation rubble piles\nRubble pile asteroids are a subject of interest. They represent an aggregation of debris formed in\ncatastrophic disruptions between \"parent\" asteroids. Each collision of this kind has the potential to\ngive rise to an asteroid family. This study employs models of collisional and dynamical evolution of at\nleast two distinct asteroid families (New Polana and Eulalia). The aim of the study is to determine the\nprobability of asteroids Ryugu & Bennu being direct remnants (\"1st generation\") bodies, from the\ninitial disruption event that created their respective families. Alternatively, it seeks to ascertain if\nthese asteroids have undergone extensive evolution and experienced multiple disruptive events over\ntheir lifetime.\nThe Near-Earth asteroids Ryugu and Bennu hold particular interests to the wider audience thanks to\nthe samples retrieved from their surfaces by JAXA's Hayabusa2 (from Ryugu) and NASA's OSIRIS-REx\n(from Bennu). This research is timely as the sample collected from Bennu is scheduled to reach Earth\non September 24th.\nThe manuscript is well written and layed out in a logical manner. However, I do raise some concerns\nregarding the assumptions made for the modeling work and the implications of the presented work.\nBelow is a list of “issues” requiring substantial attention from the authors.\nIn the calculations for collisional disruption, the authors have relied on the Q*D laws as derived by\nBallouz et al. (2020). It's worth noting that these scaling relationships have been derived from fitting\nto measurements of cratered and cracked boulders on Bennu's surface. While this approach may be\nsuitable for the initial parent body, which can be treated as a monolith, it may not hold true for\nsubsequent impacts on modified asteroids that may themselves be rubble piles. This raises the\nquestion of how the results would be affected if one were to assume a different Q*D curve for the\nsubsequent collisions. Additionally, I'm curious about the robustness of the presented results in light\nof the choice in Q*D curves.\nThe implications of this study remain somewhat ambiguous, lacking a definitive argument for its\nnovelty. For example, on line #110, the authors suggest implications for interpreting the returned\nsample; however, this point is not elaborated beyond the introduction. Given that samples from\nRyugu have already undergone analysis and resulted in numerous published papers, it begs the\nquestion of whether any parallels can be drawn.\nAt line #485, the authors have extended the same analysis used for C-type asteroids Ryugu and\nBennu to S-type asteroids Itokawa and Didymos. Nevertheless, there is presently no compelling\nevidence to suggest that the Q*D equations would be applicable to the less porous S-type asteroids.\nTherefore, additional explanation is needed or the paragraph removed altogether.\n#102: Some references are needed\n#186: There should be a comma after (Bottke et al. 2022), not full stop.\n# There are a couple of repetitions throughout the manuscript. See Lines #111 & #351", "reviewer_2": "The manuscript \"Bennu and Ryugu likely 2nd or later generation rubble piles\" investigates\na scenario in which both asteroids are 2nd or later generation fragments from an inner\nasteroid belt body. The study is motivated by apparently very different ages of asteroid\nfamilies where the asteroids are supposed to originate from and the crater retention ages.\nTo this purpose, the authors employed a novel, Monte-Carlo-based approach that simultaneously\ntreats the collisional and dynamical evolution. This is one of the first works in this\nrespect and is appropriate for the problem studied.\nIndeed, the model is based on different assumptions and approximations, but the authors\nalso investigated how sensitive the results are concerning some of the input parameters.\nFor instance, the authors show that the results do not change significantly for different\nslopes of the size-frequency distribution nor when varying the sizes of the parent bodies.\nThe manuscript is well-written, and the figures and tables are appropriate and informative.\nTherefore, it fully deserves publication.\nThere are only minor issues to be addressed:\nGenerally, if I'm not missing something, it is not described when exactly a test particle\nis assumed to be entering the NEAs region, i.e., when it is considered to be 'delivered'.\nThe nu6 and 3/1 resonances have been supposed to remove 100% of objects, right? But what\nhappens to them? I assume that a fraction from each resonance is considered to be delivered\ninto the NEA region, but I cannot find precise information. This aspect of the model should\nbe better described, and missing details should be provided.\nline 41: \"All but the largest NEAs are expected to be rubble piles\"\nThis sounds like a too bold statement. Are you claiming that also all objects below ~100 m\nare rubble piles? Even super-fast rotators? Please rephrase and/or be more specific.\nline 54: \"Models of NEA delivery find that the probability of a body\non Bennu's current orbit has a 82% and a 18% probability to originate from the ν6 or IMC,\nrespectively\"\nIt would be helpful to provide here also the probabilities and routes from the model by\nGranvink et al. (Icarus, 2018)\nline 61: Nesvorny et al. 2015 is definitely an appropriate reference for asteroid families.\nYet, it might be good to cite along a more recent review by Novakovic et al. CeMDA, 2022,\nwhich is some aspect complementary.\nline 157, Figure 1: The reference for the family ages should be added, Bottke et al. (Icarus, 2015)\nline190: \"There is some initial dispersal of orbital semi-major axis that scales with the escape velocity\nof the target bodies (see Vokrouhlicky et al. 2006), but here the very long timescale of their evolution\noverwhelms the initial scatter of initial orbital semi-major axis, so each family member starts with\nidentical semi-major axes.\"\nI agree that the initial dispersions should play a minor role compared to post-impact evolution. It\nseems,\nhowever, that it would be simple to implement this aspect into the model. At least for larger parents,\nit might be relevant because the expected ejection velocities are higher. Also, in the case of Eulalia,\nwhich is close to the 3/1 resonance, adding the initial dispersion might place some bodies on the\nother\nside of the resonance.\nline 244: \"Note that the removal efficiency is drift rate-dependent, which itself is size dependent\n(Bottke et al. 2015) and these removal calculations are valid only for the km-scale objects.\"\nMilic Zitnik & Novakovic (ApJ, 2016) studied the effect of some MMRs on the semi-major axis drift\ndue\nto the Yarkovsky effect. Though it may not be directly relevant to this work, it could be worth\nmentioning it here as well, as, in general, this analysis could be very useful in the model\nlike the one presented in the manuscript.\nline 379, Figure 6: This figure's caption says there were 300,000 initial asteroids. Was that the\nnumber\nof test particles in all simulations? Please specify. Also, this number should be stated earlier when\nthe simulation set-up is described.\nline 485: \"Flora has a similar age to New Polana...\"\nPlease add a reference for the age of the Flora family, Vokrouhlicky et al. (AJ, 2017)\nline 513: In the Yarkovsky scaling relation, it seems that the scaling with the obliquity of Bennu is\nmissing.\nOr is it omitted because it is close to 180 deg, implying that cosine is close to -1, and the sign is\nneglected?\nPlease check and fix or add missing details.", "reviewer_3": "Main belt asteroids can be disrupted by catastrophic collisions, leading to the formation of\nfragments. These fragments are referred to as a \"family\" of the original asteroid. The orbits of these\nfragments, now smaller asteroids, can be altered by the Yarkovsky effect. Additionally, these\nfragments can be further disrupted by subsequent collisions with other small bodies. When a main\nbelt asteroid fragment reaches a certain radial location, known as an orbital resonance, it can be\ndynamically transferred to the near-Earth region.\nThis study utilized numerical simulations to examine the orbit evolution of small bodies driven by the\nYarkovsky effect. It combined this with the catastrophic collisional evolution using a Monte-Carlo\napproach. The authors aimed to investigate specific asteroid families, such as New Polana and\nEulalia, and to discuss their collisional and orbital evolution. Recently, spacecrafts have visited near-\nEarth asteroids Ryugu and Bennu. These asteroids are thought to be members of the collisional\nfamilies of New Polana or Eulalia. However, the surface ages of these asteroids, as reported, are\nsignificantly younger than the estimated age of their parent family's formation. Through numerical\nsimulations, the authors determined that Ryugu and Bennu are not fragments from the initial family-\nforming collision. Instead, they are products of subsequent, additional collisions.\nIn my assessment, the findings presented are both intriguing and relevant to the current scientific\ncontext. The manuscript is largely well-written and the content is effectively conveyed. As such, I\nbelieve it meets the standards for publication in Nature Communications. However, before I can\nrecommend its acceptance, I would like the authors to provide further clarity on their numerical\nmodels. Additionally, I urge them to discuss the consistency of their findings with the returned\nsamples from Ryugu — and from Bennu, if available. Addressing these points is crucial for validating\nand bolstering the credibility of their numerical results.\nMajor comments.\n(1) In lines 125-130: I would appreciate it if the authors could first verify whether their new approach\n— encompassing both orbital and collisional evolution — are able to accurately reproduce the Size-\nFrequency Distribution (SFD) of asteroids within the main asteroid belt. Additionally, could the\nauthors discuss any updates or differences in their findings, especially in comparison to the results\npresented in Bottke et al. 2005 and Morbidelli et al. 2009?\n(2) In line 415, the authors note that collisional evolution might significantly modify the parent body\nmaterials, leading to a mixture of materials from both the impactor and the target. However, recent\nliterature, specifically Nakamura T. et al. 2022 in Science, has indicated that sample material from\nRyugu doesn't show evidence of impact shock. Furthermore, Tatsumi et al. 2021 reported that the\nfraction of exogenic materials in the Ryugu sample is comparatively small. One would anticipate that\nmultiple catastrophic impacts would both shock the materials and promote mixing between the\nimpactor and target. This presumption seems to contradict the sample analysis results from Ryugu,\nas mentioned earlier. Could the authors provide an explanation for this apparent discrepancy?\nIn line with the aforementioned concerns, I request the authors to provide specific details regarding\nthe impact conditions encountered during the collisional evolutions. Can the authors elucidate on\nthe impact velocity and the impactor-to-target size ratio as derived from their simulations? Following\nthis, I would appreciate it if the authors could estimate the consequent increase in pressure and\ntemperature resulting from the specified impact conditions. I hypothesize that if the impactor-to-\ntarget size ratio is close to unity, it would likely lead to enhanced material mixing. Could the authors\ncomment on this?\n(3) In lines 185-190: The authors have posited that collisions with background asteroids primarily\ninfluence family members. However, immediately following the initial catastrophic collision that\nresults in the formation of an asteroid family, the fragments would likely have orbits in close\nproximity to one another, given that they would share an almost identical point of impact.\nConsequently, right after the impact, there's a plausible scenario where these fragments (or family\nmembers) might collide with each other frequently, leading to further fragmentation. Could the\nauthors provide a discussion or estimation on the collision timescale under such conditions?\nMinor comments:\nHirabayashi et al. 2020; Hyodo & Sugiura 2022 also discuss the shape change and resurfacing of small\nasteroids."}, "author_response": "To whom it may concern,\nBelow we endeavor to answer all of the comments of the two reviewers and describe or list the\nchanges that we made to the manuscript. We believe that we adequately answer all of the\nconcerns and that the main findings of this work are unchanged. You will find the Reviewers\ncomments in Italics and our responses in normal font.\nReviewer 1:\nRubble pile asteroids are a subject of interest. They represent an aggregation of debris\nformed in catastrophic disruptions between \"parent\" asteroids. Each collision of this kind\nhas the potential to give rise to an asteroid family. This study employs models of\ncollisional and dynamical evolution of at least two distinct asteroid families (New Polana\nand Eulalia). The aim of the study is to determine the probability of asteroids Ryugu &\nBennu being direct remnants (\"1st generation\") bodies, from the initial disruption event\nthat created their respective families. Alternatively, it seeks to ascertain if these\nasteroids have undergone extensive evolution and experienced multiple disruptive\nevents over their lifetime.\nThe Near-Earth asteroids Ryugu and Bennu hold particular interests to the wider\naudience thanks to the samples retrieved from their surfaces by JAXA's Hayabusa2\n(from Ryugu) and NASA's OSIRIS-REx (from Bennu). This research is timely as the\nsample collected from Bennu is scheduled to reach Earth on September 24th.\nThe manuscript is well written and layed out in a logical manner. However, I do raise\nsome concerns regarding the assumptions made for the modeling work and the\nimplications of the presented work. Below is a list of “issues” requiring substantial\nattention from the authors.\nIn the calculations for collisional disruption, the authors have relied on the Q*D laws as\nderived by Ballouz et al. (2020). It's worth noting that these scaling relationships have\nbeen derived from fitting to measurements of cratered and cracked boulders on Bennu's\nsurface. While this approach may be suitable for the initial parent body, which can be\ntreated as a monolith, it may not hold true for subsequent impacts on modified asteroids\nthat may themselves be rubble piles. This raises the question of how the results would\nbe affected if one were to assume a different Q*D curve for the subsequent collisions.\nAdditionally, I'm curious about the robustness of the presented results in light of the\nchoice in Q*D curves.\nWe thank the reviewer for these very good comments. We address them in turn.\nFirst, the dependence of our results on Q* curves is a very interesting question. We\nD\nnote that the Ballouz et al. (2020) Q* scaling relationship transitions from strength- to\nD\ngravity-dominated at an asteroid diameter of 160 m (Otohime Saxum’s diameter). The\ngravity-dominated part of the Ballouz et al. (2020) is derived from observations of\ncraters on large asteroids and scaling constants derived from the hydrocode disruption\nstudies of Leinhardt & Stewart (2012). Leinhardt & Stewart (2012) compiled data from\ncatastrophic-disruption simulations from various sources and found that for (i) a variety\nof target compositions, (ii) five orders of magnitude in size (~0.4–4e5 km) and (iii) nine\norders of magnitude in impact energy, the Q* relationship could be described by a\nD\nnarrow range in scaling constants. Therefore, the Q* disruption relationship would be\nD\napplicable to gravity-dominated rubble piles regardless of their composition. Note that\nQ* relationships describe a specific impact energy. As it is detailed later in our\nD\nresponse to reviewers (page5 - paragraph starting with “Indeed, Leinhardt & Stewart\n(2012) compiled data….”) it is understood that C- and S-type might have very similar\nQ* relationships in the gravity regime (rubble piles), so a C-type would indeed require\nD\nless energy to catastrophically disrupt than a S-type of the same size, due to purported\ndifferences in their bulk densities.\nThe work presented here really only deals with collisions in the size range of\n0.5km-~10km, where the vast majority of the collisions are happening at the small end.\nIn this range the expected collisional lifetime between Bottke et al. (2005) and Ballouz et\nal. (2020) differs by no more than a factor of ~2x - with Bottke’s yielding a shorter\nlifetime. Also note that the recent scaling law from Bottke et al. (2020) is even closer to\nBallouz et al. (2020). We did a series of test runs where we simply use a 3x larger or\nsmaller lifetime than the nominal Ballouz2020 case. The results are mostly intuitive -\nthey definitely change the timescales for crossover times - and we put a short\ndescription of this in the Discussion section. Here, though we summarize - the simple\ntest we did was to re-run with the Q* law as formulated and simply use a multiplicative\nfactor of 3x to lengthen or shorten the collisional lifetime from the calculations. For the\nfiducial Polana case this lengthens the crossover time to 1430Myr or shortens it to\n260Myr. The nominal value was 450Myr, so this finds that the trend is clear, but not\nsimple and linear. The 3x lifetime case this increase in lifetime means that 1km\nasteroids can survive 1.2Gyr and 0.5km can survive for 600Myr.\nWe addressed this directly in the Discussion with a new paragraph:\n“The results here have a strong dependence on the collisional lifetime for rubble\npile asteroids. A simple experiment of lengthening or shortening the lifetime for\neach asteroid by a factor of 3x resulted in extending and decreasing the\ncrossover time in the fiducial New Polana simulation (1430Myr and 260Myr). While\nthis is not a linear dependence, it is correlated and important. Two recent\ndisruption laws built to satisfy different constraints from the Main Belt population\nestimate collisional lifetimes within a factor of ~2 in the size ranges modeled here\n(Bottke et al. 2005, Ballouz et al. 2020). Thus, within the framework where these\nworks are good estimates for collisional lifetimes, the results presented here are\nrobust. If new disruption laws emerge from future works these results very\nroughly scale as do the collisional lifetimes.”\nRef:\nLeinhardt, Z. M. & Stewart, S. T. Collisions between gravity-dominated bodies. I.\nOutcome regimes and scaling laws. Astrophys. J. 745, 79 (2012).\nThe implications of this study remain somewhat ambiguous, lacking a definitive\nargument for its novelty. For example, on line #110, the authors suggest implications for\ninterpreting the returned sample; however, this point is not elaborated beyond the\nintroduction. Given that samples from Ryugu have already undergone analysis and\nresulted in numerous published papers, it begs the question of whether any parallels\ncan be drawn.\nFair point - we should have included more direct links to the published results from\nRyugu. In fact, its really not as simple as asserted in the initial manuscript and so it has\nbeen valuable to add more.\nFundamentally the lack of shock features in Ryugu samples is very interesting. It raises\na general question on how Ryugu materials have avoided receiving substantial shock\nmetamorphism during its collisional evolution from a 100-km-sized parent body to the\ncurrent size. However, we don't think this limits the number of large-scale impacts their\nvarious parent bodies have witnessed, which is the main focus of this study. If a\nmeteorite or sample is shocked, it means the parent was shocked. However, the lack of\nshock on a single meteorite or sample does not necessarily mean that there was no\nshock and no impacts on the parent. It just means that this meteorite or sample wasn’t\naffected by it due to lack of proximity or other effects. The samples are just not perfect\ntracers of collisional history. We have tried to clarify this in the text, and we include\nsubstantially more text below for the reviewers.\nSo, there are a couple of factors that could account for the lack of shock features in\nRyugu's samples. First is that the mass fraction of impact-heating in the parent body\nis estimated to be very small. Nakamura et al. (2022) conducted systematic impact\ncalculations using physical properties of Ryugu materials measured using the actual\nreturned samples and showed that the impact-induced shockwave decays so quickly\nthat it should not leave strong shock features in most of the target volume. Fig. S47\ncaption in that paper’s Supplementary Materials indicates that the heated region up to\n>90 °C is limited within the distance of ~2 impactor diameters from the impact point.\nThis volume is ~3% of parent body. If we take 250 °C, which is the lowest stage of\nthermal metamorphism for carbonacenous chondrites (Nakamura et al. 2005) as a\nthreshold, mass fraction predicted from Fig. S47 (Supplementary Materials) becomes\neven smaller (i.e., <1%). Such rapid shockwave decay and small shock-heated volume\nis caused by the high porosity nature of Ryugu materials. If Ryugu materials are from\n99% or 97% of unshocked volume as T. Nakamura et al. (2022) suggests, then we will\nmight expect to find little shocked materials in the samples.\nSecond, highly shocked materials generally might be expected to receive higher\nejection velocities. This would help shocked materials escape from the impact site. This\nwould further reduce the mass fraction of shocked materials in re-accumulated\nfragments loosely bound with their weak self-gravity.\nFurthermore, regarding the small mass fraction of impactor materials on Ryugu, some\nof this team is writing a paper on the effect of oblique impacts (Sugiura et al. in prep).\nPreliminary results indicate that impacts with 45-deg or higher measured from the\nnormal would result in impact mass mixing ratio as low as Ryugu (~10^-6) Sugita et al.\n(2023) (LPSC abstract). Thus, small apparent mass mixing ratio observed by Sugimoto\net al. (2021a) is not a problem.\nAll that being said, we added a bit of extra content at the first mention around line 110,\nand then a short paragraph in the Discussion to expand on these topics - particularly the\nmass balance and the samples potential low likelihood of carrying shock signatures.\n“A similar mass balance argument may explain the lack of shock features\ndetected among samples from Ryugu so far. Nakamura et al. 2023 found that the\nphysical properties of Ryugu materials decays the impact-induced shockwave\nvery rapidly due primarily to the very high porosity of the material. This results in\nminimal volume being affected by the shock of each large impact and a lower\nprobability of any sample carrying the signature of the entire collisional history of\nthe asteroid.”\nAt line #485, the authors have extended the same analysis used for C-type asteroids\nRyugu and Bennu to S-type asteroids Itokawa and Didymos. Nevertheless, there is\npresently no compelling evidence to suggest that the Q*D equations would be\napplicable to the less porous S-type asteroids. Therefore, additional explanation is\nneeded or the paragraph removed altogether.\nThis is an interesting point, and we understand the concern. Unfortunately, to our\nknowledge, there also is no definitive data yet available that shows that S-type asteroids\nhave a different Q* than C-types for asteroids the sizes of Bennu/Ryugu. The limited\nD\ninformation we have instead suggests they may be similar.\nFor example, from the numerical hydrocode side of things, the closest work on this\nmight be the Jutzi et al. (2019) paper introducing and describing the P-alpha porosity\nmodel used for ~100km scale impact modeling, where they state this in their abstract:\n“Overall, the porous targets have a significantly higher impact strength (Q* ) than the\nD\nrubble-pile parent bodies investigated previously (Benavidez et al., 2012) and show a\nbehavior more similar to non-porous monolithic targets (Durda et al., 2007).” But this\nwork is specifically describing findings of pore-crushing effects on monolithic targets,\nwhereas with rubble piles much of the impact energy might be expended in simply\ngravitational dispersion of the target, which would seemingly level the playing field\nbetween different material properties.\nIndeed, Leinhardt & Stewart (2012) compiled data from catastrophic-disruption\nsimulations in the gravity regime from various literature sources and found that for (i) a\nvariety of target compositions, (ii) five orders of magnitude in size (~0.4–4e5 km) and\n(iii) nine orders of magnitude in impact energy, the Q* relationship could be described\nD\nby a narrow range in scaling constants. Therefore, the Q* disruption relationship would\nD\nbe applicable to gravity-dominated rubble piles regardless of their composition. Note\nthat Q* relationships describe a specific impact energy, so a C-type would require less\nD\nenergy to catastrophically disrupt than a S-type of the same size, due to purported\ndifferences in their bulk densities. On the other hand, C-types have a larger porosity\nthan S-types, and simulations by Martin Jutzi and others suggest this could make them\nmore difficult to disrupt. The two effects combined seem to compensate for one\nanother, at least in hydrocode impact simulations, and this suggests that S- and C-types\nmay have similar Q* functions.\nD\nA second example would be to consider cratering among main belt asteroids, which is\nrelated to Q* . Ideally, the production of a big enough crater on a world should disrupt\nD\nthe target, so cratering and disruptions have to be self consistent with each other.\nBottke et al. (2020) modeled the collisional evolution of the main belt size frequency\ndistribution (SFD) using a 1-D collisional evolution code and derived Q* function that\nD\ncould match constraints. Given that ~ 80% of the main belt is C-complex, it makes\nsense that their Q* would be most suitable to those bodies. What they found is that\nD\ntheir derived SFD did a good job reproducing the observed crater SFDs on all D > 10\nkm asteroids, including those which were S-types (e.g., Eros, Gaspra, Ida), C-types\n(e.g., Ceres, Mathilde), and other types (e.g., M-type Lutetia, V-type Vesta). So, if big\nQ* differences existed between different asteroid types, one would probably expect to\nD\nsee a modeling failure when attempting to match crater SFDs on different kinds of\nbodies.\nThe bottom line is that we believe there is no simple “correct” answer in the literature to\naddress this adequately, so instead we made an effort to more strongly caveat our\nfindings on this point. Part of the additional caveats and language added is the\nclarification that essentially all of the targets in this model are rubble piles as opposed to\nmonoliths.\nWe have also added an entirely new paragraph to the Discussion to answer a different\nreviewer’s comment, explaining that the results here roughly scale with collisional\nlifetime - so a reader now has this additional context regarding lifetimes and Q*\nD\nformulations.\nRef:\nJutzi, M., et al. (2019) Icarus 317, 215.\nBottke, W. F. and 14 colleagues 2020. The Astronomical Journal 160.\n#102: Some references are needed\nAdded\n#186: There should be a comma after (Bottke et al. 2022), not full stop.\nFixed\n# There are a couple of repetitions throughout the manuscript. See Lines #111 & #351\nYes, this was a serious repetition that has been corrected. We cut this sentence on its\nsecond appearance:\n”Thisismuchyoungerthanthe ageof thetwopotential sourcefamiliesandyoungerthan the\ncollisional lifetime for asteroidstheir size(~200Myr and~440 Myr respectively).“\nReviewer #2 (Remarks to the Author):\nThe manuscript \"Bennu and Ryugu likely 2nd or later generation rubble piles\"\ninvestigates a scenario in which both asteroids are 2nd or later generation fragments\nfrom an inner asteroid belt body. The study is motivated by apparently very different\nages of asteroid families where the asteroids are supposed to originate from and the\ncrater retention ages.\nTo this purpose, the authors employed a novel, Monte-Carlo-based approach that\nsimultaneously treats the collisional and dynamical evolution. This is one of the first\nworks in this respect and is appropriate for the problem studied.\nIndeed, the model is based on different assumptions and approximations, but the\nauthors also investigated how sensitive the results are concerning some of the input\nparameters. For instance, the authors show that the results do not change significantly\nfor different slopes of the size-frequency distribution nor when varying the sizes of the\nparent bodies.\nThe manuscript is well-written, and the figures and tables are appropriate and\ninformative.\nTherefore, it fully deserves publication.\nThere are only minor issues to be addressed:\nGenerally, if I'm not missing something, it is not described when exactly a test particle\nis assumed to be entering the NEAs region, i.e., when it is considered to be 'delivered'.\nThe nu6 and 3/1 resonances have been supposed to remove 100% of objects, right?\nBut what happens to them? I assume that a fraction from each resonance is considered\nto be delivered into the NEA region, but I cannot find precise information. This aspect of\nthe model should be better described, and missing details should be provided.\nThanks for this - yes, it should have been more clear in the manuscript. Fundamentally,\nthe setup of the model allows an indifference to the absolute or relative delivery rate\nfrom the respective resonances. The works that link the two asteroids to these two\nfamilies (Eulalia and Polana; e.g. Bottke2015) already account for that level of detail\nand the correct dynamical mixture of the near-Earth population. Here, we simply get to\nuse that knowledge of the dynamical pathway (nu6) and determine that if an asteroid\ncomes from the nu6, and its parent was Polana….then *this* is its likely collisional\nhistory. We don’t need to compare it to all of the other NEAs from all of the dynamical\npathways because we know where it came from already. Here we are just reporting on\nthe likely collisional history of simulated objects that we model to follow the same\ndynamical path. We did not establish this clearly in the original manuscript, but have\nmade the following edits - a whole new paragraph on these important details:\n“Furthermore, while the nu6 and 3:1 are effectively 100% efficient at removing asteroids\nfrom the asteroid belt, neither are as efficient at placing them onto near-Earth orbits, as\nsome are ejected from the Solar System or achieve some other fate. Here, since previous\nmodeling has found Bennu and Ryugu’s likely dynamical pathway to their current orbits,\nthere is no need to blend the relative probability of NEO delivery between these two\nremoval sources. Instead, we are exploring the collisional history of objects that follow\nthis same dynamical path. Since this effort relies on an Monte Carlo approach the total\nnumber of initial simulated asteroids in a family can be as large as computationally\nfeasible as long as it adheres to the family size distribution. The numbers and sizes of\nfragments made in each subsequent collision and disruption must follow the same rules\nand make proportionally similar numbers, while still adhering to the family size\ndistribution. Thus, nominal simulations below begin with 300,000 simulated asteroids,\nalthough the number grows throughout each simulation depending on the number and\nnature of the collisions in each.“\nline 41: \"All but the largest NEAs are expected to be rubble piles\"\nThis sounds like a too bold statement. Are you claiming that also all objects below ~100\nm are rubble piles? Even super-fast rotators? Please rephrase and/or be more specific.\nGood point - we are definitely NOT speaking about the D<200m objects, so we will\nrevise.\nIt now reads: “NEAs larger than a few hundred meters and smaller than 5-10km are\nexpected to be rubble piles, which are the reaccumulated debris from the collisional\nfragmentation of larger asteroids “\nline 54: \"Models of NEA delivery find that the probability of a body\non Bennu's current orbit has a 82% and a 18% probability to originate from the ν6 or\nIMC, respectively\"\nIt would be helpful to provide here also the probabilities and routes from the model by\nGranvink et al. (Icarus, 2018)\nHere is what we find on the Lowell website that hosts the data from the Granvik works\n(Moskovitz et al. 2022 - now referenced in the paper):\nBennu: 86% nu6 and 3% from Hungaria\nRyugu: 88% and 5%\nDidymos: 82% and 7%\nItokawa: 88% and 4%\nAnd here is the text that was added:\n“Updated NEO models, using a larger number of dynamica sources find similarly\nhigh probabilities of delivery from ν resonance where Bennu is 86% likely and\nline 61: Nesvorny et al. 2015 is definitely an appropriate reference for asteroid families.\nYet, it might be good to cite along a more recent review by Novakovic et al. CeMDA,\n2022, which is some aspect complementary.\nThanks for this. We had not seen the new review - and its good to see because a fair bit\nhas happened since 2015 that is now included.\nline 157, Figure 1: The reference for the family ages should be added, Bottke et al.\n(Icarus, 2015)\nAdded.\nline190: \"There is some initial dispersal of orbital semi-major axis that scales with the\nescape velocity of the target bodies (see Vokrouhlicky et al. 2006), but here the very\nlong timescale of their evolution overwhelms the initial scatter of initial orbital\nsemi-major axis, so each family member starts with identical semi-major axes.\"\nI agree that the initial dispersions should play a minor role compared to post-impact\nevolution. It seems, however, that it would be simple to implement this aspect into the\nmodel. At least for larger parents, it might be relevant because the expected ejection\nvelocities are higher. Also, in the case of Eulalia, which is close to the 3/1 resonance,\nadding the initial dispersion might place some bodies on the other side of the\nresonance.\nThe reviewer was totally right here, the initial velocity dispersion from the initial family\nbreakup makes a difference. We implemented this and the re-run for the nominal\nPolana case with initial velocity dispersion gets a crossover of 450Myr, instead of\n470Myr that was reported for not having velocity dispersion. That is … pretty small but\nthe Eulalia numbers drops a bit more down to 510Myr, and that is significant at the\n~10% level. We re-ran everything and updated all the plots and numbers throughout the\ntext. The big picture takeaway is not changed, but some key numbers definitely\nchanged enough to make this worthwhile.\nWe added a bit to the text:\n“There is some initial size-dependent dispersal of orbital semi-major axis that scales\nwith the escape velocity of the target bodies and this is included for each initial\nfamily and subsequent sub-family following the formulations of Vokrouhlicky et\nal. (2006). “\nAnd later to clarify that each of the sub-families, uses the same formulations:\n“Similarly, no bodies smaller than 0.5 km are added. Newly formed bodies are started\nnear the location of their recently disrupted “parent”, whose size sets the scaling for\nthe the size-dependent spread in the semi-major axis relative to the parent’s “\nline 244: \"Note that the removal efficiency is drift rate-dependent, which itself is size\ndependent (Bottke et al. 2015) and these removal calculations are valid only for the\nkm-scale objects.\"\nMilic Zitnik & Novakovic (ApJ, 2016) studied the effect of some MMRs on the\nsemi-major axis drift due to the Yarkovsky effect. Though it may not be directly relevant\nto this work, it could be worth mentioning it here as well, as, in general, this analysis\ncould be very useful in the model like the one presented in the manuscript.\nThanks for this - reference added.\nline 379, Figure 6: This figure's caption says there were 300,000 initial asteroids. Was\nthat the number of test particles in all simulations? Please specify. Also, this number\nshould be stated earlier when the simulation set-up is described.\nGood point, we should have been more detailed in the code setup description. We have\nadded this in the main text in an extra paragraph that clears up a handful of missing\ndetails about the modeling, including this.\n“Furthermore, while the nu6 and 3:1 are effectively 100% efficient at removing\nasteroids from the asteroid belt, neither are as efficient at placing them onto\nnear-Earth orbits, as some are ejected from the Solar System or achieve some\nother fate. Here, since previous modeling has found Bennu and Ryugu’s likely\ndynamical pathway to their current orbits, there is no need to blend the relative\nprobability of NEO delivery between these two removal sources. Instead, we are\nexploring the collisional history of objects that follow this same dynamical path.\nSince this effort relies on an Monte Carlo approach the total number of initial\nsimulated asteroids in a family can be as large as computationally feasible as\nlong as it adheres to the family size distribution. The numbers and sizes of\nfragments made in each subsequent collision and disruption must follow the\nsame rules and make proportionally similar numbers, while still adhering to the\nfamily size distribution. Thus, nominal simulations below begin with 300,000\nsimulated asteroids, although the number grows throughout each simulation\ndepending on the number and nature of the collisions in each “\nline 485: \"Flora has a similar age to New Polana...\"\nPlease add a reference for the age of the Flora family, Vokrouhlicky et al. (AJ, 2017)\nAdded.\nline 513: In the Yarkovsky scaling relation, it seems that the scaling with the obliquity of\nBennu is missing. Or is it omitted because it is close to 180 deg, implying that cosine is\nclose to -1, and the sign is neglected? Please check and fix or add missing details.\nIts the latter - Bennu is very close to 180deg, so its built-in to the formulations for the\nreasons you say. We added a short phrase telling the reader this - since its important!\nReviewer #3 (Remarks to the Author):\nMain belt asteroids can be disrupted by catastrophic collisions, leading to the formation\nof fragments. These fragments are referred to as a \"family\" of the original asteroid. The\norbits of these fragments, now smaller asteroids, can be altered by the Yarkovsky effect.\nAdditionally, these fragments can be further disrupted by subsequent collisions with\nother small bodies. When a main belt asteroid fragment reaches a certain radial\nlocation, known as an orbital resonance, it can be dynamically transferred to the\nnear-Earth region.\nThis study utilized numerical simulations to examine the orbit evolution of small bodies\ndriven by the Yarkovsky effect. It combined this with the catastrophic collisional\nevolution using a Monte-Carlo approach. The authors aimed to investigate specific\nasteroid families, such as New Polana and Eulalia, and to discuss their collisional and\norbital evolution. Recently, spacecrafts have visited near-Earth asteroids Ryugu and\nBennu. These asteroids are thought to be members of the collisional families of New\nPolana or Eulalia. However, the surface ages of these asteroids, as reported, are\nsignificantly younger than the estimated age of their parent family's formation. Through\nnumerical simulations, the authors determined that Ryugu and Bennu are not fragments\nfrom the initial family-forming collision. Instead, they are products of subsequent,\nadditional collisions.\nIn my assessment, the findings presented are both intriguing and relevant to the current\nscientific context. The manuscript is largely well-written and the content is effectively\nconveyed. As such, I believe it meets the standards for publication in Nature\nCommunications. However, before I can recommend its acceptance, I would like the\nauthors to provide further clarity on their numerical models. Additionally, I urge them to\ndiscuss the consistency of their findings with the returned samples from Ryugu — and\nfrom Bennu, if available. Addressing these points is crucial for validating and bolstering\nthe credibility of their numerical results.\nMajor comments.\n(1) In lines 125-130: I would appreciate it if the authors could first verify whether their\nnew approach — encompassing both orbital and collisional evolution — are able to\naccurately reproduce the Size-Frequency Distribution (SFD) of asteroids within the main\nasteroid belt. Additionally, could the authors discuss any updates or differences in their\nfindings, especially in comparison to the results presented in Bottke et al. 2005 and\nMorbidelli et al. 2009?\nWe understand the question/desire here, but find this to be well modeled in previous\nworks like Bottke et al. 2005 and Morbidelli et al. 2008 that consider the collisional\nevolution of the entire asteroid belt over very long timescales. Such models necessarily\nconsider all of the collisions for all of the asteroid belt that will globally average out to\nmatch the global asteroid belt SFD. Here, we are considering specific, singular, families\nthat represent a very small fraction of the asteroid belt. This is actually why we explored\nsuch a wide range of family SFD - to bracket a huge range of possibilities and show that\nthe selection of SFD is largely unimportant to the specific outcome here. So, we will add\nmore text explaining why we didn’t try to take on the entire asteroid belt all at once, but\nprovide here for the reviewer a more verbose explanation.\nThe entire belt is a mixture of numerous families of different size and age (see Broz et\nal. 2013, or Nesvorny et al. 2015). The entire catalog of families is uncertain, especially\nat older ages, and the sizes and ages of the clearly identified families are also both\nuncertain. The combined asteroid belt SFD is a combination of all of these families of\ndiffering size and age and all background asteroids that may not have any history with\nany family. The impact simulations of Durda 2004/2007 and Michel 2001/2003 find that\nan initial family SFD varies with the different possible impact properties (impact energy\nand angle etc.), and collisional modeling of the belt over time (Bottke 2005, see also\nVokrouhlicky et al. 2006 for an amazingly detailed work) find that a families’ SFD will\nchange over time due to the collisional evolution imparted by its interaction with the\nneighboring asteroid belt. So, collision modeling tells us that its hard or impossible to a\npriori know a families’ initial SFD and we also know their SFD will change over time. We\nalso don't know precisely when any family formed, or precisely how large it was when it\ndid form (this is why we explored such a wide range in possible parent body sizes). So,\ntrying to model the model belt over its whole history is a big task with too many\nuncertain parameters. This is exactly why we did this family-by-family with such wide,\nenveloping, ranges of parameters - to be sure that we bracketed the range out possible\noutcomes.\nThe point is well-taken though - we need to be more explicit about the motivation for the\napproach that we took. We have added text to be more clear and to more clearly\ndifferentiate this approach from the Bottke/Morbidelli SFD modeling approaches.\n“Other models have developed sophisticated routines to capture the effects of collisional\nevolution, and its large-scale effects on the size frequency distribution (SFD) of the\nasteroids in the main asteroid belt (Bottke et al. 2005, Morbidelli et al. 2009). These\nlarge models cover such large timescales and large swaths of the Main Belt that\nthey don’t take into account individual asteroid families and don’t incorporate the\nability to track individual asteroids and their changing orbits. “\nRefs:\nBroz M., Morbidelli, A., Bottke, W.~F., Rozehnal, J., Vokrouhlicky, D., Nesvorny, D.\n2013. Astronomy and Astrophysics 551.\nDurda, D.~D. and 6 colleagues 2007. Icarus 186, 498–516.\nDurda, D.~D. and 6 colleagues 2004. Icarus 170, 243–257.\nMichel, P., Benz, W., Tanga, P., Richardson, D.~C. 2001. Science 294, 1696–1700.\nMichel, P., Benz, W., Richardson, D.~C. 2003. Nature 421, 608–611.\nVokrouhlicky D., Broz, M., Bottke, W.~F., Nesvorny, D., Morbidelli, A. 2006. Icarus 182,\n118–142.\n(2) In line 415, the authors note that collisional evolution might significantly modify the\nparent body materials, leading to a mixture of materials from both the impactor and the\ntarget. However, recent literature, specifically Nakamura T. et al. 2022 in Science, has\nindicated that sample material from Ryugu doesn't show evidence of impact shock.\nFurthermore, Tatsumi et al. 2021 reported that the fraction of exogenic materials in the\nRyugu sample is comparatively small. One would anticipate that multiple catastrophic\nimpacts would both shock the materials and promote mixing between the impactor and\ntarget. This presumption seems to contradict the sample analysis results from Ryugu,\nas mentioned earlier. Could the authors provide an explanation for this apparent\ndiscrepancy?\nIn line with the aforementioned concerns, I request the authors to provide specific\ndetails regarding the impact conditions encountered during the collisional evolutions.\nCan the authors elucidate on the impact velocity and the impactor-to-target size ratio as\nderived from their simulations? Following this, I would appreciate it if the authors could\nestimate the consequent increase in pressure and temperature resulting from the\nspecified impact conditions. I hypothesize that if the impactor-to-target size ratio is close\nto unity, it would likely lead to enhanced material mixing. Could the authors comment on\nthis?\nTo the first point - this work is meant to set the baseline expectation for collisional\nhistory, from which these other findings are interpreted. We have a more verbose\nresponse to a similar question to referee #1, and we did add a paragraph to the\ndiscussion of the paper describing the Nakamura findings and some implications - we\nthink these explain this discrepancy.\nTo the second point, specific impact conditions are beyond the scope of this modeling\neffort for a few reasons that we try to describe here and with some more text in the\nmanuscript. Primarily, the specific impact conditions cannot be extracted in any\nmeaningful way because this work relies on probabilistic calculations of the collisional\nlifetime, which was derived in the works of Bottke et al. (2005) and Ballouz et al. (2020)\nthat take into account the collisional probabilities and impact conditions etc. For us, in\neffect at each timestep we are rolling the dice to see if each asteroid has an impact that\nis as energetic as Q* , if so, we break it up and move on. We are not modeling crossing\nD\norbits and tracking impact angles or impactor sizes and velocities. That being said,\nimpact speeds are relatively uniform across the asteroid belt in the 5-6km/s range,\nimpact angles fill the allowable distribution. The impactors required to reach Q* level\nD\ncollisions are relatively small compared to the size of the target, and larger ratio\nimpactors are increasingly unlikely.\nThe key paragraph where we discuss this in the Discussion, and we have added more\ntext to better address these questions, is here;\n“One of the simplifying assumptions of this model is that the only meaningful collisions\nare precisely at the energy where half of the mass is liberated from the target (Q* ). “,\nD\nand we have added this to the end of that paragraph: “While collisions substantially\nmore energetic than Q* are rare events, they could provide pathways to more\nD\nwidespread shock and mixing among the target and its reaccumulated remnants.\n“\n(3) In lines 185-190: The authors have posited that collisions with background asteroids\nprimarily influence family members. However, immediately following the initial\ncatastrophic collision that results in the formation of an asteroid family, the fragments\nwould likely have orbits in close proximity to one another, given that they would share\nan almost identical point of impact. Consequently, right after the impact, there's a\nplausible scenario where these fragments (or family members) might collide with each\nother frequently, leading to further fragmentation. Could the authors provide a\ndiscussion or estimation on the collision timescale under such conditions?\nTechnically, our model can’t account for this as it does not calculate impact probabilities\n- they are incorporated via the collisional lifetime calculation. Note that the initial velocity\ndispersion of the fragments, implemented at the request of Reviewer #2, finds that\nkm-sized fragments can have initial orbits ~0.05au away from a 100km primary. So,\nwhile all fragments will briefly share an orbital node from the disruption event, they will\nhave a large spread of orbital periods and then the action of Yarkovsky and orbital\nprecession will drive them away really quickly from that shared node. This effect was\ndiscussed in Dell’Oro et al. (2002; Icarus). From their abstract:\n“The results show the occurrence of a strong enhancement in the mutual collision rate among\nfamily members, immediately after family formation. Nevertheless, this episode lasts for a\nrelatively short time, and it does not affect too severely the overall collisional evolution of the\nfamily. The early enhancement of mutual collisions, however, may influence the cratering record\nexhibited by the surfaces of family members, possibly laying the foundation for the early\ndevelopment of a surface regolith layer.”\nMeanwhile, Bottke et al. (2020; AJ) modeled of the crater size distribution on spacecraft\nobserved asteroids, which included family members like Gaspra (with Flora) and Ida\n(with Koronis). They found nothing anomalous in the crater history of either one that\nwould suggest a major contribution from family members. Until more evidence\nbecomes available, we would argue it is reasonable to neglect the effect of intra-family\nmember collisions on our test asteroids.\nMinor comments:\nHirabayashi et al. 2020; Hyodo & Sugiura 2022 also discuss the shape change and\nresurfacing of small asteroids.\nYes, this is totally reasonable - we know YORP works out there in the Main Belt. We\nadded this (I paste the whole paragraph for context):\n“Similarly, YORP spinup and re-shaping has been attributed as the likely source of the\ndistinct top-shape that both Bennu and Ryugu have and is common among NEAs and\nespecially binary asteroids which account for 15% of the NEA population (Margot et al.\n2015). YORP spinup is possibly a key or driving element of the spectral changes\namongst S-type asteroids where movement of surface material is needed to liberate or\nexpose fresh material that has not been exposed to space weather (see Graves et al.\n2018). Therefore, YORP spinup may provide a pathway to re-surface via landslides\nand re-shaping on timescales shorter than collisional lifetimes (Walsh et al. 2008,\nHarris et al. 2009, Scheeres 2015, Hirabayashi et al. 2020, Hyodo and Sugiura\n2022). The YORP effect may also combine with the effect of tides during close flybys of\nthe terrestrial planets or the effects of seismic shaking due to small impacts, to further\ncontribute to the erasure of some cratering history on small asteroids (Richardson et al.\n1998, Richardson et al. 2020, Jawin et al. 2020;2022).”\nREVIEWERS' COMMENTS\nReviewer #1 (Remarks to the Author):\nI believe the authors have appropriately addressed my concerns and I have no further comments.\nReviewer #2 (Remarks to the Author):\nI've read Walsh et al.'s revised version of the manuscript \"Bennu and Ryugu likely 2nd or later\ngeneration rubble piles.\" I appreciate the work the authors did in response to my comments. They\nadequately addressed the concerns raised in the first report and substantially improved the\nmanuscript. Therefore, I recommend it for publication in Nature Communications.\nReviewer #3 (Remarks to the Author):\nThe authors have answered all my questions and concerns. I think the manuscript is now worthy of\npublication in Nature Communication."}]

Decision letter