Dynamical Evolution of Simulated Particles Ejected From Asteroid Bennu

In early 2019, the OSIRIS‐REx spacecraft discovered small particles being ejected from the surface of the near‐Earth asteroid Bennu.sww Although they were seen to be ejected at slow speeds, on the order of tens of cm/s, a number of particles were surprisingly seen to orbit for multiple revolutions and days, which requires a dynamical mechanism to quickly and substantially modify the orbit to prevent re‐impact upon their first periapse passage. This paper demonstrates that, based on simulations constrained by the conditions of the observed events, the combined effects of gravity, solar radiation pressure, and thermal radiation pressure from Bennu can produce many sustained orbits for ejected particles. Furthermore, the simulated populations exhibit two interesting phenomena that could play an important role in the geophysical evolution of bodies such as Bennu. First, small particles (<1 cm radius) are preferentially removed from the system, which could lead to a deficit of such particles on the surface. Second, re‐impacting particles preferentially land near or on the equatorial bulge of Bennu. Over time, this can lead to crater in‐filling and growth of the equatorial radius without requiring landslides.


Introduction
The OSIRIS-REx spacecraft arrived at the near-Earth asteroid Bennu in late 2018 . In early 2019, particles were discovered being ejected from the surface of Bennu . One surprise was the length of the lifetimes of several of the observed particles, whose orbits were estimated to last multiple days and complete many revolutions -demonstrating that some fraction of the ejected particles were put into orbits that neither immediately re-impacted the surface nor immediately escaped the system. These observations brought up many questions. What dynamical processes could lead to such orbits? How do particles launched at relatively slow speeds avoid the fate of re-impacting the surface as they come back down toward their first periapse passage? How long can ejected particles stay in orbit around Bennu? When ejected particles do re-impact, where do they land? This paper addresses these questions.
Bennu is a small near-Earth asteroid, approximately 500 min diameter, with a rubble-pile structure, a rocky surface, and a "top" shaped profile with an equatorial bulge (Barnouin et al., 2019;DellaGiustina et al., 2019;Scheeres et al., 2019). The dynamical environment of Bennu is complex due to the low gravity and non-spherical shape of this small body (Scheeres et al., 2019). This means that orbits in proximity of the body are highly perturbed by solar radiation pressure (SRP) forces and are non-Keplerian and rapidly evolving in general (Scheeres, 2016).
Most studies of orbits about small asteroids focus on stable orbits that will be useful for spacecraft exploring such bodies. Scheeres (2016) has developed an averaged theory that succinctly describes the evolution of orbits around small bodies when they are perturbed by SRP. He shows the existence of frozen orbits and stable terminator orbits, which have now been successfully flown by the OSIRIS-REx spacecraft . Many studies have advanced this work, finding specific types of orbits that exist under the SRP and solar gravity perturbations, including quasi-terminator orbits (Broschart et al., 2014), heliotropic The leading hypotheses for the cause of the observed ejection events at Bennu are thermal fracturing or micrometeorite impacts, either of which could lead to the relatively low energy ejecta seen at Bennu . There has been a significant amount of work investigating ejecta from natural and man-made impacts on small asteroids. Unfortunately, these impacts take place at high energies, meaning that much of the ejecta is at higher speeds than is of concern here. However, a few studies have looked at the low-velocity portion of the ejecta population. The understanding of the fate of impact ejecta at asteroids is discussed by Scheeres et al. (2002), which points out how ejecta at small asteroids can, in theory, enter into orbits under the effects of gravity and SRP. Specific studies of ejecta at Ida (Geissler et al., 1996) and Eros (Korycansky & Asphaug, 2004) provide interesting comparisons to the current case; however, those bodies are an order of magnitude larger than Bennu, and the dynamics are therefore more strongly dominated by gravity. Furthermore, statistical results from those studies would not directly apply here because they are conditioned on initial ejecta populations created from high-energy impacts.
Similarly, there have been studies of the fate of ejecta and debris from man-made impacts on asteroids.
In particular, studies of the expected evolution of the debris cloud after the impact of the DART mission (Schwartz et al., 2016;Yu et al., 2017;Yu & Michel, 2018) and the Hayabusa2 Small Carryon Impactor experiment (Arakawa et al., 2017;Giancotti et al., 2014) have been carried out recently. While these asteroids are more similar in size to Bennu, the source of the ejecta is again from a high-energy impact, which differs from the observed events at Bennu because they predominantly produce high-velocity ejecta that quickly escape the system. A recent study by Vetrisano et al. (2016) has provided the closest study of low-speed ejecta from a small body to predict the events at Bennu. The fate of the ejecta is strongly controlled by the effects of SRP, which has also been found by Garcia Yarnoz et al. (2014). While these works are relevant and provide valuable insight, it is crucial to include two other effects to get realistic results, especially for low-altitude particles: shadowing from the primary body which turns off SRP when eclipsed (Russell et al., 2016), and thermal radiation pressure forces from the infrared radiation leaving Bennu (Hesar et al., 2017). This paper investigates the evolutionary outcomes of populations of simulated particles ejected from the surface under conditions similar to those observed at Bennu. Such an analysis provides insight into how ejection events can influence the distribution of material over the surface of the asteroid. The results presented here are constrained by the ejection events observed in early 2019 . The estimated ejection locations, timing, and velocity ranges from  are used in this paper, as well as representative particle sizes and masses that encompass the best available data. Having said that, the point of this paper is not to produce true or estimated orbits; based on our knowledge of the particle dynamics and Bennu's properties, that can only be done reliably with trajectories estimated from observations . Rather, this paper explores the influence of the parameters of the dynamical system and the particle initial conditions to understand the larger issues regarding how particles could move around in this system. We seek to balance the accuracy of the dynamics with computational speed, given the uncertainties still in the models (e.g., from gravity, albedo, and unmodeled dynamics), to keep computational speed tractable such that we can produce large numbers of simulations to understand the trends within a population of ejected particles. Thus, the real value in the results presented here is in the range of behaviors that can result from an ejection event. The population evolution that we simulate indicates that if ejection events occur often enough, they can play an important role in the geophysical properties of Bennu.

