# The Sands of Phobos: The Martian moon’s eccentric orbit refreshes its surface

###### keywords:

Solar System Science, Asteroids, Phobos, Regolith Science, Space Weathering^{†}

^{†}journal: Nature Geoscience

## Abstract

The surface of the Martian moon Phobos exhibits two distinct geologic units, known as the red and blue units. The provenance of these regions is uncertain yet crucial to understanding the origin of the Martian moon and its interaction with the space environment. Here we show that Phobos’ orbital eccentricity can cause sufficient grain motion to refresh its surface, suggesting that space weathering is the likely driver of the dichotomy on the moon’s surface. In particular, we predict that blue regions are made up of pristine endogenic material that can be uncovered in steep terrain subject to large variations in the tidal forcing from Mars. The predictions of our model are consistent with current spacecraft observations which show that blue units are found near these regions.

## 1 Introduction

Previous dynamical analyses of Phobos’ decaying orbit have shown that changes to its tidal environment lead to surface mobility and evolution, such as the formation of mass-wasting features Shi2016 () and linear grooves Hurford2016 (). These studies ignore the eccentricity of the moon’s orbit, as this would have a small effect over the dynamical timescales set by the orbital decay rate of the planetary satellite ( years). Here, we model the effect of Phobos’ eccentric orbit at its current distance, and demonstrate that time-varying forces can lead to substantial mass motion and surface erosion over much shorter timescales (- years). We find that librations induced by the eccentricity of the Martian moon can vary dynamical slopes by up to 2 per orbital period (7hr 39 min). These variations are sufficient to trigger a slow erosion process in high-slope regions that experience the largest changes in the tidal environment.

Using direct numerical simulations of particle dynamics, we demonstrate this new mechanism for surface mobility on Phobos for the first time. In addition, we find that these regions of high surface mobility coincide with one of two distinct ”color” units revealed by multispectral and hyperspectral imaging Fraeman2014 (): the blue unit. The blue unit is characterized by a relatively high albedo and increasing visible to infrared continuum slope. In contrast, the red unit is defined by a slightly lower albedo and relatively steeper spectral slope. It has been unclear whether these two geologic units are compositionally different, or if the spectral differences are due to space weathering, as there is an observed lack of strong absorption features in spectral observations of both unitsBasilevsky2014 (); Pieters2014 (); Murchie2015 (). Through a new eccentricity-driven mechanism for surface mobility, we develop a new model for regolith development on Phobos that can explain the relationship of the two color units. We show that this process can transform ”red” space-weathered Hapke2001 (); Vernazza2009 (); Kaluna2016 () regolith by exposing ”blue” sub-surface material in a process analogous to the tidally-induced refreshing of asteroid surfaces Binzel2010 ().

Our finding provides new perspectives on the space weathering process for airless bodies in the Solar System. For example, since the proposed mechanism makes predictions on the rate at which fresh material can be exposed, it is possible to place new constraints on the space weathering timescales of Phobos by comparing current and future observations of its surface. Additionally, a relatively fast surface-refreshing process on Phobos has implications on the origin of the Martian moons. An accretion scenario (following a giant-impact or during the formation of Mars) or a capture of an inner solar system body implies the presence of mafic materials on Phobos Fraeman2014 (); Thomas2011 (). Since our analysis suggests that the blue unit represents pristine endogenic material, the observed lack of mafic mineral absorptions in this unit Fraeman2012 () raises challenges for a giant impact scenario. Most of these challenges are likely to remain unsolved until dedicated missions to Phobos and Deimos will visit these remote objects. The Mars Moons eXploration missionKawakatsu2017 () (MMX), which will visit Phobos and return samples from its surface, will conclusively determine the origin of this enigmatic body and shed light on the history of the Martian system.

## 2 Results

We perform a two-step analysis to investigate the dynamical stability of grain particles on the surface of Phobos. In the first step, we model the gravitational attraction of the planetary satellite via a constant density polyhedron model WernerScheeres1996 () and calculate the evolution of its dynamical slopes (Fig. 1a) in the framework of the elliptical restricted three body problem (ERTBP, Broucke1969 ()). Despite the tiny eccentricity of Phobos’ orbit (), this model produces time-varying effects that change the local acceleration across the whole surface of Phobos. In particular, these fluctuations result in dynamical slope changes up to 2 per orbit that may trigger mass-wasting events near Stickney and the anti-Mars point (Fig. 1b). In the second step, we attempt to quantify the erosion rate in these regions.

In order to evaluate the magnitude of this effect on transforming Phobos’ surface features, we consider the analytical theory of erosion for transport-limited downslope flowCulling1960 (). For a planetary surface with loose regolith (with little to no cohesion) where the flow of grains is controlled by the transportation rate rather than the regolith supply or production rate, the flow rate, , can be defined as a function of the local surface slope, , as

(1) |

where is the downslope flow constant in units of volume flux per unit time. In previous studies of regolith mobility on small bodiesRichardson2014 (), the value of was estimated for impact-induced disturbances using a Newmark slide-block modelRichardson2005 (), and the value of the flow rate was shown to be linear with the slope only when circa Richardson2014 (). For higher slopes, fast granular flows are better modeled as a non-linear function that takes the critical slope, , into consideration Roering1999 (). Here, we perform local simulations of grain dynamics in order to estimate the value of for a variety of surface slopes and periodic fluctuations, i.e., . Using a soft-sphere discrete element code, PKDGRAVRichardson2011 (); Schwartz2012 (), we set up simulations of a 1.5 m by 1.2 m regolith bed with a porosity of 45% settled in Phobos gravity ( mm/s) and made up of grains with radii between 1.4 and 1.6 cm. The material properties of the grains are similar to sand of medium hardnessJiang2015 (). The simulated grains have a measured critical angle of repose of . In order to allow the grains to flow freely in the event of an avalanche trigger or creep-induced motion, we impose periodic boundaries in the directions perpendicular to the initial local gravity vector. The simulation setup is illustrated in Fig. 2a.

With the goal of determining a conservative estimate for expected volume flux across different regions of Phobos, we perform simulations for initial surface slopes between 5 and 32 with . The highest surface slope achieved across all of the simulations is 34.5. The output of the simulations display creeping motion of small grain displacements over the course of a Phobos orbital period. In particular, we find that a granular bed that initially starts at a high slope with significant slope changes will experience the most mass flow. However, even granular beds that are in a sub-critical state and are subject to relatively small variations can lead to surface mobility. We term this a cold flow process as opposed to the “fast” flow of an actual landslide.

The cold flow process is illustrated in Fig. 2b with an example of local simulations for the case of a regolith bed that has a mean slope of 31 deg and between and . Here, we highlight the largest particle displacements (anything greater than 3 cm) by showing the initial (black) and final (yellow) positions of the particles along with their tracks (dotted red lines). For each simulation, we measure the total displacement of particles in order to calculate the volume flux across the simulation area over one orbit. For the phase space studied here, which consists of regolith resting at sub-critical angles, we calculate volume flux rates between 10 to 4 particles/m/orbit. We show the results of all simulations in Fig. 2c.

We find that across low and high slopes, the expected volume flux due to eccentricity-driven surface slope changes, , is best described by,

(2) |

Based on the simulations conducted in this study, we find that the magnitude of the particle flow can be described by two linear regimes with different values of the flow constant . The cold flow prescription depends on the proximity of the initial slope to the critical angle of repose and the magnitude of the surface slope-change. We find that the value of switches from particles/m/orbit to particles/m/orbit when , which corresponds to a value of when . Therefore, equation (2) provides a description for the cold-flow mechanism on the surface of Phobos induced by its eccentric orbit. We use this prescription to compare the expected erosion rates in different regions of Phobos with observational data and identify areas where the cold flow mechanism may be dominating the geomorphology of the Martian moon.

Observations by the HiRISE camera on NASA’s Mars Reconnaissance Orbiter (MRO) spacecraft and the High Resolution Stereo Camera (HRSC) on ESA’s Mars Express spacecraft have produced detailed images showing the heterogeneity of the area in and around Stickney crater (Fig. 3a,d). While high albedo features are typically indicative of mass-wasting motion, the area east of Stickney has a relatively smooth surfaceBasilevsky2014 (), free of streaks that typify granular flow on Phobos Shi2016 (). There have been some suggestions that the color variation in this area is a result of a catastrophic landslide that began in Stickneyâs western wall and transitioned to the sub-Martian point by motion of a thin surface layer of material in a long run-out landslide ShingarevaKuzmin2001 (); CollinsMelosh2003 (). However, these studies either assume small and unrealistic values for the friction properties of the Phobos surface (friction coefficient of 0.01-0.09ShingarevaKuzmin2001 ()), or require the development of acoustic fluidization CollinsMelosh2003 (). Another interpretation of this region is that it is asymmetric ejecta deposit from StickneyThomas1998 (). Here, we offer a novel interpretation.