Dynamic Modeling
Effects that are typically thought of as small perturbations from the perspective of classical astrodynamics around planets become extremely important around small bodies owing to the weak gravity. The dynamics considered in this work are shape model-based gravity, solar tides, SRP including shadowing, and shape model-based thermal/albedo radiation pressure. In the course of this work and previous studies, it is found

Gravity
Although small, the main source of orbital dynamics is still the gravitational forces caused by the asteroid. We use the constant density polyhedral gravity model (Werner & Scheeres, 1996) to simulate the gravity field from the v20 Bennu shape model constructed from data obtained by the OSIRIS-REx spacecraft (Barnouin et al., 2019) and the estimated Bennu density of 1.19 g/cc (Scheeres et al., 2019). Particular parameters used for these models are given in section 4.
The top shape of Bennu produces a gravity field that is primarily dominated by the even zonal harmonics, especially J 2 and J 4 (McMahon et al., 2018). The body is relatively symmetric at a global scale with respect to the pole and the equator (Barnouin et al., 2019), meaning that the odd zonal and tesseral harmonics are less significant but do exist and are captured by the polyhedral gravity model. Because Bennu does not exhibit any significant wobble in its rotational pole (Barnouin et al., 2019;, the main effect of the non-spherical gravity potential is to precess a particle orbit's angular momentum and eccentricity vectors (Scheeres, 2016).
The other important gravitational effect which must be considered is the effect of solar tides, which are modeled as where Sun is the gravitational parameter of the Sun, r Sun/p is the vector pointing from the particle to the Sun, and r Sun/Ast is the vector pointing from the center of the asteroid to the Sun. Solar tides will also primarily have the effect of torquing a particle's orbit to precess the angular momentum and eccentricity vectors. On a longer timescale, the solar tides can lead to the Kozai effect trading inclination and eccentricity for non-equatorial orbits (Rieger et al., 2018); however, this secular effect is often interrupted for the particles considered in this work given the rapid evolution of orbits from the other dynamics acting in the system.

Solar Radiation Pressure
After gravity, SRP is the most important force acting on the ejected particles. The most widely used model for SRP is the so-called cannonball model, which captures the primary component of the acceleration in the anti-Sun direction. The particular version of the SRP model used here is shown in equation (2).
where H(r) is the shadowing function that takes a value of 0 if the particle is positioned (where r is the particle's position with respect to the asteroid) behind Bennu such that the Sun is occulted, and 1 otherwise. We do not model any partial shadowing/penumbra effects. In our code, and as shown in equation (2) we approximate the distance from the Sun to the particle as |r Sun/Ast |, as the difference between these is minimal. The same is true forr Sun∕Ast , which is the unit vector from Bennu (as opposed to the particle) to the Sun. The minus sign makes the SRP acceleration act in the anti-Sun direction. P 0 is the solar pressure constant, which has a value of 1 × 10 14 kg km/s 2 , is the reflectivity, or albedo, of the particles, and A∕m is the area-to-mass ratio. The simulations only use these values in ratio, although we do define the individual values from an assumed spherical shape of constant density (see section 4). The 4∕9 factor that appears with the reflectivity comes from the assumption that the particle is a sphere (on average) that reflects light in a diffuse Lambertian pattern.
It is important to understand the assumptions that are embedded in using this model for SRP. The name cannonball implies that the particles are spherical. This assumption is commonly used because an object of any shape, if it is tumbling, will experience an SRP acceleration away from the Sun on average. Specifically, if an object is tumbling such that (1) its rotational rate is much faster than the mean motion of the orbit and (2) there is an equal probability of the body being at any inertial attitude in time, then the SRP model will average out to being in the anti-Sun direction. The interpretation of the area-to-mass ratio being from a spherical particle of constant density is an easy way to compute realistic and representative area-to-mass ratios. Because the particles in reality could be closer to a tumbling plate-like shape (Rizk et al., 2019), the relationship of area-to-mass ratio to density and reflectivity should be taken with some uncertainty as it is an averaged dynamical quantity. Two further assumptions are embedded in this model: (1) Any reflected light is reflected in a purely diffuse Lambertian manner; (2) absorbed light that is re-emitted as infrared radiation cause any acceleration on the body because the small sizes and assumed tumbling motion lead to the particles being isothermal.
In order to produce realistic orbital evolution, it is crucial to include shadowing as represented by the H(r). This fundamentally changes the effects of SRP on an orbit. For example, without shadowing, SRP on average does not change the semimajor axis of the orbit. However, when shadowing is taken into account a change of semimajor axis can occur. The details of our implementation of a fast shadowing algorithm are discussed in section 3.

Thermal Radiation Pressure
Thermal radiation pressure (TRP) from the radiation emanating from the asteroid is generally much smaller than SRP. However, in this scenario, all particles necessarily spend time near the surface, where the TRP forces can approach or even exceed SRP. Therefore, it is crucial to include these forces in the dynamical models simulated.
The TRP model used is from Hesar et al. (2017) but simplified for a cannonball particle instead of a complex spacecraft shape as in that work. The acceleration can be computed as where the summation goes over the number of facets of the shape model, N F , whose positions are referenced on the body by the position of their centers, r i . There can be a reflection of the incident radiation based on an infrared albedo, ; however, we treat this parameter as zero in this work given the isothermal assumption discussed in section 2.2. P i is the infrared pressure coming from facet i, which is defined as where i is the visibility function of the surface element i with respect to the sunlight; that is, i is equal to 1 if that surface element is lit by the sunlight and 0 otherwise. Θ is the angle between the facet normal and the incident sunlight. Ast is the albedo of Bennu, which is defined as the fraction of the shortwave radiation reflected from the surface of the body to the incident shortwave solar radiation. Here we assume a constant albedo across the entire surface of the body of 4% . G R is the solar flux at the distance R = |r Sun/Ast | from the Sun (=1,368 J s −1 m 2 at 1AU), and c is the speed of light. A i is the surface area of facet i. is the angle between the facet normal and the vector connecting the particle and the facet center. This determines the visibility, and if < 0, this facet does not contribute to the total TRP at this time. is the surface emissivity of Bennu, and B is the Stefan Boltzmann constant.
T i is the temperature of the facet, which is determined by the Advanced Thermophysical Model (ATPM) of Rozitis and Green (2011, 2012 using the thermophysical properties of Bennu derived by DellaGiustina et al. (2019). The hottest region on the asteroid is in the mid-afternoon. The ATPM takes into account topography and thermal inertia effects such that the temperatures are not symmetric, and the TRP acceleration at a given location will vary with the spin state of Bennu. This variation shrinks as altitude increases such that it is insignificant by around 1 km, but at low altitudes the variation can be 5% to 10% of the total TRP. The temperature map is computed at one specific Bennu orbit distance, so that the temperature used is scaled by the relationship where R 0 and T i,0 are the distance to the Sun and the facet temperature at the epoch location, respectively. As with the SRP model, this model assumes that the particle is rapidly rotating such that its area-to-mass ratio averages to an effective constant value represented by the sphere in this work. The final term in equation (4) becomes extremely large as a particle approaches the surface such that |r − r i | → 0. This is not physical, but rather is an artifact of the discretization of the asteroid surface with finite facets. Thus, we implement a limit in our simulations such that A i ∕|r−r i | can never be larger than 1. Although this is not physically exact,

10.1029/2019JE006229
it captures the main behavior without requiring us to switch to a higher-resolution shape and temperature map, which would not significantly change the results.

Numerical Methods
The main simulation is written in Matlab, using the variable-step Runge-Kutta 45 integrator od4e5. This integrator performed well in this scenario once a normalization scheme was implemented to improve the numerics. The normalizing length is chosen to be the minimum radius of the shape model used, r = 214.68 m. This has the effect that a normalized position vector of length <1 is guaranteed to be inside the body. The normalizing time is then computed based on the mean motion at this distance, which is t = √r 3 ∕ = 1, 421.51 s, and the associated normalizing velocity is computed as the circular speed at the reference length, which is thenv = √ ∕r = 15.1 cm/s. This results in a normalized = 1. Using this normalization scheme allows us to use reasonable tolerances: a relative tolerance of 1 × 10 −3 and an absolute tolerance of 1 × 10 −6 .
Several other important components of the simulation implementation allow for fast execution. The polyhedral gravity mode, which is by far the most computationally complex portion of the dynamics, is coded in C and interfaced through a MEX function. The TRP model is written in Matlab but is formulated to take advantage of Matlab's sparse matrix capabilities to speed up the dot products that are computed for every facet of the shape model, which has produced a significant speed increase.
Finally, the shadowing model can be another computational bottleneck if ray-tracing is used. To avoid this, the shadowing algorithm is based on approximate limbs of Bennu represented by a convex hull defined by the maximum radius at every 12 • of latitude. This can then be represented with 30 pie-shaped triangular facets connected to the center of the shape model. This set of facets is used to check for shadowing and/or re-impact by projecting a particle's position vector onto the terminator plane and testing whether it resides within any of these facets; if so, then it can be determined whether it is in shadow or has impacted the body by looking at the total radius and comparing to the limb radius at that latitude. Our testing has shown that, while this approximation may be too rough for fitting precise measurement data, the dynamics produced do not differ meaningfully from a more precise model, and so the general trends presented in this work do not change substantially.

Ejection Event Simulation Parameters
The simulation results presented here are constrained by the measured quantities of Bennu and the particle ejection events. We investigate the evolution of particles based on the first three largest observed ejection events, which occurred on 6 January, 19 January, and 11 February 2019 . Various parameters used in the simulations are given in Table 1. The second and third events have well-estimated ejection locations on the body, which are used here. The 6 January event, however, has some uncertainty in the ejection location, which results in two possible ejection locations, which are referred to as Site A and Site B in this work (near and far solutions, respectively, in . Thus, we simulate four ejection events, one for each site/date combination as shown in Table 1. In this work the v20 shape model of Barnouin et al. (2019) was downsampled to a vertex spacing resolution of approximately 12.58 m with 12,288 facets and 6,534 vertices, which provides a good balance between accuracy for topography and gravity for a reasonable computational load. The radius for each event location in Table 1 is computed from where the indicated latitude and longitude intersect this shape model, so these values may differ slightly from reality at that location. The temperature model uses the same shape model resolution but is updated from the v13 shape model used in DellaGiustina et al. (2019) to the v20 shape model used here.
We made some approximations and assumptions to simplify certain aspects of the simulation without sacrificing the understanding of the general behavior of the ejected particles. First, particles are all modeled with reflectivity = 0.04, which is the mean Bennu albedo. Particles are modeled as spheres, such that the area-to-mass ratio varies as where d p is the particle density and r part is particle radius. In this work we used an assumed constant particle density of d p = 2 g/cm 3 which is similar to Bennu's bulk density and consistent with meteorite analogs (Hamilton et al., 2019;. This value is within the range of densities found in ; however, as discussed in section 2, the area-to-mass ratio controls the SRP and TRP accelerations, thus trading density and particle size can result in equivalent trajectories for different particle models. The SRP acceleration is also modified by the (4∕9) term in equation (2), which means that changing the reflectivity will also influence the dynamics, albeit with a weaker effect than the area-to-mass ratio. Overall, these values are based on the best information to date, but the population explored covers a range of area-to-mass ratios to try to encompass any expected variation.
Two other approximations are made to simplify the simulation environment. Bennu's spin pole is assumed to be perfectly retrograde with respect to its orbit angular momentum, when in fact there is a small obliquity difference (Barnouin et al., 2019). However, the maximum error in this assumption is only 2.55 • over Bennu's orbit (determined using the Bennu ephemeris and estimated pole available from the Osiris-rex naif repository, 2020); thus, this approximation should have only a small effect. Second, as discussed previously, the gravity is based on a constant density assumption with a finite-resolution shape model. While there are some indications that there is an inhomogeneous density distribution (Scheeres et al., 2019), the differences in the gravity field seen so far indicate that the constant density assumption is a reasonable first approximation, especially given that we do not know the true density distribution at this point. The same reasoning indicates that the chosen shape model gives a representative gravity field, especially at altitudes more than a few meters from the surface.

10.1029/2019JE006229
Given the above parameters, there are four degrees of freedom left to sample to simulate a population of ejected particles: the three dimensions of the launch velocity vector and the area-to-mass ratio. The launch velocity vector is the initial velocity vector with respect to the Bennu surface at which a particle is launched. The vector is parameterized by the magnitude and two directions: an azimuth angle measured from local East and an elevation angle measured from the plane of the shape model facet where the ejection event is located. The observations of the three ejection events show initial velocities ranging from 7 to 330 cm/s . In order to understand the possible orbital evolutions, we create populations of particles that sample all directions in the hemisphere above the ejection facet. The azimuth is simulated in discrete steps of 30 • , while the elevation is simulated in steps of 15 • . The velocities simulated range from 10 to 30 cm/s (note that all particles launched faster than 30 cm/s escape immediately, as shown below), in steps of 2 cm/s. Finally, to explore the area-to-mass ratio, the particle radius is varied from the set of 0.1, 0.5, and 1 to 20 cm. All told, this results in a grid of 11 velocities, 7 elevations, 12 azimuths, and 22 particle radii and area-to-mass ratios for a total of 17,666 simulations from each event/site (the azimuth does not come into play at an elevation of 90 • ).
Because the particle velocities are sampled from a Bennu-relative grid, the initial velocity used for simulation must be expressed in the inertial frame: This means that the initial inertial velocity will be skewed with an eastward component that grows in magnitude for sites closer to the equator of Bennu. Thus, westward (azimuth around 180 • ) cases can have initial inertial velocity magnitudes less than 10 cm/s, while eastward particles can be greater than 30 cm/s.

Results
Given the set of initial conditions and parameters discussed above, the 17,666 test particles were simulated for each of the four event times and locations (6 January, Site A; 6 January, Site B; 19 January; and 11 February). The following sections present some illustrative orbits to demonstrate the complex dynamical environment with the focus on understanding the general trends seen within the populations for all of the simulated scenarios. In cases where results for one scenario are representative of all simulated scenarios, we show only the results for one.

Orbit Evolution
The simulated particles demonstrate the rich, complex dynamical environment near Bennu. The non-Keplerian dynamics must quickly modify these orbits such that the particle will not impact the surface within its first revolution. Figure 1 shows the initial conditions from the grid discussed in section 4 mapped to a subset of the initial orbit elements. For every set of initial conditions in orbit element space (or in position/velocity space), there are 22 cases for the different particle sizes, as particle size does not change the initial state. In any given subset of initial conditions, there can be more cases at the same combination of initial conditions, as multiple launch velocities can lead to common orbit elements. Thus, Figure 1 does not intend to quantify the outcomes but indicates how the strongly non-Keplerian dynamics can result in very different evolutionary outcomes for the same or similar initial orbits.
Each simulated trajectory is grouped into one of four outcomes: suborbital, direct escape, escape, or orbital. A suborbital case is where the particle re-impacts the asteroid before passing through periapse, thus completing less than one revolution. A direct escape case is where the particle escapes the system before passing through periapse. Escape from the Bennu system is defined by a particle reaching a distance of 35 km from Bennu, which is roughly its Hill radius. An escape case is a particle that eventually escapes but first passes through one or more periapses. Finally, an orbital case is one which passes through one or more periapses and eventually either re-impacts with Bennu or, in a small number of cases, continues orbiting for a full Bennu year (437 days). These classifications roughly correspond to the classification proposed by Scheeres et al. (2002): Suborbitals are Class I; direct escapes are Class V; escapes are Class IV; and orbital cases encompass both Classes II and III.
There are several interesting conclusions to be drawn from Figure 1. First, many particles that are launched on what should be hyperbolic orbits (e > 1 and/or a < 0) do not escape immediately. Most escape eventually, but they often come back toward Bennu before escaping. These particles are usually launched toward the Sun, and SRP has enough time and strength to reverse the direction of motion such that the particles return toward Bennu and then fly by to a subsequent escape. Second, most particles that are launched with e < 1 are suborbital and do not make it past their first periapse; however, the orbital cases can begin with a wide variety of semimajor axes and very low periapse radii (all cases pictured have periapse radii less than the equatorial radius of Bennu)-indicating that the non-Keplerian dynamics can greatly change the trajectory to prevent immediate re-impact. The suborbital fate likewise dominates the low-energy (small a) trajectories, as would be expected. A third observation is that there are some cases where particles launched on very high trajectories (a ≃ ±20 km) enter orbit. These trajectories also typically move toward the Sun, which allows SRP to remove a significant portion of their orbital energy such that they can be in a lower energy state upon their first periapse passage.
To further demonstrate the non-Keplerian environment experienced by the particles, Figures 2 and 3 show time histories of the orbits and orbit elements for two particles that remained in orbit for the maximum simulation length of one Bennu year. These two particles had the same launch velocities-magnitude of 24 cm/s, azimuth of 150 • , and elevation of 45 • and differed only in their sizes, which were radii of 5 and 7 cm. The rapid variations in the orbit elements over the course of the year illustrate the complex dynamical environment.

Population Evolution
A grid study such as is presented here is best used to understand the general behavior of the overall populations of ejected particles. To this end, we wish to understand how the population for each ejection event evolves with time. It is of particular interest to understand what portions of the initial conditions lead to the four fates discussed in the previous section. This is pictured for one event in Figure 4; the other simulated events follow very similar trends. The population quickly drops with nearly half of the particles re-impacting the surface of Bennu within the first day, most of which are the suborbital cases. Interestingly, all direct escape cases last more than 1 day, meaning it takes at least that long for any particle to reach the Hill sphere. Most of the population has either re-impacted or escaped within 10 days. However, there is a small subset of the population which survives for much longer. Even a finite number of direct escape cases remain within the Hill sphere for nearly 50 days. In this set of 17,666 particles, approximately 20 particles survive between 50 and 437 days, with 10 particles still in orbit after one Bennu year.   One important aspect to understand about the particle lifetime is where the particles exist at a given time. To first order, Figure 5 answers this question by showing how many particles survive in a given radius band for the first 10 days after an ejection event. The population is grouped into three radii groups: <1 km, which is the near-surface environment; 1 to 5 km, which, for OSIRIS-REx, is of particular interest because this is where the spacecraft operates for most of the mission; and finally, 5+ km. The final line shown is the rest of the population, which has already returned to the surface or escaped. This plot is very similar for all four ejection events. It shows that the near-surface environment quickly loses most of its population, with less than 1% of particles spending time in this region after 1 day. The middle radius region also reaches 1% after around 2 days. More than 95% of particles re-impact or escape after 10 days. Finally, many simulated particles reside for long periods of time at high altitudes with respect to the asteroid; roughly half of the particles are beyond 5 km from Bennu 1 to 2 days after the ejection event, with many taking several more days to either escape or return to the surface. The population is not restricted to low altitudes. Figure 6 shows the relationship between launch velocity, area-to-mass ratio, particle size, and particle energy to the probability of escape. This figure demonstrates why we limited the grid search to be between 10 and 30 cm/s; all particles below 10 cm/s return to the surface, while all above 30 cm/s escape. Three main results can be drawn from these relationships. First, all particle sizes and area-to-mass ratios tested have a higher probability of escaping the system than re-impacting, but this is especially true for sub-centimeter particles. SRP can quickly add significant energy to these small particles, causing them to escape from lower initial velocities and energies. Second, and unsurprisingly, the latitude of the ejection event site plays an important Figure 6. Percentage of population that escapes the system as a function of (a) particle radius/area-to-mass ratio, (b) launch velocity, and (c) launch kinetic energy. role in the chance of escape; the lower-latitude events provide more velocity to the particles from Bennu's spin, and thus, particles at lower launch velocities can escape, but also those at higher velocities launched westward move slower and do not escape. Third, the relationship with launch energy is interesting because there is a sweet spot in terms of maximizing the chance to re-impact. The lowest energies are associated with the smallest particles (due to their small mass), and thus, they predominantly escape, while the largest energies also mostly escape due to the fact that they are launched at the highest velocities. In between, the interplay between mass, velocity, and launch geometry makes for a non-monotonic relationship.

Mass Migration
From a geophysical perspective, the most important aspect of the dynamics of ejected particles pertains to the particles that re-impact the surface. Where do they go? Is their distribution random? We gain insight into these questions through mapping the simulated re-impact locations from the four ejection event scenarios that we modeled, as shown in Figures 7 and 8, where the re-impact locations for each event are binned by latitude and longitude over the surface of Bennu.
The highest concentration location in each case is roughly west of the launch sites. This corresponds to a large number of suborbital particles that do not leave the surface for very long, simply letting Bennu rotate under them for some period before coming back to the surface. Not all suborbital particles follow this pattern, however, as some can reach high altitudes above the surface before coming back down, allowing much more movement. Next, in terms of longitude, although each individual event displays some preferences, the pattern is not systematic across all event scenarios tested. This makes sense: As with the high suborbital cases, the particles that enter orbit for a finite period of time can have their orbits drastically changed, and, along with the variable lifetime, this allows these particles to land at random longitudes. It is noted that there are not strong patterns in terms of the local time at landing, other than the fact that the short period suborbital particles land within a few hours of the ejection local time. Longer lived particles can land at a random local time given their assorted longitudes and lifetimes.
Latitude, however, is different. There is clearly an overall excess of ejection conditions that lead to re-impact at low latitudes. The 19 January and 11 February cases show a strong concentration near the equator. The 6 January cases are not concentrated as strongly near the poles but still show a bias in landing locations at lower latitudes than their launch locations. This can be explained by the shape of Bennu, whose radius is largest near the equator and tapers toward the poles, and therefore has a higher chance of catching a particle at a low portion of its trajectory in this region. Overall, the re-impacting particles appear to be migrating toward the equator.
The results shown in Figures 7 and 8 were totaled over all launch conditions to obtain a global view of the outcomes from a uniform ejection event. However, given the uncertainty surrounding the detailed physics of the ejection process creating the initial velocities , there could be a preferential Figure 9. Map of the re-impact locations for the 6 January, Site A launch site case for the azimuthal direction sensitivity study, along with the associated latitude and longitude histograms. The sketch indicates how the four azimuth cases are determined by projecting the launch velocity into the facet plane-in this example this case falls within the east grouping. direction of launch. To initially investigate this, we study two cases: an azimuthal preference versus an elevation preference for the launch velocity.
In the azimuthal study, the launch velocity directions are defined in cones, such that all initial velocities projected onto the facet are within ±45 • of the local cardinal direction included in that case-north, south, east, or west. The results of this study for one ejection event are shown in Figures 9 and 10. We note a longitudinal preference in re-impact locations between the different cases, with the East and North cases favoring a westward location, the West cases moving even further westward to include the opposite side of the body, and the South cases wrapping around and covering the eastward motion. We again see a trend of particles moving to lower latitudes-while this may be expected for such a high-latitude launch site, it was already shown in Figure 8 that lower-latitude launch sites are even more strongly biased toward low-latitude landings. This result is interesting because regardless of the direction, much of the material ends up downhill of the ejection site, even if it does not reach the equator (see Scheeres et al., 2019, for details of Bennu's low-latitude region being at a lower potential than higher latitudes). It is also noteworthy in Figures 9 and 10 that the eastward cases appear to follow a ground-track-type pattern with a maximum latitude around that of the launch site, which reinforces the fact that cases launched to the East are more likely to enter orbits that precess for some period before re-impacting than those launched in other directions.
In the elevation study, the cases are put into three bins: near horizontal (< [30 • ), near vertical ( >] 60 • ), and middle elevation between those two. Results for the 11 February case are shown in Figures 11 and 12. Here we see that the near-vertical cases move the least in longitude, while the near-horizontal cases move the farthest. All three cases show a fairly strong bias toward landing near the equator, which is partly due to this ejection event starting near the equator. However, events starting in this region do not show a preference for migrating to higher latitudes.

Discussion
The simulation results presented in section 5 demonstrate several interesting phenomena that may be taking place around Bennu based on the ejection events seen in early 2019. Figure 11. Map of the re-impact locations for the 11 February launch site case for the elevation direction sensitivity study, along with the associated latitude and longitude histograms.

Figure 12.
Re-impact locations for the 11 February launch site case for the elevation direction sensitivity study, with number of particles (indicated by colorbar) binned in 10 • by 10 • latitude-longitude bins.

Observed Outcomes of the Simulated Populations
The combination of dynamical processes acting on ejected particles can result in many particles not only surviving for multiple revolutions but also potentially surviving for more than one heliocentric orbit around asteroids. The grid of initial conditions explored here was fairly rough and by no means exhaustive. Thus, the fact that conditions that lead to orbits that survive for multiple months exist in all four ejection scenarios studied here implies that there is a non-negligible chance for long-lived orbits to occur in nature. Depending on how regularly such ejection events take place, and how many particles are released at these events, it is possible that some particles are in orbit around Bennu for significant periods of time. The ejection of particles and their subsequent motion also allows for mass movement at small near-Earth asteroids both across the surface, and leaving the system.
The range of particles studied here indicates that, over our grid space, a given particle has a greater than even chance of escaping the system. Those odds dramatically increase for smaller particles with high area-to-mass ratios. This implies that when particles are ejected from the surface, there is a deficit of smaller particles among those that return to the surface. If the ejection process also plays a role in creating small particles, there may be a general lack of sub-centimeter particles on the surface of Bennu. Similarly, if the ejection process is lofting particles that already exist on the surface, then over time, this process could clean the surface of free, small particles. Overall, the population of small surface particles will depend on the relative rates of their creation and subsequent removal through the ejection process.
These results also show that particles that return to re-impact the surface have significant mobility across the body. In all cases, re-impacting particles land preferentially at lower latitudes. A main reason for this is simply because Bennu has a larger radius near its equator. We do not consider here the dynamics of re-impact; however, it has already been established that the rotational Roche lobe for Bennu intersects the body around ± 20 • in latitude (Scheeres et al., 2019). Thus, particles that travel to this region are more likely to remain captured than those that re-impact at higher latitudes, which could further exacerbate the trend seen here. Importantly, this finding indicates that there could be a self-reinforcing mechanism at play: Once an equatorial bulge is established, ejected material is more likely to land there, thus increasing the radius of the bulge (and if material is coming from higher latitudes, decreasing the radius there), thereby exaggerating the "top" shape. Detailed simulation investigating how such a process might work in coupling the change in shape with the dynamics of ejected particles will be explored in future work.
This mass movement also provides a previously unconsidered mechanism which can contribute to crater erasure, especially at lower latitudes. Landslides are thought to be the main mechanism for crater erasure (Miyamoto et al., 2007), which should leave evidence of directional mass motion. Erasing craters through in-fall of ejected particles may not leave such prominent directional evidence, given that material can come from a variety of directions based on the variety of orbits and trajectories that can be established. However, considering the preferential loss of smaller particles through ejection, craters filled in this manner should preferentially contain larger particle sizes.

Dynamical Implications
For an ejected particle to survive in orbit for more than one revolution, there must first be a mechanism to raise the particle's periapse altitude before its first periapse passage. There are two ways to increase the periapse radius: either increase the semimajor axis (and thus the energy) or decrease the eccentricity.
The basis for understanding the rapid evolution of orbits around small bodies is given by Scheeres (2016), which accounts for the effects of the point mass gravity and SRP. That work shows that averaged over an orbit, SRP does not change the semimajor axis of an orbit, but it can change the eccentricity and the angular momentum in a coupled manner. Thus, SRP alone can increase survivability by lowering the eccentricity of some ejected particles. Furthermore, when a particle passes behind Bennu and is shadowed for some portion of its orbit, the SRP perturbation disappears. This changes the averaging results and can lead to a net gain in energy over these orbits.
However, Scheeres' theory cannot fully explain all of our simulated results. Our simulations show that Scheeres' theory describes the main evolution of particle orbits that are far from the surface (on the order of 1 km and above), for periods where the semimajor axis does not vary substantially. However, at lower altitudes, the non-spherical gravity and TRP provide significant perturbations that cause different evolution. TRP, in particular, can cause significant perturbations during low-altitude portions of the orbit, including at the initial stages of an orbit. The dominant component of the TRP acceleration is always in the radial direction away from the body, which can modify the eccentricity and, during some portions of an elliptical orbit, can lead to an energy change. Furthermore, because asteroids such as Bennu have a hot spot in the afternoon that is hottest at the equator, depending on the orientation of an orbit with respect to this hot region, there can be a net gain or loss in orbital energy as the particles fly past.
Beyond modifying the semimajor axis and eccentricity of the orbit, reorientation of the orbit plane and periapse location can also extend the orbital lifetime in two ways. First, if the location of periapse is moved to higher latitudes, the periapse altitude is increased because Bennu has a smaller radius at higher latitudes. Second, there can be a resonance between the precession of the orbit and the inertial precession of the thermal hot spot. The hot spot is always located at the same Bennu local time, but that location varies in inertial space as Bennu moves in its orbit about the Sun. If an orbit is oriented such that this hot spot adds energy through TRP, this relationship can be kept for many revolutions if the precession rates of the orbit line up appropriately. Orbital precession is caused by non-spherical gravity, third body gravity, and SRP (and to a lesser degree by TRP); thus, there is a complicated coupling between the various dynamical processes that can lead to a higher periapse and a longer orbit.
It is also pertinent to point out how the dynamics affect the escape speed of ejected particles. It has previously been noted that due to the significant spin rates and the complex shapes of small asteroids, the escape speed is not constant over the surface of the asteroid as is the case for a planetary body (Scheeres, 2016). Escape speeds are higher from potential lows on the surface, and particles can more easily achieve the escape speed if they are launched in the direction of surface motion (to the east typically), whereas they would have to be launched faster relative to the surface to achieve escape when launched in the direction opposite surface motion. However, SRP makes this even more complex and dependent on the area-to-mass ratio of the particles. Standard results from the literature indicate that SRP does not change orbital energy of unshadowed orbits, but this argument is based on treating SRP as a small perturbation and performing orbital averaging (Scheeres, 2016). In this scenario, these assumptions do not hold. Particles launched toward the Sun will lose energy and thus may not escape even though they are launched with a velocity above the local escape speed, and vice versa for those launched away from the Sun. Particles that do not escape will often subsequently approach close to the surface where other perturbations are significant enough to interfere with the averaging process. These effects become more severe as the particle area-to-mass ratio increases.
In short, a small asteroid ejecting particles is a rich and complex dynamical environment, and we have only explained some of the main mechanisms here. A detailed discussion and theoretical derivation to build upon current theories will be left to future work.

Limitations of the Presented Study
While our inferences are well supported by the simulations presented in this work, further investigation should be carried out to ensure these results are robust given the assumptions that have been made. Care should be taken in extrapolating these results for statistical interpretations because they are conditioned on a uniform grid across the input parameters. Furthermore, the population statistics presented here may be skewed by the range of parameters used, in particular with regard to particle size, which could exist at smaller sizes than we simulated. The simulations also only investigated particle dynamics associated with the three observed ejection events (four possible ejection sites) documented in , which share a late afternoon local time of launch and occurred relatively close to Bennu's perihelion. Finally, our simulated populations do not include very slow or very fast particles, which will clearly produce suborbital and direct escape trajectories, respectively. Therefore, in order to apply the results here in a statistical sense based on some distribution of launch conditions, the results must be weighted accordingly to account for particles outside the range used here.
Several other dynamical effects may be acting on these particles that are not included here. In particular, the particles could be shedding mass or outgassing after their release, creating an effective thrust, and possibly changing their area-to-mass ratio over time (Clark et al., 2004). Treating the particles as effective spheres for SRP and TRP modeling may also be inaccurate, and accommodations for the time-varying effects of a rotating flat plate may result in SRP acting in a slightly different direction, which would influence the results (Rosengren & Scheeres, 2014). Electrostatic forces are also not considered here but could be important near the surface (Hartzell & Scheeres, 2013;Hartzell, 2019), effectively modifying the launch conditions, what happens on low-altitude periapse passages, and the details of the landing locations. Finally, gas drag could play an important role at low altitudes, although the navigation team has determined it is insignificant at 1 km radius . Further investigation of these effects is warranted in the future.

Conclusion
We simulated the dynamical evolution of populations of particles similar to those that were ejected from Bennu in events observed by OSIRIS-REx in early 2019. We showed that the combined effects of gravity, solar radiation pressure, and thermal radiation pressure from Bennu can cause the orbits of many simulated particles to last for months or longer. Furthermore, the simulated populations exhibit two interesting phenomena that could play an important role in the geophysical evolution of bodies such as Bennu. First, small particles (<1 cm radius) are preferentially removed from the system, which could lead to a deficit of such particles on the surface. Second, re-impacting particles preferentially land near or on the equatorial bulge of Bennu. Over time, this can lead to crater in-filling and growth of the equatorial radius without requiring landslides.