By analyzing the spatial distribution of the surface slopes in the region east of Stickney (Fig. 3b) and the predicted variation due to Phobos’ eccentric orbit, we predict a relatively large amount of regolith flow compared to nearby regions (Fig. 3c). This area, labeled A in Fig. 3a, should have the most mobile regolith due to the unique combination of steep terrain and high surface-slope changes. In contrast, the relatively steep southern wall of Stickney (labeled B) does not experience large variations in the tidal environment of Mars and–in agreement with our model–shows little indication of blue material (Fig. 3d).

Through this qualitative comparison of Phobos color data with our predicted volume flux, we propose that blue geologic units are composed of sub-surface material that is being continually exposed by the eccentricity-driven librations of Phobos. This suggests that the blue units are the footprints of “fresh”, pristine material, as opposed to the “old”, space weathered regolith that is found in the red unit of the Martian moon. To evaluate this possibility, we compare our model’s predicted excavation rate to the likely space weathering rate at Phobos. The results of this analysis are discussed in the next section.

## 3 Discussion

In order to determine the efficiency of this process in uncovering fresh material, we consider the possible timescale and depth of space weathering on Phobos. Since we have a poor understanding of Phobos’ composition Murchie2015 () and the specific mechanisms that may influence the maturation of its regolith PietersNoble2016 (), we consider a range of space weathering rates to provide some perspective on the possible surface refreshing rate. We use estimates for space-weathering based on previous spectral analysis of main-belt asteroids. Studies on space weathering of S-complex asteroid dynamical families Vernazza2009 (); Binzel2010 () reveal that space-weathering rate for these bodies is relatively rapid ( years), suggesting that solar wind ion implantation dominates the process. In contrast, analysis of C-complex asteroids Kaluna2016 () suggests that space weathering is a long term process, with spectral slopes being modified on timescales closer to years. Using the proposed space-weathering rates Kaluna2016 () and slope changes in CRISM reflectance measurements between blue and red geologic units Fraeman2012 (), we estimate that space weathering on Phobos can occur on timescales with a lower limit of 10 years.

Due to the lack of in-situ exploration of small-body subsurface properties, we rely on analysis of lunar core samples from the Apollo missions to obtain an estimate on the depth of space-weathered material on Phobos. Measurements of the regolith maturity parameter ( ratio) of 12 Apollo core samples show that the regolith typically transitions between mature weathered material to relatively fresh material around 20-50 cm below the surface Lucey2006 (). Combining this estimate of the depth of weathered material with our assumptions on the space-weathering rate and the predicted excavation rates derived from our numerical simulations, we analyze the efficiency of the eccentricity-driven cold-flow process in exposing fresh material. We find that, while the calculated excavation rates are relatively small (only a few grains typically move from a 1 m high-slope region every Phobos orbit), the accumulated effect over the course of Phobos’ geologic time can be quite substantial. In Fig. 4, we show the excavation rates for the cold-flow of mm-size particles required to expose material at certain depths (represented by the solid curves) for different timescales. It is important to note that we make a steady-state assumption for regolith as it becomes exposed. In reality, the regolith can mature as it is uncovered. Nevertheless, this analysis provides an order of magnitude estimate of the necessary excavation rates required to expose fresh material.

From our numerical simulations of the cold-flow mechanism, we show that the simulation-derived mean excavation rates (dashed lines) for the labeled regions in Fig. 3a. The values of and for these regions are determined based on our dynamical calculation (Fig. 1), and are fed into the prescription of mass flux, equation (2). We find that the regions with the highest spatial distribution of blue material in the HRSC data (region A) has the highest mean excavation rate, followed by the western wall of Stickney (region D), and both the southern wall (region B) and the crater floor (region C) have very low predicted excavation rates. Overall, we find that blue regions represent areas where grain-mobility is high enough to expose fresh material down to depths of 20 cm at a rate faster than the expected space weathering time-scale (grayed-out region). Red regions have predicted mean excavation rates that may only expose material a few cm deep before years, which is the age where our steady-state assumption may fail, as some asteroidsBinzel2010 () are found to rapidly redden in this time-scale. Therefore, we conclude that eccentricity-driven mass-wasting is an efficient process for uncovering pristine un-weathered material on Phobos, represented by blue geologic units found near Phobosâ sub-Mars point. We propose a model for regolith development where i) Phobos is inherently blue, ii) space weathering reddens the surface material, iii) impacts and Phobos’ decaying orbit can expose fresh or weathered material depending on crater depth and the characteristic space weathering depth , iv) this freshly exposed material is reddened over time, and v) blue regions are continuously refreshed by mass-wasting that occurs in high-slope regions that experience significant variations in tidal forces due to Phobosâ eccentric orbit.

## 4 Perspectives

We have shown that Phobos’ small eccentricity can drive changes in the surface slope that mobilize grains on the top-layers of high-slope regions. While the rate of mass flow is slow, the accumulated effect over time is enough to expose fresh subsurface material, giving a plausible explanation for the dichotomy between red and blue geologic units on the surface. The analysis described here relies on a simplistic model of space weathering based on other airless bodies, and does not take into account Phobos’ unique environment, which is not well understood. However, since the predicted surface refreshing is tied to Phobos’ orbital motion, we can place new constraints on the space weathering timescale for Phobos itself, which may then allow us to better understand this process for other airless bodies in our Solar System. With a better understanding of the surface and subsurface conditions on the Martian moon, the Phobos surface can become a “Rosetta-stone” for understanding space weathering elsewhere in the solar system. MMXKawakatsu2017 () (led by JAXA, with contributing instruments from CNES and NASA) will explore Phobos and obtain a core sample from its surface at a depth of up to 10 cm. The combination of a returned sample, new remote sensing observations of the Phobos surface, and a better understanding of surface-refreshing process on Phobos will revolutionize our understanding of the interaction between the surface of an airless body and the space environment.

Finally, one of MMX’s main goals is to provide an answer to the question of Phobos’ enigmatic origin. It is unclear whether Phobos is a captured Solar System body Burns1978 (); Hansen2018 (), or if it formed in an accretion disk following a giant impact with Mars Craddock1994 (); Hyodo2017 (). Current theories of a giant impact or in-situ formation origin predict that Phobos should be made up of a relatively large proportion of Martian material and imply the presence of mafic minerals such as olivine and pyroxene. The observed lack of diagnostic absorption features in spectral data suggests that either: a) all of Phobos’ surface is heavily space weathered, removing spectral features, b) the surface grains are extremely fine, or c) the regolith is substantially mixed with opaque materials such as carbon. Our analysis suggests that the blue units on Phobos represent pristine endogenic material. Therefore, the lack of mafic absorption features in the spectral data of blue units Fraeman2012 (), which are relatively featureless, would suggest that the first scenario (a) is unlikely, and that it is possible that Phobos has a sparse abundance of Martian material. Therefore, we find that this cold-flow mechanism which can refresh the Phobos surface indicates that Phobos is not spectrally like Mars; or, the observed lack of mafic absorptions must be explained in some other scenario, such as those listed above.

## Methods

### The dynamical environment on Phobos’ surface

A point mass in the vicinity of Phobos is subject to the gravitational attraction of both the planetary satellite itself and Mars. Accordingly, the acceleration felt by particle grains on the surface of the Martian moon can be approximated by the equations of the ERTBP Broucke1969 (). Since the distance between the grain and the barycenter of Phobos is much smaller than the distance between the grain and the barycenter of Mars, the equations of motion can be further simplified by expanding the gravitational term of Mars in a Taylor series up to the first order. The resulting system of equations of motion is known as the Hill approximation of the ERTBP and is better understood in a pulsating rotating reference frame centered on Phobos and such that is constantly align with an imaginary line connecting the two primaries, is parallel to the orbital angular momentum of the moon, and . Following the length and time scaling described in Scheeres2002 (), the equations read as

(3) |

where denotes differentiation with respect to , i.e., the true anomaly of Phobos, represents the state of the mass particle, is the eccentricity of Phobos, and is the gravitational attraction of the Martian moon. The acceleration is obtained in closed form via a constant density polyhedron model consisting of vertices, facets, and density equal to km/mWerner1996 (); Willner2014 (); Murchie2015 ().

In this paper, we are interested in monitoring the variation of the surface slopes over one orbital period of Phobos around Mars. Accordingly, let and denote the normal vector and barycenter position vector of the i-th facet, respectively. The surface slope is obtained as the supplement of the angle between and the acceleration vector evaluated in at , i.e.,

(4) |

While Phobos rotates at a constant rate of rev/day, the pulsating reference frame rotates about the -axis with a time-changing angular velocity of . Consequently, as seen from the rotating reference frame, the polyhedral shape of Phobos describes librations of approximately deg about the axis. These librations cause periodic changes in the tidal acceleration due to Mars that may be resposible for triggering the cold flow mechanism in high-surface slopes regions of Phobos.

### Granular Dynamics Simulation

For the local simulations of granular flow, we used the N-body code PKDGRAV Richardson2011 () which is capable of accurately simulating the complexity of grain-grain and grain-boundary interactions through a soft-sphere discrete element method Schwartz2012 () (SSDEM). In SSDEM, collisions of spherical grains are resolved by allowing them to slightly overlap (typically % of their radii) and then applying multi-contact and multi-frictional forces, including static, rolling and twisting friction. Modeling grain friction accurately is a critical component for high-fidelity granular physics simulations. For PKDGRAV, a new rolling friction model has recently been implemented Zhang2017 () that allows for a more accurate modeling of grain shape and angularity, through a shape parameter .

The mechanical properties of Phobos’ surface are poorly known; however, the range of plausible material types are limited, and some remote sensing observations can help constrain some important parameters that are required for accurate simulations. Radar albedo measurements Busch2007 () provide an estimate of the near surface porosity (%). Local slopes on the surface of Phobos can be as high as 45; however, these are likely regions that are in a super-critical state and may not accurately reflect the material’s true angle of repose. We chose a set of friction parameters for the grains such that they emulate the static behavior of sand of medium hardness as determined by laboratory experimentsJiang2015 (), which have a corresponding angle of friction of 34. The grain size of Phobos regolith is highly uncertain. Early measurements from the Viking spacecraft suggested that the typical grain size of Phobos regolith is 50-100 mSasaki1995 (). More recent work that used thermal inertia measurements combined with a thermal conductivity model for regolith suggests a larger average grain size of 1.1 mmGundlachBlum2013 () ; however, there are large uncertainties ( mm). We performed a few high-resolution simulations with grains that have radii of = 0.9-1.1 mm; however, the combination of particle number (50,000 particles), simulation time (7 h 39.2 m), and small timestep (0.05 s) make the computation very expensive. Therefore, to speed up the calculation, we used a few high-resolution simulations as a baseline, and the majority of our simulations were lower-resolution cases ( = 1.4-1.6 cm, N=7,000). This was necessary to sweep a larger parameter space in starting slope, , and slope change amplitude, . Table 1 summarizes the properties of grains and regolith bed for our simulations.

Property | Low Resolution | High Resolution |
---|---|---|

1.4-1.6 cm | 0.09-0.11 cm | |

34 | 34 | |

3.1 g cm | 3.1 g cm | |

1.7 g cm | 1.7 g cm | |

7,000 | 50,000 | |

Sim. Area | 1.5 1.2 m | 0.60 0.4 m |

Since we are interested in measuring the possibility of surface grain motion due to these small variations in the surface slope, and not complete hillslope-failure, we restrict the vertical scale of our simulations to 5 grain diameters. We set up simulations with initial slopes ranging from 5 to 32 and amplitudes of variation between 0.5 and 1.5. This is done by performing an initial setup simulation where the direction of the gravity vector is changed till the desired gravity slope is achieved. The change in the gravity vector is done gradually (over 1 hr in simulated time) as a precaution so as to prevent granular flow from triggering before the time-varying slope change is simulated. The grains are then allowed to settle till all small motions are dissipated. Only after this setup is complete do we perform a full simulation, varying the gravity slope of the inclined granular bed over 1 Phobos orbit (7hr 39 minutes).

For the limited number of high-resolution cases, we found that, for the same and , the particle flux (particle/m/orbit) was similar to low resolution cases. Therefore, we expect the volume flux (m/m/orbit) to depend on the actual grain size on the surface of Phobos. For the same period of time, a surface made up of small particles will excavate material at a given depth at a slower rate than a surface with larger particles. Future work will focus on testing this finding in order to develop scaling laws for the volume flux as a function of grain size.

While the grains we simulate here are cohesion-less, the actual grains on Phobos are likely to have some non-zero cohesion that will influence the grain dynamics Scheeres2010 (). Cohesive bonding between grains may act as an inhibitor to particle motion, increasing an individual particle’s shear strength; however, it is not obvious how grains may actually bond in the Phobos environment. In particular, a size distribution of grains may lead to small grains being bonded to larger ones, and being carried away by the eccentricity-driven mass-wasting process, enhancing the expected excavation rates. Furthermore, as the sizes of the grains become smaller, the relative importance of cohesion is magnified, and processes such as grain loftingHartzellScheeres2011 (), driven by charging from the solar wind or from passing coronal mass ejectionsFarrel2017 (), may be present on Phobos’ surface. Therefore, the contribution of cohesion to this process may be non-trivial, and we will investigate its influence in the future with a recently implemented cohesion-model in PKDGRAVZhang2018 ().

### Comparing excavation rates to space weathering rates

We calculate the necessary volume flux required to uncover fresh material (Fig. 4) by equating an excavation rate with a space weathering rate to a certain depth. We define the excavation rate, , to be the volume, , out of an area, , for a given time period, ,

(5) |

and we define the space weathering rate, , as the timescale, necessary to mature regolith up to a depth of .

(6) |

We calculate the equilibrium values for excavation and weathering to different depths by equating equations (5) and (6), setting the area to units of 1m, to one Phobos orbit, and the grain radius, to the volume of a mm-size particle, such that the number of particles that must flow out of a 1m region in 1 Phobos orbital period, , to excavate that region to a depth of in years is given by:

(7) |

### Space weathering timescales for Phobos from spectral slopes

We derive spectral slopes from CRISM spectra of representative blue and red regions on Phobos corrected to I/F with incidence and phase angles equal to 30 degrees and emission angle of 0 degrees Fraeman2012 (). We measure spectral slopes following Kaluna2016 () using the normalized reflectivity gradient S’ defined by JewittMeech1986 ():

(8) |

where is the slope of the reflectance between wavelengths m and m and slopes are normalized to reflectance at m. We find m and m, for an overall slope change from blue to red of m. Slope change rates due to space weathering have been estimated for several asteroid families of different compositions: Kaluna2016 () observed a slope change of m over 2.3 Gyr for C-complex asteroids, while Lazzarin2006 () determined a rate of AUmMyr for C-complex and AUmMyr for S-complex asteroids. Vernazza2009 () determined that asteroids may weather relatively quickly, with a slope change of m on the order of a million years. Although the space weathering environment may be very different at Phobos than for asteroids, and Phobosâ composition is not well constrained, these space weathering rates can be used to obtain an order of magnitude of the effect expected for Phobos. Using these slope change rates due to space weathering of asteroid families and the slope difference of m between blue and red regions, we estimate a space weathering age for Phobosâ red regions on the order of yr. While the upper time range based on slow C-type asteroid reddening is physically infeasible, it is the lower time limit that is relevant to our study, in order to determine whether the cold flow mechanism is efficient enough to compete with rapid space weathering.

## References

## References

- (1) Shi, X., Oberst, J., Wilner, K. Mass wasting on Phobos triggered by an evolving tidal environment. Geophysical Research Letters. 43, 12, 371-12, 379 (2016).
- (2) Hurford, T.A., et al. Tidal disruption of Phobos as the cause of surface fractures. J. Geophys. Res. Planets 121, 1054-1065 (2016).
- (3) Fraeman, A.A., Murchie, S.L., et al. Spectral absorptions on Phobos and Deimos in the visible/near infrared wavelengths and their compositional constraints. Icarus 229, 196-205 (2014).
- (4) Basilevsky, A.T., Lorenz, C.A., et al. The surface geology and geomorphology of Phobos. Planetary & Space Science 102, 95-118 (2014).
- (5) Pieters, C.M., Murchie, S., Thomas, N., Britt, D. Composition of surface materials on the Moons of Mars. Planetary and Space Science 102, 144-151 (2014).
- (6) Murchie, S.L., Thomas, P.C., Rivkin, A.S., Chabot, N.L. in Asteroids IV (eds. Michel, P., DeMeo, F.E., Bottke, W.F.) 451-467 (University of Arizona Press, 2015).
- (7) Hapke, B., Space weathering from Mercury to the asteroid belt. Journal Geophysical Research 106, 10,039-10,073 (2001).
- (8) Vernazza, P., Binzel, R.P., Rossie, A., Fulchignoni, M., Birlan, M. Solar wind as the origin of rapid weathering of asteroid surfaces. Nature 458, 993-995 (2009).
- (9) Kaluna, H.M., Masiero, J.R., Meech, K.J. Space weathering trends among carbonaceous asteroids. Icarus 264, 62-71 (2016).
- (10) Binzel, R.P., et al. Earth encounters as the origin of fresh surfaces on near-Earth asteroids. Nature 463, 331-334 (2010).
- (11) Thomas, N., et al. Spectral heterogeneity on Phobos and Deimos: HiRISE observations and comparisons to Mars Pathfinder results. Planetary and Space Science 59, 1281â1292 (2011).
- (12) Fraeman, A.A., et al. Analysis of disk-resolved OMEGA and CRISM spectral observations of Phobos and Deimos. Journal Geophysical Research 117, E00J15 (2012).
- (13) Kawakatsu, Y., et al. M. Mission concept of Martian Moons eXploration (MMX). 68th International Astronautical Congress, IAC-17-A3.3A.5 (2017).
- (14) Werner, R.A., & Scheeres, D.J. Exterior gravitation of a polyhedron derived and compared with harmonic and mascon gravitation representations of asteroid 4769 Castalia. Celestial Mechanics and Dynamical Astronomy 65, 313-344 (1996).
- (15) Broucke, R. Stability of periodic orbits in the elliptic, restricted three-body problem. AIAA Journal 7, 6, 1003-1009 (1969).
- (16) Culling, E.H. Analytical theory of erosion. The Journal of Geology 68, 336â344 (1960).
- (17) Richardson, J.E., & Bowling, T.J. Investigating the combined effects of shape, density, and rotation on small body surface slopes and erosion rates. Icarus 234, 53-65 (2014).
- (18) Richardson, J.E., Melosh, H.J., Greenberg, R.J., OâBrien, D.P.. The global effects of impact-induced seismic activity on fractured asteroid surface morphology. Icarus 179, 325â349 (2005).
- (19) Roering, J.J., Kirchner, J.W., Dietrich, W.E. Evidence for nonlinear, diffusive sediment transport on hillslopes and implications for landscape morphology. Water Resource Research 34, 853-870 (1999).
- (20) Richardson, D.C., Walsh, K.J., Murdoch, N., Michel, P. Numerical simulations of granular dynamics: I. Hard-sphere discrete element method and tests. Icarus 212, 427-437 (2011).
- (21) Schwartz, S.R., Richardson, D.C., Michel, P. An implementation of the soft-sphere discrete element method in a high-performance parallel gravity tree-code. Granular Matter 14, 363-380 (2012).
- (22) Jiang, M. , Shen, Z. , Wang, J. A novel three-dimensional contact model for granulates incorporating rolling and twisting resistances. Computational Geotechnology 65, 147â163 (2015).
- (23) Patsyn, V., et al. Spectrometric characteristics of the surface of Phobos from data obtained by HRSC on Mars Express. European Planetary Science Congress, Madrid (2012).
- (24) Shingareva, T.V., & Kuzmin, R.O. Mass-wasting processes on the surface of Phobos. Solar System Research 35, 431-443 (2001).
- (25) Collins, G.S., & Melosh, H.J. Acoustic fluidization and the extraordinary mobility of sturzstroms. Journal of Geophysical Research 108, 2473 (2003).
- (26) Thomas, P.C. Ejecta emplacement on the Martian satellites. Icarus 131, 78â106 (1998).
- (27) Pieters, C.M., Noble, S.K. Space weathering on airless bodies. Journal of Geophysical Research 121, 1865-1884.
- (28) Lucey, P., et al. in New Views of the Moon (eds. Jolliff, B.L., Wieczorek, M.A., Shearer, C.K., Neal, C.R.) 83-219 (Mineralogical Society of America, Washington D.C., 2006).
- (29) Burns, J.A. The dynamical evolution and origin of the Martian moons. Vistas in Astronomy 22, 193-208 (1978).
- (30) Hansen, B.M.S. A dynamical context for the origin of Phobos and Deimos. Monthly Notices of the Royal Astronomical Society 475, 2452-2466 (2018).
- (31) Craddock, R.A. The origin of Phobos and Deimos. Lunar and Planetary Science Conference 25, 293 (1994).
- (32) Hyodo, R., Genda, H., Charnoz, S., Rosenblatt, P. On the impact origin of Phobos and Deimos. I. thermodynamic and physical aspects. The Astrophysical Journal 845, 125-133 (2017).
- (33) Scheeres, D. J., Marzari, F. Spacecraft dynamics in the vicinity of a comet. Journal of the Astronautical Sciences 50, 1, 35-52 (2002).
- (34) Werner, R. A., Scheeres, D. J. Exterior gravitation of a polyhedron derived and compared with harmonic and mascon gravitation representations of asteroid 4769 Castalia, Celestial Mechanics and Dynamical Astronomy 65, 3, 313-344 (1996).
- (35) Willner, K., Shi, X., Oberst, J. Phobos’ shape and topography models. Planetary and Space Science 102, 51-59 (2014).
- (36) Zhang, Y., et al. Creep stability of the proposed AIDA mission target 65803 Didymos: I. Discrete cohesionless granular physics model. Icarus 294, 98-123 (2017).
- (37) Busch, M.W., et al. Arecibo radar observations of Phobos and Deimos. Icarus 186, 581-584 (2007).
- (38) Sasaki, S. Surface properties of Phobos/Deimos and Formation of self-sustained Martian dust torus. Lunar & Planetary Science Conference 26, 1219 (1995).
- (39) Gundlach, B., & Blum, J. A new method to determine the grain size of planetary regolith. Icarus 223, 479-492.
- (40) Scheeres, D.J., Hartzell, C.M., Sánchez, P. Swift, M. Scaling forces to asteroid surfaces: The role of cohesion. Icarus 210, 968-984 (2010).
- (41) Hartzell, C., & Scheeres, D.J. Dynamics of levitating dust particles near asteroids and the moon. Journal of Geophysical Research 118, 116-125 (2013).
- (42) Farrel, W.M., et al. Anticipated electrical environment at Phobos: nominal and solar storm conditions. Advances in Space Research, In Press (2017).
- (43) Zhang, Y., et al. Rotational Failure of rubble-pile bodies: Influences of shear and cohesive strengths. The Astrophysical Journal 857, 15-35 (2018).
- (44) Jewitt, D., Meech, K.J. Cometary grain scattering versus wavelength, or ‘What color is comet dust’? The Astrophysical Journal 310, 937â952 (1986).
- (45) Lazzarin, M., et al. Space weathering in the main asteroid belt: The big picture. The Astrophysical Journal 647, L179-L182 (2006).

## Acknowledgements

R.L.B acknowledges support from JAXA’s Aerospace Project Research Associate Program. N.B. conducted this work as a JSPS International Research Fellow. S.T.C. was supported by the JAXA International Top Young Fellowship Program. The authors also thank Patrick Michel for constructive feedback on the results and implications of this work. Grain dynamics simulations were calculated on YORP cluster run by the Center for Theory and Computation at the Department of Astronomy at the University of Maryland. For data visualization, the authors made use of the freeware, multi-platform, ray-tracing package, Persistence of Vision Raytracer.

## Author Contributions

R.L.B. conceptualized the study, designed and performed the local simulations of granular dynamics, and led the research. N.B. initiated the project through a study of the three body elliptical problem on Phobos, performed the gravitational dynamics calculations, and contributed to the analyses. S.T.C. provided geophysical and geomorphological expertise, remote sensing analysis, and constructed the model for regolith development on Phobos. Y.K. provided guidance and advice on the formulation and scope of the research. M.F. provided guidance and discourse on the implications of the results. Y.K. and M.F. provided expertise in small body exploration and contextualized the research in the frame of the MMX mission. All authors contributed to the interpretation of the results and preparation of the manuscript.

## Competing interests

The authors declare no competing interests.

## Data Availability

The datasets generated and analyzed during the current study are available from the corresponding author on reasonable request.

## Code Availability

The code used to generate the datasets are available from the corresponding author on reasonable request.