ABSTRACT
Title of dissertation: NUMERICAL SIMULATIONS
OF GRANULAR PHYSICS
IN THE SOLAR SYSTEM
Ronald-Louis G. Ballouz, Doctor of Philosophy, 2016
Dissertation directed by: Professor Derek C. Richardson
Department of Astronomy
Granular physics is a sub-discipline of physics that attempts to combine prin-
ciples that have been developed for both solid-state physics and engineering (such
as soil mechanics) with 
uid dynamics in order to formulate a coherent theory
for the description of granular materials, which are found in both terrestrial (e.g.,
earthquakes, landslides, and pharmaceuticals) and extra-terrestrial settings (e.g.,
asteroids surfaces, asteroid interiors, and planetary ring systems). In the case of our
solar system, the growth of this sub-discipline has been key in helping to interpret
the formation, structure, and evolution of both asteroids and planetary rings. It is
dicult to develop a deterministic theory for granular materials due to the fact that
granular systems are composed of a large number of elements that interact through
a non-linear combination of various forces (mechanical, gravitational, and electro-
static, for example) leading to a high degree of stochasticity. Hence, we study these
environments using an N -body code, pkdgrav, that is able to simulate the gravita-
tional, collisional, and cohesive interactions of grains. Using pkdgrav, I have studied
the size segregation on asteroid surfaces due to seismic shaking (the Brazil-nut ef-
fect), the interaction of the OSIRIS-REx asteroid sample-return mission sampling
head, TAGSAM, with the surface of the asteroid Bennu, the collisional disruptions
of rubble-pile asteroids, and the formation of structure in Saturn's rings. In all of
these scenarios, I have found that the evolution of a granular system depends sensi-
tively on the intrinsic properties of the individual grains (size, shape, sand surface
roughness). For example, through our simulations, we have been able to determine
relationships between regolith properties and the amount of surface penetration a
spacecraft achieves upon landing. Furthermore, we have demonstrated that this
relationship also depends on the strength of the local gravity. By comparing our
numerical results to laboratory experiments and observations by spacecraft we can
begin to understand which microscopic properties (i.e., grain properties) control the
macroscopic properties of the system. For example, we can compare the mechan-
ical response of a spacecraft to landing or Cassini observations of Saturn's ring to
understand how the penetration depth of a spacecraft or the complex optical depth
structure of a ring system depends on the size and surface properties of the grains
in those systems.
NUMERICAL SIMULATIONS OF
GRANULAR PHYSICS IN THE SOLAR SYSTEM
by
Ronald Ballouz
Dissertation submitted to the Faculty of the Graduate School of the
University of Maryland, College Park in partial fulllment
of the requirements for the degree of
Doctor of Philosophy
2016
Advisory Committee:
Professor Derek Richardson, Chair/Advisor
Professor Doug Hamilton
Professor Jessica Sunshine
Dr. Thomas Statler
Professor Wolfgang Losert
Dr. Olivier Barnouin
c
 Copyright by
Ronald Ballouz
2016
To my wife Kellie, my parents Eppie & Georges, my brothers Samir, Tony, and
David, and to Maddie & Mickey.
ii
Acknowledgments
I owe so much gratitude to my advisor, Derek Richardson, for helping me
throughout the last 5 years that we have been working together. I came into the
Department of Astronomy at the University of Maryland thinking that I would
pursue a career in galactic astronomy. When I came for a prospectives visit to
the University of Maryland, one of the rst people I talked to was Derek, and his
kindness and intelligence left an impression on me. After he gave the rst ever
colloquium I attended at the department, I knew that I was going to switch to
planetary sciences. Derek has been an engaging and supportive advisor. Above all
else, I'd like to thank him for his patience.
I would also like to thank Patrick Michel for his hospitality during my months
in Nice, and whose valuable insights into the eld have shaped much of my thinking.
I also owe him and Derek so much thanks for setting an excellent example, and
developing a wonderful team of seasoned and burgeoning scientists who somehow
all got interested in the mechanics of tiny grains on far-away worlds. Much of
this thesis has come about from discussions with team members that I've had the
luck and pleasure to coincide with: Kevin Walsh, Steve Schwartz, Yang Yu, Clara
Maurel, and Yun Zhang.
I am grateful for my fellow graduate students, who have helped me survive and
appreciate graduate school. Thank you Johnny, Vicki, Krista, Gabriele, Mahmuda,
and Maggie, for sticking around, sharing a laugh, and encouraging me when I felt
hopeless.
iii
I also have to thank the members of my thesis committee who have stuck with
me since the very beginning of my PhD work. Thank you Jess for reminding me
to keep the science in focus. Thank you Doug for showing me how fun planetary
dynamics can be. Thank you Tom for teaching me the importance of being cognizant
of the scientic world outside one's own little bubble.
Finally, I'd like to thank my family. My parents for trusting me to make my
own choices. My brothers for their encouragement. My wife for her constant love
and patience, and our dogs for their unbridled enthusiasm.
iv
Table of Contents
List of Tables viii
List of Figures ix
List of Abbreviations xi
1 Introduction 1
1.1 Evidence of Regolith and its Mobility on Asteroid Surfaces . . . . . . 5
1.2 Numerical Simulations of Granular Physics . . . . . . . . . . . . . . . 12
1.3 pkdgrav . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 12
1.3.1 Computational Resources . . . . . . . . . . . . . . . . . . . . 16
1.3.2 Code Improvements: Fixed Search Ball . . . . . . . . . . . . . 17
1.4 Thesis Outline: Granular Physics in the Solar System . . . . . . . . . 21
1.4.1 Regolith Mobility on Asteroids . . . . . . . . . . . . . . . . . 21
1.4.2 Spacecraft-Regolith Interactions . . . . . . . . . . . . . . . . . 22
1.4.3 Catastrophic Disruption of Rubble-Pile Asteroids . . . . . . . 23
1.4.4 Structure Formation in Saturn's B Ring . . . . . . . . . . . . 24
1.4.5 Outlook and Future Work . . . . . . . . . . . . . . . . . . . . 27
2 Size Segregation on Asteroid Surfaces: Brazil-nut Eect 28
2.1 Chapter Preface . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 28
2.2 Background . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 30
2.3 Methodology . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 34
2.3.1 pkdgrav: PBC . . . . . . . . . . . . . . . . . . . . . . . . . . 34
2.3.2 Friction parameters . . . . . . . . . . . . . . . . . . . . . . . . 36
2.3.3 Simulation parameters and walls setup . . . . . . . . . . . . . 39
2.3.4 Initialization process and main simulations . . . . . . . . . . . 42
2.4 Results . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 43
2.4.1 Earth-gravity cases: Walls versus PBC . . . . . . . . . . . . . 43
2.4.2 Earth-gravity cases: static and rolling friction . . . . . . . . . 50
2.4.3 Low-Gravity Cases . . . . . . . . . . . . . . . . . . . . . . . . 55
2.5 Discussion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 66
v
2.5.1 Driving mechanisms . . . . . . . . . . . . . . . . . . . . . . . 66
2.5.2 Interplay between gravity and rolling friction . . . . . . . . . . 73
2.6 Conclusion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 76
3 Low-speed granular impacts:Modeling Spacecraft-Regolith Interactions 80
3.1 Chapter Preface . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 80
3.2 Background . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 81
3.3 Methodology . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 85
3.3.1 Inertial walls . . . . . . . . . . . . . . . . . . . . . . . . . . . 86
3.3.1.1 TAGSAM specications . . . . . . . . . . . . . . . . 88
3.4 Impact Experiments . . . . . . . . . . . . . . . . . . . . . . . . . . . 92
3.4.1 Material Property . . . . . . . . . . . . . . . . . . . . . . . . . 93
3.4.2 Penetration Depth Experiments . . . . . . . . . . . . . . . . . 97
3.4.3 Gravity Dependence . . . . . . . . . . . . . . . . . . . . . . . 100
3.5 OSIRIS-REx TAGSAM . . . . . . . . . . . . . . . . . . . . . . . . . . 107
3.5.1 Simulations . . . . . . . . . . . . . . . . . . . . . . . . . . . . 107
3.5.2 Penetration Depth . . . . . . . . . . . . . . . . . . . . . . . . 113
3.5.3 TAG Reconstuct . . . . . . . . . . . . . . . . . . . . . . . . . 117
3.5.3.1 Force Proles . . . . . . . . . . . . . . . . . . . . . . 117
3.5.3.2 S/C vertical speed . . . . . . . . . . . . . . . . . . . 120
3.5.3.3 Spring Compression Time . . . . . . . . . . . . . . . 122
3.6 Conclusions & Future Work . . . . . . . . . . . . . . . . . . . . . . . 122
4 Catastrophic Disruption of Gravitational Aggregates 126
4.1 Chapter Preface . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 126
4.2 Background . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 127
4.3 Methodology . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 132
4.3.1 Rubble-Pile Model . . . . . . . . . . . . . . . . . . . . . . . . 132
4.3.2 Simulation Parameters and Collision Geometries . . . . . . . . 135
4.4 Results . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 136
4.4.1 Head-on Equatorial Collisions . . . . . . . . . . . . . . . . . . 137
4.4.2 Oblique Equatorial Collisions . . . . . . . . . . . . . . . . . . 143
4.4.3 Head-on Collisions with  < 90 . . . . . . . . . . . . . . . . . 146
4.5 Discussion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 147
4.5.1 The \Universal" Law for Catastrophic Disruption . . . . . . . 147
4.5.2 Rotation Dependence of Catastrophic Disruption . . . . . . . 155
4.5.3 Dependence on Material Properties . . . . . . . . . . . . . . . 159
4.6 Conclusions & Future Work . . . . . . . . . . . . . . . . . . . . . . . 161
5 Large-scale Structure in Saturn's B ring 171
5.1 Chapter Preface . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 171
5.2 Background . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 171
5.3 Methodology . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 175
5.3.1 Periodic Boundary Conditions & Patch Properties . . . . . . . 175
5.3.2 Diagnostic Tools . . . . . . . . . . . . . . . . . . . . . . . . . 179
vi
5.3.2.1 Viscosity Calculation . . . . . . . . . . . . . . . . . . 179
5.3.2.2 Physical Optical Depth Calculation . . . . . . . . . . 182
5.4 Results . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 187
5.4.1 Correlation between Friction and Viscosity . . . . . . . . . . . 190
5.4.2 Correlation with Density . . . . . . . . . . . . . . . . . . . . . 192
5.4.3 Physical Optical Depth Measurements . . . . . . . . . . . . . 194
5.4.4 Comparison with Cassini UVIS occultations . . . . . . . . . . 196
5.5 Conclusions & Future Work . . . . . . . . . . . . . . . . . . . . . . . 198
6 Conclusions & Future Work 203
6.1 Conclusions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 203
6.2 Future Work . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 207
6.2.1 MASCOT lander . . . . . . . . . . . . . . . . . . . . . . . . . 208
6.2.2 The Plateaus in Saturn's C ring . . . . . . . . . . . . . . . . . 209
Bibliography 213
vii
List of Tables
1.1 High-Performance Computing Usage Summary . . . . . . . . . . . . . 17
2.1 Material Property Parameters . . . . . . . . . . . . . . . . . . . . . . 38
2.2 Values of test statistics: Table 1 . . . . . . . . . . . . . . . . . . . . . 48
2.3 Values of test statistics: Table 2 . . . . . . . . . . . . . . . . . . . . . 48
2.4 Regression analysis results . . . . . . . . . . . . . . . . . . . . . . . . 49
3.1 The physical properties of real granular materials. . . . . . . . . . . . 94
4.1 Summary of Mass Loss Results: Table 1 . . . . . . . . . . . . . . . . 165
4.1 Summary of Mass Loss Results: continued . . . . . . . . . . . . . . . 166
4.1 Summary of Mass Loss Results: continued . . . . . . . . . . . . . . . 167
4.1 Summary of Mass Loss Results: continued . . . . . . . . . . . . . . . 168
4.1 Summary of Mass Loss Results: continued . . . . . . . . . . . . . . . 169
4.2 Summary of 2 analysis . . . . . . . . . . . . . . . . . . . . . . . . . 170
5.1 Summary of ring patch parameters and simulation results . . . . . . . 188
5.1 Summary of ring patch parameters and simulation results, continued 190
viii
List of Figures
1.1 Mosaic of small Solar System Bodies . . . . . . . . . . . . . . . . . . 3
1.2 Granular Solid Liquids & Gas . . . . . . . . . . . . . . . . . . . . . . 4
1.3 Regolith on the Surface of Itokawa . . . . . . . . . . . . . . . . . . . 7
1.4 Regolith Mobility on Asteroids: Surface Accelerations . . . . . . . . . 11
1.5 Saturn Optical Depth Structure . . . . . . . . . . . . . . . . . . . . . 25
2.1 Initial set up for a periodic boundary simulation. . . . . . . . . . . . 35
2.2 Angle of repose example: Gravel . . . . . . . . . . . . . . . . . . . . . 38
2.3 Least square regression of mean rise speed: Earth Gravity Cases . . . 45
2.4 Regression line of mean rise speed: s and r in Earth Gravity . . . . 51
2.5 Convection Pattern . . . . . . . . . . . . . . . . . . . . . . . . . . . . 56
2.6 Rise Speeds in 102 g . . . . . . . . . . . . . . . . . . . . . . . . . . . 58
2.7 Regression lines of mean rise speeds: 104 g . . . . . . . . . . . . . . 60
2.8 Regression lines of mean rise speeds: 1 mm oscillation comparison . . 62
2.9 Height evolution of Brazil-nut: PBC vs. Walls . . . . . . . . . . . . . 67
2.10 PBC vs Walls: Clustering at dierent jump heights . . . . . . . . . . 70
2.11 PBC vs. Walls: bed height dilation at dierent local gravity . . . . . 74
3.1 TAGSAM diagram . . . . . . . . . . . . . . . . . . . . . . . . . . . . 89
3.2 TAGSAM diagram . . . . . . . . . . . . . . . . . . . . . . . . . . . . 90
3.3 Diagram of a uni-axial compression test . . . . . . . . . . . . . . . . . 95
3.4 Young's modulus measurement . . . . . . . . . . . . . . . . . . . . . 96
3.5 Young's Modulus and Poisson ratio as a function of friction properties 98
3.6 Comparison with Impact Experiments onto Dry Ottawa Sand . . . . 101
3.7 Comparison with Impact Experiments onto Glass Beads . . . . . . . 102
3.8 Gravity Dependence of the Penetration Depth . . . . . . . . . . . . . 103
3.9 The dependence of the drag coecient on gravity . . . . . . . . . . . 106
3.10 TAGSAM montage . . . . . . . . . . . . . . . . . . . . . . . . . . . . 108
3.11 Changes in spacecraft and TAGSAM separation . . . . . . . . . . . . 110
3.12 Force Prole on TAGSAM and resulting spacecraft dynamics . . . . . 111
3.13 TAGSAM penetration depth as a function of friction parameters . . . 114
3.14 Density maps of TAGSAM penetration . . . . . . . . . . . . . . . . . 116
3.15 Force prole for moderate friction cases . . . . . . . . . . . . . . . . . 119
ix
3.16 Force prole for moderate friction cases . . . . . . . . . . . . . . . . . 120
3.17 Spacecraft vertical speed as a function of material property . . . . . . 121
3.18 TAGSAM spring compression timings . . . . . . . . . . . . . . . . . . 123
4.1 Schematic of two collision scenarios . . . . . . . . . . . . . . . . . . . 134
4.2 Mass Loss for head-on Equatorial Impacts . . . . . . . . . . . . . . . 138
4.3 Map projection of mass loss . . . . . . . . . . . . . . . . . . . . . . . 140
4.4 Mass loss distribution with impact longitude . . . . . . . . . . . . . 141
4.5 Mass loss distribution with latitude . . . . . . . . . . . . . . . . . . . 143
4.6 Mass loss distribution for oblique equatorial impacts . . . . . . . . . . 144
4.7 Mass of the largest remnant vs. impact speed for pole-on impacts . . 148
4.8 Comparison with universal law for head-on collisions . . . . . . . . . 149
4.9 2 analysis for determining  and A . . . . . . . . . . . . . . . . . . . 153
4.10 QAR=Q
?
RD with adjusted interacting mass . . . . . . . . . . . . . . . . 156
4.11 Semi-Analytic t for Q?RD . . . . . . . . . . . . . . . . . . . . . . . . 160
4.12 Q?RD for dierent material properties . . . . . . . . . . . . . . . . . . 162
5.1 Diagram of a sliding-patch model for simulating a ring system . . . . 177
5.2 Viscosity Calculation: Comparison with Literature . . . . . . . . . . 183
5.3 Dependence of Optical Depth on Ring Geometry . . . . . . . . . . . . 185
5.4 Ring Simulation Snapshots . . . . . . . . . . . . . . . . . . . . . . . . 189
5.5 Viscosity and Toomre Q dependence with s . . . . . . . . . . . . . . 193
5.6 Dependence of Toomre Q on particle density and friction . . . . . . . 194
5.7 Physical optical depth variation with observing angle and friction . . 197
5.8 Comparison with Cassini UVIS occultations . . . . . . . . . . . . . . 199
6.1 MASCOT bouncing . . . . . . . . . . . . . . . . . . . . . . . . . . . . 210
6.2 Radial prole of Saturn's C ring . . . . . . . . . . . . . . . . . . . . . 211
x
List of Abbreviations
General:
pkdgrav parallel k-d tree gravity code
NASA National Aeronautics and Space Administration
ESA European Space Agency
JAXA Japan Aerospace Exploration Agency
ApJ Astrophysical Journal
MNRAS Monthly Notices of the Royal Astronomical Society
P&SS Planetary & Space Sciences
Chapter 1:
NEAR { Shoemaker Near Earth Asteroid Rendezvous { Shoemaker
NEA Near-Earth Asteroid
SSDEM Soft-Sphere Discrete Element Method
OSIRIS-REx Origins, Spectral Interpretation, Resource Identication, Security,
Regolith Explorer
TAGSAM Touch and Go Sample Acquisition Mechanism
Chapter 2:
BNE Brazil-nut Eect
PBC Periodic Boundary Condition(s)
Chapter 3:
S/C Spacecraft
Chapter 4:
SSSBs Small Solar System Bodies
Chapter 5:
UVIS Ultraviolet Imaging Spectrograph
VIMS Visual and Infrared Mapping Spectrometer
Chapter 6:
MASCOT Mobile Asteroid Surface Scout
DLR Deutches Zentrum fur Luft - und Raumfahrt
CNES Centre National d'Etudes Spatiales
xi
Chapter 1: Introduction
Granular material, in the form of regolith, is found in the uppermost layer of
most solid bodies in our Solar System, from planets and their satellites to asteroids
and comets. Regolith is loosely dened as broken up particulates that are present
on the surface of a moon or other small solar system body. In the planetary science
community, the terminology of sedimentary deposits is used to describe the range
of sizes of these particulates. The terms pebble, cobble, and boulder describe indi-
vidual particles of sizes ranging from 4 mm to 6.4 cm, 6.4 cm to 2.6 m, and > 2.6
m, respectively (Miyamoto et al., 2007). Before the rst spacecraft encounter with
small solar system bodies, it was generally assumed that most of these bodies were
monolithic rocks with bare surfaces (although, there was some discussion of alterna-
tives: see Dollfus et al. 1977; Housen et al. 1979; Michel et al. 2001). Because of the
low surface gravity on asteroids, it was assumed that regolith formation would not
be possible since small grains that were created during a cratering event on a small
body would then escape easily its gravitational in
uence (Chapman, 1976). How-
ever, the Galileo, NEAR and Hayabusa space missions revealed the asteroids Gaspra,
Ida, Eros (Sullivan et al., 2002; Robinson et al., 2002) and Itokawa (Fujiwara et al.,
2006) were covered entirely by a layer of regolith (see Fig. 1.1). Furthermore, the
1
regolith on the surface of these and other small bodies seemed to exhibit com-
plex dynamics. The crater morphology on the asteroid Lutetia as observed by the
Rosetta spacecraft is thought to be aected by dry granular 
ows and avalanches
on the surface (Vincent et al., 2012). Meter-size boulders resting on the surface of
tiny Itokawa pose a conundrum for the asteroid's surface evolution. Furthermore,
ground-based observations of near-Earth asteroids (NEAs) have revealed that the
weathered surfaces of some asteroids appear to be \freshened" possibly by close
encounters with planets (Binzel et al., 2010; Nesvorny et al., 2010). Astronomical
observations of small bodies are, by their nature, limited to only surveying a body's
surface. Therefore, understanding how the surface behaves is critical for interpret-
ing observations and for formulating coherent theories of their overall structure and
evolution.
Granular physics is a relatively new sub-discipline of physics that attempts
to combine principles that have been developed for both solid-state physics and
engineering (such as soil mechanics) with 
uid dynamics in order to formulate a
coherent theory for the description of granular materials, which are found in both
terrestrial (e.g., earthquakes, landslides, and pharmaceuticals) and extra-terrestrial
settings (e.g., asteroid and terrestrial planet surfaces, asteroid interiors, and plan-
etary ring systems). Granular matter behaves dierently than any other familiar
form of matter, such as solids, liquids or gases (Fig. 1.2). While the individual grains
are solid, collectively granular matter can also behave as liquids or gases. Although
constitutive equations exist for granular interactions on Earth, the inferred scaling
to the gravitational and environmental conditions on other planetary bodies such
2
Figure 1.1: Images of asteroids visited by spacecraft. Clockwise from the
top left: Eros, Itokawa, Vesta, Lutetia (not to scale). These asteroids
span a wide range in surface gravity. Due to their elongated shapes
and fast spins, small asteroids have surface gravities that can vary by
nearly an order of magnitude. Compared to Earth gravity, g, the mean
surface gravity on the small bodies shown here are 2.5102g, 5103g,
6104g, and 105g for Vesta, Lutetia, Eros, and Itokawa, respectively.
Source: Murdoch et al. (2015).
3
Figure 1.2: A pile of steel beads exhibiting avalanching behavior. Grains
in the lower right portion of the image are static, behaving as solids. A
phase transition occurs near the surface where grains begin to 
ow like
liquids. Slightly above the surface, ballistic grains can attain gas-like
qualities. source: Forterre & Pouliquen (2008).
as asteroids is currently untested. In the case of our Solar System, the growth of
this sub-discipline has been key in helping to interpret the formation, structure, and
evolution of environments that exhibit granular processes.
In this thesis, we apply the principles of granular physics to four dierent
planetary science scenarios. We study the size segregation of granular material in
an attempt to explain the presence of large boulders in spacecraft observations of
asteroids. We explore low-speed impacts onto the surface of asteroids, in order to
4
better understand how a spacecraft might be able to land or interact with the surface
of a small body. We study the collisions of rubble-pile asteroids, bodies that are
a gravitationally bound granular system, in contexts relevant to planet formation.
Finally, we study Saturn's B Ring, a unique granular system where the combined
eects of dissipative collisions and self-gravity lead to the emergence of large-scale
structures.
1.1 Evidence of Regolith and its Mobility on Asteroid Surfaces
Spacecraft visits to asteroids have discovered evidence of regolith on small
asteroids. Furthermore, it is likely that this surface material is able to be mobilized.
Surface images provide evidence of this motion, and understanding the underlying
forces behind these phenomena allows us to gain insight into how a small body
interacts with the space environment and how the asteroid itself functions as a
physical system. Some examples of environmental factors that can in
uence an
asteroid's structural evolution are micro-meteorite impacts, solar radiation pressure,
and tidal interactions with larger bodies.
Close-up images of Itokawa by the Hayabusa mission show that the surface is
covered with unconsolidated particles (Fig 1.3). These grains are typically found
piled on each other without being buried by a ne layer of sediment. Furthermore,
the position and orientation of the rocks appear to be stable against local gravity.
Also, small impact craters on the surface are found to have partially disrupted rims,
and the debris from these rims seems to drain towards the eective gravitational
5
potential lows of the asteroid (see the region highlighted by the white circle and ar-
rows in Fig. 1.3). These features indicate that grains on the surface of the asteroid
were likely reallocated after they were formed or reaccumulated on the surface. This
implies that the surface has been subjected to global vibrations. The source of these
vibrations could be micrometeoroid impacts (Garcia et al., 2015), tidal encounters
with larger objects in the Solar System (Yu et al., 2014), or thermal 
uctuations
due to solar radiation (Delbo et al., 2014). Focusing on impacts, we can use analyt-
ical formulations for the expected surface accelerations due to an impact to better
understand the phase space of possible seismic activity that can drive regolith mo-
bility. Following Miyamoto et al. (2007), the maximum half-cycle amplitude, A, and
seismic frequency, f , of a vibrating rock dene the average seismic strain energy per
unit volume of rock, e:
e = a
2f 2A2; (1.1)
where a is the bulk density of the asteroid. The maximum acceleration magnitude
that the rock experiences is:
a = 42f 2A: (1.2)
The total seismic energy density in the system is given as 2e (Lay & Wallace, 1995);
therefore, the total seismic energy, Es, of an asteroid with radius Da (assuming a
spherical shape), is given as:
Es =
aD
3
aa
2
48f 2
: (1.3)
For a seismic event that is generated by an impact, the kinetic energy, Ei, imparted
6
Figure 1.3: The Hayabusa mission's orbital survey showed regions on
the surface of Itokawa that were very smooth, occupied by cm-size par-
ticles, and regions of rough terrain, occupied by boulders. Here, we have
a closeup image of Itokawa's surface. The white circle indicates an im-
pact crater, which has a collapsed rim. The material seems to drain
towards the asteroid's eective gravitational potential low (indicated by
the white arrows).
7
by the event is given by (Richardson et al., 2004)
Ei =
1
12
iD
3
i v
2
i ; (1.4)
where i is the density of the impactor, Di is the diameters of the impactor, and
vi is the impact speed. Only a small fraction of the kinetic energy is transferred to
seismic energy. Therefore, the actual seismic energy can be written as
Es = Ei; (1.5)
where  is the seismic eciency factor. Using Eq. 1.3, 1.4, and 1.5, the maximum
acceleration can be written as
a = 2vif
s

i
a

Di
Da
3
: (1.6)
The gravitational acceleration on the surface of a spherical body, g, is given as
g =
2
3
GDaa; (1.7)
where G is the gravitational constant. Therefore, the maximum acceleration, a,
relative to the local gravity can be written as
a
g
=
3fvi
G
s

i
3a
D3i
D5a
: (1.8)
While this equation may be useful for estimating the magnitude of surface accelera-
tions due to seismic shaking for a monolithic spherical body, actual small bodies are
likely rubble piles, collections of gravitationally bound boulders. A rubble-pile as-
teroid can have micro- and macro-porous interior structure Richardson et al. (2002),
meaning their interiors have signicant void spaces. The asteroid Itokawa is thought
8
to be very porous as it has a very low bulk density of 1.9 g/cm3 compared to its
constituent rocks (which typically have densities > 2.7 g/cm3). Such an interior
structure would cause seismic waves to attenuate as they can be partially re
ected
and refracted o boulder-boulder and boulder-void boundaries. Therefore, a better
estimate of the surface accelerations on a small asteroid should include the eects of
seismic energy attenuation. A precise approach would require knowledge of the size
distribution of the constituent \rubble" of a rubble-pile. Analytically, an estimate
can be obtained by applying diusive scattering theory. Richardson et al. 2005 pro-
vide a prescription for including this eect, giving the seismic attenuation, At, of a
wave as it propagates across a body:
At = exp
fD2a
KQ

; (1.9)
where K is the seismic diusivity and Q is the seismic quality factor. Therefore,
the surface acceleration on the surface of a body due to an attenuated seismic wave
can be written as:
a
g
=
3fvi
G
s

i
rho3a
D3i
D5a

 exp
fD2a
KQ

: (1.10)
Most of the parameters in Eq. 1.10 are dicult to determine, especially since the
internal structure of any asteroid has never been directly measured through experi-
ments. However, some general properties of asteroids that may be directly measured,
such as size and bulk density, can give an indication of the general interior proper-
ties of an asteroid (i.e., whether it is likely to be a rubble pile or a solid monolithic
body).
9
In order to get a sense of how surface accelerations driven by impacts might
vary as a function of asteroid size and density, we adopt some nominal values for the
impact and seismic parameters (Richardson et al., 2004). Fig. 1.4 shows the surface
acceleration relative to gravity due to impacts (Eq. 1.10) as a function of the asteroid
diameter, for dierent combinations of f , Di and a. The shaded regions show
the magnitude of surface acceleration that is required to drive granular convection
(which occurs at values over unity, Jaeger et al. 1996) or to destabilize the surface
(which occurs at values over 0.2, Lambe & Whitman 1979). We see that the regolith
on smaller asteroids can be mobilized by smaller impactors. For an Itokawa-sized
body, an impactor that is 10 cm in size imparts sucient seismic energy to drive
granular convection on the asteroid's surface.
The analysis here has assumed seismic properties of asteroids that might not
actually re
ect reality. In order to gain a better understanding of the plausibility
of regolith motion by small impactors, a more accurate model of asteroid interiors
needs to be obtained. Garcia et al. (2015) performed simulations of seismic wave
propagation for dierent asteroid interiors and found that a naive approach, such
as the one described here, can underestimate the seismic accelerations by a factor
of  50. While this encourages the picture of seismically driven regolith motion,
direct simulations need to be performed to determine whether these seismic events
are suciently energetic to mobilize regolith such that there are observable changes
on the surface (this is the topic of Chapter 2).
10
Figure 1.4: We use nominal values of the seismic and impact parameters
to plot curves that may be representative of the expected surface accel-
erations for dierent asteroid sizes. We use values of Q = 200, K=0.25,
 = 104, vi = 4 km/s, and i = 2.5 g/cm3. Eros and Itokawa have bulk
densities of 2.7 g/cm3 and 1.9 g/cm3, respectively. We see that small im-
pactors (down to 10 cm) can cause regolith to be mobilized on Itokawa's
surface. For the assumed seismic and impact parameters, regolith on the
surface of Eros is unlikely to be mobilized by an impact.
11
1.2 Numerical Simulations of Granular Physics
It is dicult to develop a complete deterministic theory for granular materials
due to the fact that granular systems are composed of a large number of elements
that interact through a non-linear combination of various forces (mechanical, grav-
itational, and electrostatic, for example) leading to a high degree of stochasticity.
A tool for the study of these complex system has been numerical simulations that
model the dierent dynamical interactions (collisions, electrostatic forces, gravita-
tional forces, etc) directly. For a study of Solar System objects, numerical simula-
tions are invaluable as we can recreate extra-terrestrial conditions that would be very
dicult or impossible in a laboratory setting, including micro-gravity and vacuum
conditions. We approach our study with a combination of computational modeling
and analytical tools. Where appropriate, we use observational data to constrain
the very large parameter space involved with granular dynamics. The goal of each
of these projects was to study the behavior of granular material in extreme condi-
tions so that we may enhance our understanding of both astrophysical and granular
systems.
1.3 pkdgrav
All the simulations presented in this thesis were conducted with the parallel N -
body gravity tree code pkdgrav (Stadel, 2001), adapted to model collisions between
solid spherical grains (Richardson et al., 2000, 2009, 2011). Originally collisions
12
in pkdgrav were treated as idealized single-point-of-contact impacts between rigid
spheres. A soft-sphere option was added (Schwartz et al., 2012); with this option,
particle contacts can last many timesteps, with reaction forces dependent on the
degree of overlap (a proxy for surface deformation) and contact history. This allows
us to model multi-contact and frictional forces. The code uses a 2nd-order leapfrog
integrator, with accelerations due to gravity and contact forces recomputed each
step.
The soft-sphere discrete element method (SSDEM) implementation in pkdgrav
uses a spring/dash-pot to model the collisional forces between particles. In this
model, a spherical particle overlapping with a neighboring particle feels a reaction
force in the normal and tangential directions determined by spring constants (kn,
kt). The normal restoring spring force, FN , is generated according to Hooke's law,
FN;restoring = knxn^; (1.11)
where x is the amount of overlap between the two particles, and n^ is the unit vector
that points from a particle's center to its collider's center. The tangential restoring
force is given by
FT;restoring = ktS; (1.12)
where S is the tangential displacement from an equilibrium contact point, which is
tracked over the entire duration of an overlap. The value of kn determines the speed
of propagation of a sound wave through a granular system. A \correct" value for
kn would depend on the sound speed of the material being simulated. In practice,
a value of kn is chosen to ensure that overlaps do not exceed 1% of the smallest
13
particle in a collision. There are two regimes that constrain this determination: one
where the kinetic energy of the particles dominate, and another where the conning
pressure of surrounding particles dominates. For the rst regime, kn is determined
by considering the scenario where all of the energy of a collider with massm traveling
at a speed vmax is input into a single spring. Limiting the fractional overlap, xmax,
we obtain the following prescription for kn:
kn = m
vmax
xmax
2
: (1.13)
For the second regime, where the particle is in a dense static medium under the
in
uence of a global potential (such as gravity), the potential energy of the entire
medium is taken into consideration. In order to prevent overlaps exceeding values
of xmax, we require a spring with a strength that can resist the overburden pressure.
For a system with height H under the in
uence of a gravitational acceleration g,
and composed of particles with size s and density  that are packed so that the
porosity of the system is , kn is given by
kn =
gHs2
xmax
: (1.14)
In situations where both kinetic energy and conning pressure play a role, the largest
value of kn is the one chosen. Conventionally, the value of kt is typically set to be
 2
7
kn. The SSDEM implementation also allows optional damping eects. The
damping parameters (Cn, Ct) are related to the conventional normal and tangential
coecients of restitution used in hard-sphere implementations, "n and "t, by the
spring constants kn and kt (see Schwartz et al. 2012 for the full procedure). The
14
normal and tangential damping forces are given by
FN;damping = Cnun; (1.15)
FT;damping = Ctut; (1.16)
where un and ut are the normal and tangential components of the relative velocity,
respectively. We also include the eects that impose static, rolling, and/or twisting
friction. For the studies presented in this thesis, only the static and rolling friction
eects were used as these eects were found to dominate the behavior of particles in
the relevant simulations. The static and rolling friction components are parameter-
ized by dimensionless coecients s and r, respectively. The coecient of static
friction determines the maximum amount of tangential force that can be supported
by the contact point of two particles. This tangential force scales linearly with the
normal force at the contact point. In the event that the tangential force exceeds the
maximum value given by the static friction coecient, the particles slip past one an-
other, and the tangential force then is determined by the tangential spring constant
and the tangential damping coecient. The coecient of rolling friction is used to
determine the resulting torque due to two particles rolling against one another. Two
particles are said to be rolling with respect to each other if their relative velocities
are zero, but they have some relative rotational motion. The torque induced by the
rotation scales linearly with the coecient of rolling friction. In practical terms, a
non-zero coecient of rolling friction codies some non-spherical shape information
to our spherical particles. Plausible values for these various parameters are obtained
15
through comparison with laboratory experiments (Section 2.2).
This numerical approach has been validated through comparison with labora-
tory experiments; e.g. Schwartz et al. (2012) demonstrated that pkdgrav correctly
reproduces experiments of granular 
ow through cylindrical hoppers, specically the

ow rate as a function of aperture size, Schwartz et al. (2013) demonstrated suc-
cessful simulation of laboratory impact experiments into sintered glass beads using
a cohesion model coupled with the soft-sphere code in pkdgrav, and Schwartz et al.
(2014) applied the code to low-speed impacts into regolith in order to test asteroid
sampling mechanism design.
1.3.1 Computational Resources
Our simulations were run on high-performance supercomputing facilities that
are available to the Department of Astronomy at the University of Maryland, College
Park. These include the yorp cluster operated by the Department of AStronomy,
deepthought and its successor deepthought2 operated by the University of Mary-
land , and the MARCC bluecrab cluster administered jointly by the University of
Maryland and John Hopkins University and operated by John Hopkins University.
Table 1 summarizes the usage of each cluster according to simulations done for
each chapter by the total number of computational hours (service units, or SU's).
In general, simulations that required a low number of cores (< 32) or simulations
that were done to prototype new code implementations were performed on the yorp
cluster. Studies that required a large number of cores (> 32) or a large parameter
16
space sweep (such as the cases presented in chap. 2) were performed using one of
the supercomputing clusters (deepthought, deepthought2, or MARCC bluecrab).
Table 1.1: Usage summary by simulations done in each chapter, in kSUs (approxi-
mate).
yorp deepthought deepthought2 MARCC bluecrab
Chapter 2 300 200 600 -
Chapter 3 300 1200 2200 -
Chapter 4 500 400 - -
Chapter 5 100 - 1000 2000
1.3.2 Code Improvements: Fixed Search Ball
In this brief section, we detail code changes implemented during this thesis that
resulted in faster run times for simulations. In pkdgrav, each particle maintains a
list of its nearest neighbors in order to check for potential colliders. Rather than
performing a collision check for each other particle in a simulation, which would take
N2 time, a particle only performs a collision check on its nearest neighbors. The
tree code implemented in pkdgrav provides nearest-neighbor lists in order N log N
time. The number of nearest neighbors that a particle has to search for in order to
eectively determine whether it is in overlap with another particle largely depends
on its size with respect to other particles. For equal-sized particles, the theoretical
maximum number of neighbors that a single particle may be in contact with is 12.
(Note: this is only true if the particles are \just touching" and not overlapping
signicantly.) In general, for a size distribution of particles, the maximum number
17
of overlaps that a single particle may experience is given by
MAX NUM OVERLAPS PER PARTICLE = 8 + 4

Rmax
Rmin
2
; (1.17)
where Rmax and Rmin are the maximum and minimum particle radii in the simula-
tion, respectively.
Of vital importance to detecting neighbors, and therefore particle contacts,
is the use of a search ball by each particle. The search ball of a particle denes
the limited region around it where it will look for nearest neighbors. In the default
implementations of pkdgrav, the search ball is dynamically constructed each time
step. Years of practical use of the code have shown that this might not be the most
ecient way to nd nearest neighbors in all cases. Methods that generate a xed
search ball for each particle based on some prescription make the code run faster as
the search ball does not need to be re-constructed each time step. Three dierent
xed-search-ball strategies were implemented in order to speed up processing time
and ensure proper neighbor nding for certain pkdgrav scenarios. We outline these
methods below, and best practices for using them.
 Maximum diameter: The size of the search ball of each particle is set to the
diameter of the largest particle in the simulation. This option is optimal for
equal-size particles, or a size distribution of particles with a relatively small
size ratio (Rmax
Rmin
. 3). The advantage of this option is the speed-up in having
a xed search ball, while maintaining a relatively small memory requirement.
 Sum of particle radius and mean radius: The size of each particle's
search ball is set to its own radius plus the mean radius of all particles in
18
the simulation. This allows each particle to have a unique search ball that
depends on its own size. This option is especially useful for scenarios where a
single very large particle interacts with a large number of small particles (of
roughly equal size). This allows the large particle to eciently detect all the
small particles as nearest neighbors, while also limiting the size of the small
particles' search ball. This reduces the number of nearest neighbors that each
particle needs to check for overlaps. However, the disadvantage is the large
memory requirement. Each particle needs to carry around an overlap list that
can contain MAX NUM OVERLAPS PER PARTICLE number of particle IDs.
 Sum of particle radius and next-largest particle radius: The size of
each particle's search ball is set to its own radius plus the radius of the next-
largest particle in the simulation. This option is useful when there is a large
and continuous size distribution of particles, or in the case of a single large
intruder interacting with a continuous size distribution of particles. This op-
timizes the number of nearest neighbors a particle needs to check for overlaps.
The disadvantage is the large memory requirements.
In all cases (including the default) it is also necessary to sort the particles in
size order, either in increasing order for the default or maximum-diameter cases,
or decreasing for the other two. This is because to save memory, a given particle
only stores SSDEM contact information for neighbors with higher order numbers,
and neighbors are determined by measuring distances between particle centers, not
surfaces. For the increasing-size-order cases, the sort strategy ensures there is never
19
more than 12 possible contacts per particle (assuming overlaps are always small),
keeping memory usage low. For the decreasing-size-order cases, the sort strategy
ensures no contacts are missed (a large particle will see all possible neighboring small
particles), at the expense of a larger memory requirement.
The most obvious advantage of a xed search ball is that processors spend
less time re-computing search balls each time step. The cost of using a xed search
ball is larger memory requirements for certain cases. Practically, we found that
using a xed search ball speeds up the computational time by a factor of 2 for cases
where the simulated particles are of roughly equal size (within 10% in radius). For
simulations with a large size dierence between the smallest and largest particle
(e.g., for a bed of target particles that are 1/30th the size of an impactor), we nd
an order of magnitude or more speed up.
The advantage of using the default search ball option is that it is general, and
is low memory use (only 12 colliders ever need to be stored into memory). The
disadvantage to this method is that the search ball needs to be regenerated every
time step, and, for a very large size dierence in particles, each particle will have
to loop through MAX NUM OVERLAPS PER PARTICLE nearest neighbors to
check for overlaps, which increases as the square of the size ratio. These two factors
require a lot of processing time, and scale poorly with total number of particles.
20
1.4 Thesis Outline: Granular Physics in the Solar System
We used the code described above to study various scenarios related to Solar
System dynamics that are connected by the fact that they all exhibit granular
processes. Each study builds from and upon the others. Although these studies take
place at rather dierent scales, they illustrate the importance of granular physics
throughout the eld of planetary science.
1.4.1 Regolith Mobility on Asteroids
In Chapter 2, we present a study related to regolith motion on Asteroids. Size
segregation on asteroid surfaces can be driven by the the Brazil-nut Eect (BNE).
The BNE is a process whereby repeated vibrations of a granular system causes large
particles, embedded in a medium of smaller particles, to begin to rise against the
direction of the local eective gravity in the system. This eect has been invoked
in the past to explain the presence of large boulders on the surface of asteroids
(Matsumura et al., 2014). The BNE is a physical process that has been studied
extensively by granular physicists in the lab, and recently by planetary scientists
through simulations. The BNE study is presented rst as it is a relatively simple
problem to conceptualize that provides a good introduction to the complexities of
granular physics and its importance to the eld of planetary science. The work
presented here extends the current knowledge of the BNE by testing its eciency
for dierent boundary, gravity, and vibration conditions that closely resemble those
found on asteroids.
21
1.4.2 Spacecraft-Regolith Interactions
In Chapter 3, we discuss impacts onto a granular medium, and studies related
to the Touch-and-Go Sample Acquisition Mechanism (TAGSAM) on the OSIRIS-
REx spacecraft. In the last two decades, the exploration of small-solar-system bodies
has been advanced greatly by spacecraft missions that 
yby (e.g., Rosetta at Steins
and Lutetia), orbit (e.g., Dawn at Vesta and Ceres), touchdown on the surface (e.g.,
NEAR at Eros), and/or return samples (e.g., Hayabusa at Itokawa). In particular,
with the success of the Hayabusa mission, sample return has become an increas-
ingly viable and popular option for the study of asteroids. Besides the incredible
scientic value of a returned surface sample, the ability of such missions to directly
interact with the surface provides a unique opportunity to test and validate mod-
els of small-body surface geotechnical properties (Biele et al., 2015). While there
have been several missions that have touched down on a small-body surface (NEAR,
Hayabusa, and Rosetta's Philae lander), we still have a poor understanding of how
regolith behaves in low-gravity environments. The study of hypervelocity impacts
on asteroids (impacts at speeds of  1 km/s) has a rich literature (Michel et al.,
2001; Housen & Holsapple, 2011); however, very little is known about the response
of an asteroid surface to low-speed impacts ( 1 m/s). This dearth of knowledge is
due to the fact that there are no natural processes that produce  1 m/s impacts
on asteroids. Therefore, in Chapter 3, we study low-speed impacts onto a granular
medium by attempting to extend the current theory of impacts in Earth-gravity.
Then, we focus our attention on the wholly anthropocentric activity of landing
22
on a small solar-system body. Specically, we model the future touchdown of the
OSIRIS-REx spacecraft on the surface of the asteroid Bennu.
1.4.3 Catastrophic Disruption of Rubble-Pile Asteroids
In Chapter 4, we move from discussions of cm-size particles on the surface
or near-surface of small bodies to exploring the collisions of  1 km-sized rubble
piles, gravitational aggregates of meter size objects. There is strong observational
(Pravec et al., 2005) and theoretical (Richardson et al., 2002) evidence that some
asteroids in our Solar System are rubble piles rather than solid monolithic bodies.
Unlike monolithic objects, rubble piles have almost no tensile strength. A rubble
pile does have shear and compressive strengths. The relative importance of each of
these components depends on the surface properties of the constituent \particles"
of the gravitational aggregate. What little tensile strength a rubble pile has comes
from the cohesive forces between individual boulders. This cohesion likely comes
from short-range Van der Walls and electrostatic forces between the smaller particles
that make up the regolith layer on the surface of the boulders. The relatively recent
discovery of the rubble-pile structure of asteroids has had important implications for
understanding the evolution of asteroids undergoing physical processes that depend
on the strength of these objects. In particular, computational studies of rotational-
shape evolution of rubble piles (Cotto-Figueroa et al., 2015; Zhang et al., 2016) are
changing our understanding of how asteroids spin-up through thermal re-radiation
(i.e., the YORP eect, see Bottke et al. 2006) and how they might break apart
23
through rotational ssion (Walsh et al., 2008; Sanchez & Scheeres, 2012). Further-
more, the internal structure of a planetesimal (km-size objects found in the early
Solar System) aects the rate at which it can grow through mergers, an impor-
tant step for the formation of planets. Numerical simulations are indispensable
tools for probing the range of possible collisional outcomes (Leinhardt et al., 2000;
Leinhardt & Stewart, 2012). In Chapter 4, we study the binary collisions of rubble
piles in order to determine how the collisional parameters of the system aect the
outcome. We focus on how the rotational properties of the system aect the rate
at which planetesimals catastrophically disrupt (i.e., the rate at which a collision
would result in an object losing more than half of its initial mass). Since a rubble-
pile is inherently a gravitationally bound granular system, its behavior is dependent
on the physical properties of its constituent particles. Therefore, accurately mod-
eling rubble-pile dynamics requires carefully considering the material properties of
the individual particles. Hence, we use knowledge of material properties developed
by our study of granular physics on small scales (Yu et al., 2014) to better model
larger-scale systems. Symbiotically, we can use large-scale systems as astronomical
laboratories, where we can probe the physics of granular 
ow in extreme conditions.
1.4.4 Structure Formation in Saturn's B Ring
In Chapter 5, we continue our study of larger-scale granular systems by study-
ing Saturn's B ring. The rings are made up of meter-size particles composed mostly
of ice. The B Ring is one of Saturn's main rings, beginning  92,000 km away from
24
Figure 1.5: Optical depth pattern (top panel) and Cassini Image of
Saturn's main rings. The optical depth is derived from measurements
of stellar occultation showing the A,B, C rings and Cassini Division (10
km radial resolution). source: Colwell et al. (2009).
Saturn and extending outwards for  26,000 km. It is also the brightest of Saturn's
rings as it has the highest optical depths. Fig. 1.5 shows an image of Saturn's main
rings (upper panel) and their optical depth structure (lower panel).
The B ring is divided into four subrings, from the innermost B1 ring to the
outermost B4 ring. The normal optical depths range from  0.7 in the B1 ring to 
3 in the B2 ring, and transmitted light is completely blocked in the the B3 ring. The
Cassini mission has provided a wealth of high-resolution data that has revealed rich
and varied structure in the rings (such as propellers and ridges, see Tiscareno et al.
2010). These structures vary in scale, from thousand-km structure clearly visible in
images to small sub-km wake-like structure inferred through observations of stellar
25
occultations (Colwell et al., 2007). Wake-like structures are repeated variations in
the optical depth with radial distance. The rst indication of their presence was
from the observing-geometry dependence of the B ring optical depth as the observ-
ing angle changed. This occurs because low opacity regions are obscured by the
wake-like structures at certain observing angles. Two dierent mechanisms have
been suggested to explain the presence of wakes: self-gravity and viscous oversta-
bility. Self-gravity wakes are produced by the continuous gravitational collapse and
subsequent tidal disruption of ring-particle aggregates, causing the formation of
unique extended structures that are typically tilted with respect to the direction of
orbital motion. Viscous overstability results from the balance between the dissipa-
tion of energy due to ring particle collisions and the input of energy due to Keplerian
shear. This process leads to the formation of overdense structures that are perfectly
aligned with the direction of orbital motion. In Chapter 5, we present simulations
of the formation of km-scale structure in the B-ring, and attempt to clarify the
conditions necessary for the emergence of self-gravity wakes or viscous overstability
structure. We nd that the presence of each mechanism depends on the specic
physical properties (density, surface properties) of the ring particles. By comparing
with observational results, we are able to place constraints on the properties of the
actual ring particles. This allows us to place some constraints on the mass of the B
ring.
26
1.4.5 Outlook and Future Work
In Chapter 6, we conclude the thesis by summarizing the important ndings
of each of the granular-physics-related problems that we have studied. We also
give a brief overview of other collaborative projects that are underway that extend
and compliment the work presented here. Finally, we present perspectives for future
work and detail the improvements in simulation software and data analysis necessary
to build upon the results presented in this thesis.
27
Chapter 2: Size Segregation on Asteroid Surfaces: Brazil-nut Eect
2.1 Chapter Preface
The work presented in this chapter was published in the Monthly Notices of
the Royal Astronomical Society (MNRAS) as Maurel et al. (2017). The work was
done jointly by Clara Maurel (a PhD student now at the Massachusetts Institute
of Technology) and myself. Each of us gave a roughly equal amount of eort in
performing the simulations, analyzing the results, and preparing the manuscript for
publications. Clara's eorts focused on performing simulations in Earth-gravity and
conducting statistical analyses of these results. My contribution was updating the
code to function properly for Periodic Boundary Conditions (PBC)|see discussion
below|then performing simulations in low gravity. Furthermore, by performing
post-simulation analysis on the output data, I was able to show that dierent mech-
anisms drove the Brazil-nut Eect (BNE) for dierent boundary conditions, and
that these mechanisms depend on the local gravity. Roughly speaking, the rst half
of the chapter (pages 23{44), related to the procedures and generation of results
in Earth gravity, was Clara's work, while the latter half (pages 45{67), focusing on
analysis and discussion of results, was my work.
In this chapter, we investigate numerically size segregation among regolith
28
grains, which is likely to occur after repeated shaking events induced by micro-shocks
arising from micrometeorite impacts or other disturbances on asteroid surfaces. In
particular, we are interested in the so-called Brazil-nut eect (BNE), i.e., the migra-
tion of a large intruder toward the top of a vertically shaken granular system. We
go a step forward in simulating this segregation eect by implementing horizontal
periodic boundary conditions (PBC) in the N-body code pkdgrav, with the aim of
making the simulations more representative of the expected asteroid environment.
Since previous simulations used walls as a boundary condition, which don't exist on
asteroids. We study the in
uence of the PBC on the outcomes of a shaking simula-
tion in Earth gravity and compare them to those obtained within a walled container.
For both congurations, we also study the in
uence of static and rolling fric tion on
the BNE. We observe the well-known convection mechanism that drives the BNE
in the presence of walls. However, we nd that a dierent mechanism, consisting of
void lling, is responsible for the segregation of the \Brazil nut" with PBC, and we
discuss its relevance in light of previous granular physics studies. By running similar
simulations in 104g, we show that this void-lling mechanism remains relevant to
explain the BNE with PBC in a low-gravity environment. However, depending on
the gravity, the friction properties of particles in
uence the rise speed of an intruder
dierently. An increase in the rolling friction properties of the grains diminishes
the void-lling mechanism in Earth gravity, while an increase enhances it in 104g.
We speculate that this is likely due to a change in the granular 
ow timescales at
dierent gravity regimes.
29
2.2 Background
Particle segregation is one of the processes that can play a critical role in the
evolution of asteroid surface properties, and possibly internal ones (Murdoch et al.,
2015). Main Belt collisional evolution models suggest that a majority of asteroids
of size ranging from a few hundreds of meters to about 50 kilometers are of second
generation at least, i.e., are born as a result of the collision between two parent
bodies (Bottke et al., 2005, 2015). Numerical simulations of catastrophic disrup-
tion (Michel et al., 2001, 2015b) indicate that those second-generation bodies are
generally gravitational aggregates, formed by reaccumulation of ejected fragments,
and thus can be modeled as granular systems (Richardson et al., 2002). Moreover,
the vast majority of their surfaces seem to be covered in regolith, another kind
of granular medium, as suggested by thermal inertia measurements (Delbo' et al.,
2007; Delbo et al., 2015)) and direct images (e.g., on Eros, Asphaug et al. 2001, and
Itokawa, Miyamoto et al. 2007). The granular properties of regolith are still poorly
understood, as well as the diversity of these properties among the asteroid popula-
tion. For instance, the two near-Earth asteroids visited by spacecraft so far, i.e., the
17-km size Eros (NASA-Near-Shoemaker) and the 320-meter size Itokawa (JAXA-
Hayabusa), possess thoroughly dierent regolith properties, despite their same S
taxonomic type: a thick layer of ne grains was found on Eros, while a shallow layer
of gravel-like grains was found on Itokawa (Yoshikawa et al., 2015). In addition to
the asteroids' dierent dynamical and collisional histories, their dierent gravita-
tional environments likely have a strong in
uence on the variety of their surface
30
properties. However, both Eros and Itokawa show big boulders on their surfaces
(e.g., Barnouin-Jha et al. 2008) surrounded by smaller grains. This characteristic is
particularly striking on Itokawa, whose gravity is so small that the common invo-
cation of gravitational capture of ejecta after impact cratering can hardly explain
the presence of large rocks at the surface. In fact, ejection speeds of excavated ma-
terial are expected to exceed the escape speed of the asteroid (a few tens of cm/s).
Michel & Richardson (2013) proposed that the boulders at the surface of Itokawa
result from the asteroid's formation process: after the disruption of its parent body,
Itokawa could have formed by gravitational re-accumulation of ejected fragments.
However, regolith segregation and more specically the Brazil-nut eect or
BNE (Rosato et al., 1987) may be another process responsible for the presence of
big boulders at the surface of asteroids. The BNE consists of large particles rising
up through a granular medium composed of smaller grains, when the latter is re-
peatedly shaken. On an actual asteroid, the occurence of this eect could be the
result of seismic shaking induced by thermal cracking (Delbo et al., 2014) or mul-
tiple micro-impacts experienced by the asteroid over its lifetime (e.g., Cheng et al.
2002; Richardson et al. 2004). Seismic shaking, and the resulting motion of regolith
material, was invoked to explain the paucity of small craters on Eros and Itokawa
(Richardson et al., 2004; Michel et al., 2009) as well as the formation of ponds on
Eros (Cheng et al., 2002). In addition they may also cause particle segregation.
This could explain the presence of big boulders at the surface, and also imply that
an asteroid's interior could be composed of smaller particles than those observed at
the surface, suggesting a variation of macroporosity with depth. Therefore, parti-
31
cle segregation and BNE might have important implications in the arrangement of
asteroid surfaces and interiors.
Two issues arise with the BNE. First, the mechanism driving the BNE is still
not fully established and doubts remain about its unicity (Kudrolli, 2004). In par-
ticular, granular convection is a process often invoked as the origin of granular segre-
gation at the surface of asteroids (Miyamoto et al., 2007; Asphaug, 2007). However,
granular convection is strongly dependent on the gravitational acceleration and it
has been shown that a weak gravitational acceleration may reduce the eciency of
particle size segregation (Gray & Thornton, 2005; Murdoch et al., 2013). Moreover,
recent numerical simulations and parabolic 
ight experiments show that the rising
speed of a large intruder through a granular material decreases with the external
gravity (Tancredi et al., 2012; Guttler et al., 2013; Matsumura et al., 2014). There-
fore, asteroids with their low gravity should not be the most suitable place for such
a process. However, as the typical lifetime of asteroids of similar or larger size than
Itokawa is greater than tens of millions of years (e.g., Bottke et al. 2005), provided
that the meteoroid 
ux is ecient enough to induce repeated shaking events, particle
segregation may have time to occur despite the low-gravity conditions. The second
issue is that the BNE is known to occur in a conned, walled environment (such as
the one used by, e.g., Tancredi et al. 2012 and Matsumura et al. 2014), which does
not seem to be the best representation of an asteroid interior or near-surface. A
numerical study by Perera et al. (2016) showed that a rubble-pile asteroid, made up
of a population of particles of two dierent sizes, will show size segregation for suf-
ciently high impact energies. Their study showed that size segregation can occur
32
on a global scale without the need of a conning wall; however, their simulations
were limited to showing segregation of large 40-m and 80-m size particles. Thus,
although these experimental and numerical studies greatly helped to improve our
understanding of the size segregation process and its sensitivity to various parame-
ters (e.g., gravity, friction coecients), the conclusions need to be conrmed in an
environment more consistent with the one expected on an actual asteroid surface.
The aim of this chapter is to accomplish a step towards this goal by numeri-
cally investigating the BNE in an unconned environment, using periodic boundary
conditions, or PBC. Particle size segregation is, in this case, not in
uenced by the
connement caused by walls of the container and the potential occurrence of the
BNE will therefore only be due to the interaction among particles. If BNE occurs
under such conditions, the possibility of an extrapolation to the asteroid environ-
ment will be strengthened. We investigate the dierences in behavior of the large
intruder when wall boundaries are either present or absent, and also specically
focus on the in
uence of the static and rolling friction properties of the grains. We
also get closer to realistic asteroid environments by performing the same simulations
in low gravity, such as those encountered on Itokawa-like asteroids. Understanding
how a granular medium with a size distribution may segregate in an unconned
environment can allow us to better assess the porosity and grain distribution on
asteroids.
Today, we only have direct access to global bulk densities through space mis-
sions (e.g., Yeomans et al. 1997; Fujiwara et al. 2006), study of the Yarkovsky eect
(Chesley et al., 2003) or binary asteroid observations (Merline & Chapman, 2001).
33
This new aspect of investigation may be of great value for indirectly estimating
the density heterogeneity within asteroids. Considering the paucity of direct mea-
surements, the community needs other means to derive information from indirect
measurements, through, e.g., the modeling of dierent dynamical processes possibly
contributing to the evolution of asteroids' internal structure. This chapter oers a
contribution to this eort. In Section 5.3 we provide details on the various param-
eters investigated, with and without PBC. Results are presented and analyzed in
Section 5.4, and are followed by complementary analysis, interpretation and discus-
sion in Section 4.5. Conclusions and proposed future work are found in Section 2.6.
2.3 Methodology
2.3.1 pkdgrav: PBC
In our study, we make use of periodic boundary conditions (PBC). Doing so, we
can study particle segregation for an ensemble of particles that are freely vibrating
without a surrounding physical boundary. Periodic boundary conditions are imposed
on a patch of particles by replicating the patch in the horizontal directions (x-y),
perpendicular to the oscillating direction (z). Each replicated patch contains ghost
particles that match the relative positions of the original particles. Particles near the
boundary edges may interact in collision with ghost particles. Figure 2.1 illustrates
an initial set up with PBC.
34
Figure 2.1: Initial set up for a periodic boundary simulation. Central
pattern is surrounded by 8 replicas made of ghosts particles. Lines were
added to help visualize the replicated patterns. The colors of the real
particles are varied for better visualization.
35
2.3.2 Friction parameters
One aspect of this study is to analyze how material properties in
uence the
outcomes of the simulations. Given the strong uncertainties concerning the granular
properties of asteroids, we performed simulations using a broad spectrum of friction
coecients. Each grain in the system was simulated with the same set of coecients.
We made the choice to focus on materials of medium-to-high levels of static and
rolling friction, starting in the range where the BNE is expected to occur with wall
boundaries under Earth gravity, and pushing to levels where friction might impede
this eect.
Matsumura et al. (2014) studied the in
uence of friction on the rise of the
\Brazil nut" for static and rolling friction parameters, respectively s and r, be-
tween 0.0 and 0.9, using pkdgrav and cylindrical containers. Their results showed
that certain values of friction parameters were more favorable to the BNE: s & 0:5
and r . 0:2. They obtained these thresholds with coecients of restitution
n = 0:5, t = 0:5 and for an amplitude and frequency of oscillation A = 1 cm
and ! = 93:9 rad s1, respectively. We based the lower limit of the investigated
range of friction on a set of parameters identied by Matsumura et al. (2014) as
leading to the rise of a big intruder in an oscillating bed of small grains: s = 0:8,
r = 0:2 and medium dissipation of the grains with n = 0:5, t = 0:5. Note that
this set may not necessarily re
ect properties of a material existing on Earth. As
an upper limit, we dened a gravel-like material (Yu et al., 2014), characterized by
s = 1:3, r = 2:0, and again n = 0:5, t = 0:5. This very high friction aims
36
to represent the behavior of irregularly shaped grains that resist 
ow even though
we are using spherical particles. If we were able to simulate true angular shapes,
these coecients would likely be smaller compared to their current values. We also
implemented a \control" material, for which we did not expect to see the intruder
rise. For this we used a set of parameters mimicking the behavior of glass beads,
a low-friction material (s = 0:43, r = 0:1) with a very small amount of dissipa-
tion (n = 0:95, t = 1:0). This set of material parameters was calibrated against
avalanche experiments by Richardson et al. (2012).
These three sets of material parameters were compared to the literature by
carrying out angle-of-repose simulations. Grains of 0.5-cm radius were initially
contained, at equilibrium, in a rectangular box. Then, by removing one side of
the box, they were allowed to fall through an extended container, paved with a
layer of similar grains glued to the 
oor (to create a rough surface). We recorded
the evolution of the slope formed by the pile of particles and calculated the resulting
slope angle, which converged to the value of the angle of repose once the equilibrium
state was reached. Figure 2.2 shows the equilibrium state for the gravel-like material
and the corresponding angle of repose. The values of the angles of repose for the
dierent materials are summarized in Table 2.1. The angles obtained belong to
the ranges that can be found in the literature: () 2 [21,25] for glass beads (e.g.,
Barabasi et al. 1999) and () 2 [32,40] for gravel-like material (e.g., Orlova 1962).
We performed a series of simulations for equally distributed values of friction
parameters ranging between moderate and the gravel-like values: s 2 [0:8; 1:3] in
steps of 0.1 and r 2 [0:2; 2:0] in steps of 0.2, with n = 0:5, t = 0:5 throughout.
37
Figure 2.2: High-friction \gravel" material system at equilibrium after
the release of the grains. Here the nal achieved angle of repose is 37
(see Table 2.1).
Parameters Glass beads Moderate friction Gravel-like
n 0.95 0.5 0.5
t 1.0 0.5 0.5
s 0.41 0.8 1.3
r 0.1 0.2 2.0
 () 19.5 27.5 37.0
Barabasi et al. (1999) [21,25] - -
Orlova (1962) - - [32,40]
Table 2.1: Coecients of restitution and friction of the simulated glass beads, and of
the materials representing the lower and upper limits of the range of static and rolling
friction investigated. These are compared to the ranges found in the literature.
38
Each of these simulations was repeated 10 times (3 times for the cases in a low-
gravity environment|see Section 5.4 for details), with walls and PBC, each time
starting with a new random initialization to ensure the independence of the simu-
lations. In summary, we had a total of 60 simulations repeated 10 (or 3) times, for
both walls and PBC cases, in addition to the control simulation with glass beads.
Note | The rolling friction model used for this study is the same as that in
Yu et al. (2014) and Matsumura et al. (2014). The latest version of pkdgrav has
an optional new rolling friction model that uses a rolling force spring (Zhang et al.,
2016). This new model is particularly important to treat quasi-static problems. In
our case, however, the BNE simulations are strongly dynamic and the conclusions
should not suer from a change in the rolling friction model. We tested this by con-
ducting a limited number of runs comparing intruder rise speeds for the two models,
using roughly equivalent material parameters and found no substantial dierences.
We refrain from using this model for the majority of our simulations as it adds to
the computational cost of a simulation.
2.3.3 Simulation parameters and walls setup
In our simulations, two key parameters are the normal spring constant kn (see
Section 2.3.1) and the simulation timestep t. The value of kn was chosen so that
the maximum overlap of particles was no more than 1% of the smallest particle
radius for the entire simulation. The timestep was then chosen so that collisional
loading and unloading of the spring takes at least 30 steps. This means the timestep
39
is always very small compared to the gravitational dynamical time of the system. In
the present simulations, the acceleration of the bottom wall is several times the local
gravity, and constrains the required timestep and normal spring constant necessary
to properly resolve particle-particle and particle-wall collisions. The values of kn
and t are given by:
kn  m

vmax
xmax
2
; t 

m
kn
1=2
; (2.1)
where m is the typical mass of the particles, vmax is the maximum expected particle
speed and xmax is the maximum overlap when two particles collide, 1% of the smallest
particle radius in this case (see Schwartz et al. 2012). The maximum 1% overlap
criterion is sucient to properly resolve collisions when the maximum speeds in
the simulations are less than the sound speed of the simulated material. For our
simulations, the maximum expected speeds are less than a third of the sound speed
of the material. Therefore, by using a 1% maximum overlap criterion, we are able
to resolve collisions correctly and in a computationally ecient manner.
The small bed particles in our simulations are spheres with a normal distribu-
tion of radii of mean rp = 0.5 cm, with a standard deviation 0.05 cm, and a cuto at
 1 standard deviation. The Brazil nut has a radius of rBN = 1.5 cm. These values
were chosen for their convenience regarding computational time. Although BNE
may involve meter-size boulders on actual asteroids, we expect the eect to scale to
larger size, up to the limit where self-gravity becomes non-negligible. The density
is P = 2700 kg m
3 for each particle, meaning the large intruder is 27 times more
massive than a small particle in the granular bed. In order to estimate kn and t, we
40
used the mean bed-particle mass as the typical value for m and their mean size to
evaluate xmax. Also, vmax is the speed of the oscillating wall (i.e., vmax = A! = 93:9
cm s1; see below). From this we derived the following values for kn and t:
kn  5 105 kg s2 ; t  4 106 s: (2.2)
With pkdgrav, particles are allowed to interact with each other as well as
with walls. One of our main objectives is to compare walls cases and PBC cases. In
order to make walls and PBC setups as similar as possible, every wall uses the same
coecients of restitution and friction as the grains. In every walls case, lateral walls
form a parallelepiped open on the top of a 10  10 cm base, which is high enough
to prevent any escape of the particles once oscillations of the bottom plate start.
The choice of a rectangular container is dierent from the one of Matsumura et al.
(2014), who performed all their simulations with a cylindrical container. However,
as PBC require rectangular boundaries, the walls cases have to be performed with
a rectangular container to ensure a relevant comparison. In order to evaluate the
in
uence of the container's geometry, we carried out pairs of simulations aiming to
compare, under dierent conditions of friction, the outcomes of a simulation inside
a cylinder and inside a parallelepiped. We found that in all cases, both geometries
led to the same result (either no rise of the Brazil nut at all or full rise) within a
comparable amount of time. Consequently, we assumed that for our purposes, the
use of a parallelepiped container had no signicant in
uence on the outcomes of the
simulations compared to the use of a cylinder.
Finally, we set the oscillation of the bottom wall (and lateral walls if they
41
exist) to a frequency ! = 93:9 rad s1 and an amplitude A = 1 cm. These values
match the ducial choice of Matsumura et al. (2014) and represent a dimensionless
acceleration of  = A!
2
ag
 9 for simulations under Earth gravity (ag is the grav-
itational acceleration taken equal to 1g). This value of  will be kept the same
throughout the study, even for simulations under low gravity, as another means of
ensuring the validity of our comparisons.
2.3.4 Initialization process and main simulations
The initialization process was the same for the walls and PBC simulations.
In order to avoid any bias due to a uniform placement of the grains, we generated
randomly the particles inside an ellipsoidal volume large enough to avoid contacts
between the grains. This volume was suspended above the parallelepiped container
and particles were allowed to fall inside the container under Earth gravity. The
Brazil nut was positioned a few millimeters above the 
oor of the container and
small grains fell from a few tens of centimeters. This way, the large intruder was
at the bottom of the system at the end of the lling phase. Once the particles
were at equilibrium inside the container, we could either start walls cases or turn on
the PBC and remove the lateral walls. The central patch of particles, in the latter
case, was thus surrounded by a layer of replicas. The particles were then allowed to
re-equilibrate before starting the main simulations.
As mentioned in Section 2.3.3, we gave the bed particles a small size distribu-
tion. This prevented articial packing of the particles while lling the container or
42
during a main simulation. Fill-ups were performed separately for each conguration
of (s, r).
2.4 Results
In this section, we present the results of the various simulations described in
the previous section. We rst report on the results obtained under Earth gravity,
by analyzing the dierences in outcomes between PBC and walls cases, emphasizing
on the in
uence on the BNE of static and rolling friction levels. In the second part,
we present the results obtained in a low-gravity environment. Further analysis and
discussions about the results are found in Section 4.5.
2.4.1 Earth-gravity cases: Walls versus PBC
The principal observation that can be made while running identical simulations
performed with walls and PBC in Earth gravity is that nal outcomes of both
congurations are relatively similar. If the Brazil nut rises up with a walls case, it
also rises in the associated PBC case, but if it does not rise in certain conditions, it
will not rise with PBC either. For example, as already shown by Matsumura et al.
(2014) with walls, the simulation involving glass beads resulted in no net ascension
of the intruder, independent of the boundary conditions. Conversely, with grains
of moderate friction (s = 0:8 and r = 0:2), the BNE occurred with walls as well
as with PBC. Therefore, we conrm here that BNE does not necessarily require a
conned environment to occur. However, although the system achieves a similar
43
nal state in both cases, the mechanism driving the BNE is dierent. Here we
present the relevant results, which are analyzed and contextualized in Section 4.5.
As a means of comparison, we focus in this section on the rise speed of the Brazil nut.
Doing so, we can include in the analysis the cases where the intruder starts rising
but does not reach the top within the time of the simulation (arbitrarily xed). The
rise speed is measured in cm per oscillation (cm/osc). Recall from Section 2.3.2 that,
to provide better statistics, we performed 10 times independently the 60 simulations
derived from the chosen range of s values (from 0.8 to 1.3) and r values (from 0.2
to 2.0), each time with an independent initialization. Thus, for each value of s, we
had 10  10 = 100 cases (one for each value of r, for 10 simulations), and 10  6
= 60 cases for each value of r. In total, we had 600 simulations for both walls and
PBC congurations.
In this section, we are interested in comparing the global behavior of the walls
cases to that of the PBC cases. To do so, we focus on the evolution of the rise speed
with s for all values of r mixed up (i.e., independent of r), and similarly, the
evolution of the rise speed with r, for all the values of s mixed up. The question
of the in
uence of r for a given s and vice versa will be addressed in Section 2.4.2.
Figure 2.3 shows the least-squares regression lines obtained by tting the 600 values
of rise speed as a function of s (left panel) and r (right panel), for PBC and walls
cases. For reasons of clarity, the 1200 individual scatterpoints are not plotted, but
shaded areas show the 95% condence intervals for the ts.
At rst glance, the left panel of Fig. 2.3 suggests that, in both walls and
PBC cases, the higher the value of s, the higher the mean rise speed of the Brazil
44
Figure 2.3: Least squares regression lines showing the evolution of the
mean rise speed of the Brazil nut with s (left) and r (right), for walls
cases (dashed line) and PBC cases (solid line). The lines are obtained
by tting the data of the series of 10 simulations, all values of r and s
mixed up, respectively for left and right panel. Shaded areas show the
95% condence intervals for the ts.
45
nut. In addition, this growing trend is similar in both cases. On the right panel of
Fig. 2.3, the opposite trend appears. The mean rise speed of the intruder decreases
with an increasing value of the rolling friction coecient r. This time, a steeper
slope is evident for PBC cases, meaning that with PBC, the rise of the intruder
is more aected by the level of rolling friction than with walls. Finally, looking at
both panels, the mean rise speed of the Brazil nut with PBC seems to be globally
higher than with walls. This last observation would suggest that either the same
mechanism is driving the BNE in both cases, but it is more ecient with PBC,
or two dierent mechanisms are actually leading the BNE with walls and PBC.
The dierence in steepness of the slopes in the right panel of Fig. 2.3 would actually
support the second hypothesis. However, before elaborating on this, we performed a
series of statistical tests on the raw data to strengthen our forthcoming conclusions.
For this, we used a variation of the Student t-test adapted to the study of the slope
of regression lines. The objective was threefold: 1) ensure that the trends visible
on Fig. 2.3 are not artifacts due to a linear regression made on dispersed data; 2)
check whether the two slopes in the left panel of Fig. 2.3 are statistically equal and
whether the two slopes of the right panel are signicantly dierent; 3) assess whether
the intercepts of the PBC and walls lines are statistically dierent from each other.
We will rst present the three tests and their results successively, before elaborating
on the actual information we can get from these results.
Test 1) requires comparing the four slopes of linear regression in Fig. 2.3 to
zero. If the corresponding test statistics are higher than a critical value, the null
hypothesis is rejected, and the slopes would indeed have a signicant steepness. To
46
compare a slope to zero, the following test statistic can be used:
t =
j b j
sb
; (2.3)
where b is the slope of the regression line considered, and sb is a function of the
standard deviation of the abscissa values (s or r) and of the standard error of
the estimate, which measures the accuracy of predictions made by the regression
line. This t-value is then compared to a critical value tcrit found on Student's law
table, which depends on the number of available data and on the desired condence
interval. Here, for 600  2 available data and a 95% condence interval, Student's
law table gives tcrit  1:96. The four calculated values of t are summarized in
Table 2.2 (two for the evolution of the mean rise speed in walls and PBC cases with
s, and the two others with r). Each t-value of Table 2.2 being greater than tcrit,
this test ensures that statistically, the mean rise speed of the Brazil nut is indeed
aected by an increasing level of static and rolling friction, whether walls or PBC
are turned on. However, Test 1) does not compare walls and PBC cases to each
other.
Test 2) requires comparing the slopes of the PBC and walls regression lines to
each other. The test statistic that can be used is given by
t =
j bw  bpbc jq
s2bw + s
2
bpbc
: (2.4)
Here, bw=pbc is the slope of the corresponding regression line (walls or PBC), and
sbw=pbc is the corresponding value of sb. The critical value remains tcrit  1:96, and
47
tcrit = 1.96
Walls cases PBC cases
RS = f(s) RS = f(r) RS = f(s) RS = f(r)
t-value 6.04 9.39 3.34 9.92
Table 2.2: Values of the test statistic obtained from Eq. 2.3. From left to right, the
two rst columns correspond to walls cases, and the other to PBC cases. The left
column of each set of two (labeled \RS = f(s)") gives the t-value corresponding
to the slope of the regression line shown in the left panel of Fig. 2.3. This t-value
shows how signicantly the walls and PBC cases are aected by an increase of the
static friction level. The right column (labelled \RS = f(r)") gives the t-value
corresponding to the lines in the right panel of Fig. 2.3, which characterize this time
the eect of rolling friction level on the walls and PBC cases. The critical t-value
for 95% condence, tcrit, is given at the top of the table. The t-values are underlined
if they are above tcrit.
tcrit = 1.96
RS = f(s) RS = f(r)
t-value 0.17 3.52
Table 2.3: Values of the test statistic obtained from Eq. 2.4. These values show how
signicantly the evolution of the mean rise speed with respect to s (left column)
and r (right column) diers from PBC to walls cases. The critical t-value for 95%
condence, tcrit, is given at the top of the Table. The t-values are underlined if they
are above tcrit.
Table 2.3 shows the two t-values that compare the evolution of mean rise speed with
s (left) and r (right) between walls and PBC cases. As expected from Fig. 2.3,
this test ensures that statistically, the mean rise speed of the Brazil nut is similarly
aected by the level of static friction (the slopes are statistically equal), whether
walls of PBC are turned on (ts < tcrit). On the other hand, the eect of r on the
mean rise speed of PBC cases is signicantly dierent for walls and PBC cases.
To complete the set of statistical verications, we want to assess whether the
intercepts of the regression lines in each case (left and right panels of Fig. 2.3)
are signicantly dierent. For this, we subtract one regression line from the other
48
tcrit = 1.96
RS = f(s) RS = f(r)
Condence interval [0:05 ; 0:02] [0:05 ; 0:002]
Table 2.4: For each panel of Fig. 2.3 (\RS = f(s)" and \RS = f(r)") is given
the condence interval for the value of the intercept of the line corresponding to the
walls-cases regression line subtracted from the PBC-cases regression line. The value
zero belonging to both intervals, we cannot conclude that the intercepts of the walls
and PBC cases are signicantly dierent.
and nd the condence interval of the resulting predicted value when x = 0 (y-
intercept). If 0 belongs to the condence interval, it means that the intercept of the
line resulting from the subtraction can be zero, and therefore that the two intercepts
cannot be signicantly distinguished. Table 2.4 shows the condence intervals for
each panel of Fig. 2.3. Here, the visual interpretation was actually misleading, and
we cannot say that the mean rise speeds of PBC cases are statistically higher than
that of walls cases.
To summarize, Fig. 2.3 visually suggests that BNE could, under Earth gravity,
be driven by a dierent mechanism with PBC compared to with walls. Our results
suggest that whatever the mechanism is that drives the BNE for PBC, it has a
similar dependence on the static friction as the mechanism driving the BNE with
Walls. However, the BNE mechanism involved in PBC has a stronger dependence
on the rolling friction than the BNE mechanism involved with Walls. The previous
statistical tests, however, help us to qualify these rst conclusions. So far, we know
49
that, for both PBC and walls cases, the mechanism driving the BNE is sensitive to
the level of static friction (the mean rise speed increases with s) and the level of
rolling friction (the mean rise speed decreases with r). Moreover, with PBC, the
in
uence of r is signicantly more important than with walls, whereas the in
uence
of s is at rst glance similar. Finally, we cannot conclude from the present analysis
that the BNE mechanism with PBC is making the intruder rise faster than with
walls, as the intercept of the regression lines on Fig. 2.3 are actually not signicantly
dierent.
At this point, further analysis is needed to discriminate the mechanisms driving
BNE with PBC and walls, provided that they are dierent. For this, we no longer
focus on the mean rise speed, but on the in
uence of the static and rolling friction
at given values of r and s respectively.
2.4.2 Earth-gravity cases: static and rolling friction
Figure 2.4 shows the regression lines obtained by tting the data of the series
of 10 simulations for PBC cases (top) and walls cases (bottom). However, unlike
Fig. 2.3 where all values of friction parameters are mixed up, each line represents
a certain level of rolling friction (left panels) and static friction (right panels). For
clarity, the 95% condence intervals are omitted from these plots. For PBC cases,
these intervals are globally  0.02 for the top left panel and  0.01 cm/osc for the
top right panel. For Walls cases on bottom left panel, we nd a maximum of  0.04
cm/osc for condence interval of the r = 0:2 tted line, while for all other lines we
50
0.8 0.9 1.0 1.1 1.2 1.3
Static friction coefficient, µs
0.00
0.02
0.04
0.06
0.08
0.10
0.12
0.14
R
is
e
 s
p
e
e
d
 (
cm
/o
sc
)
PBC
µr  = 0.2
µr  = 0.4
µr  = 1.0
µr  = 1.4
µr  = 2.0
0.5 1.0 1.5 2.0
Rolling friction coefficient, µr
0.00
0.02
0.04
0.06
0.08
0.10
0.12
0.14
R
is
e
 s
p
e
e
d
 (
cm
/o
sc
)
PBC
µs  = 0.8
µs  = 1.0
µs  = 1.3
0.8 0.9 1.0 1.1 1.2 1.3
Static friction coefficient, µs
0.00
0.02
0.04
0.06
0.08
0.10
0.12
0.14
R
is
e
 s
p
e
e
d
 (
cm
/o
sc
)
Walls
µr  = 0.2
µr  = 0.4
µr  = 1.0
µr  = 1.4
µr  = 2.0
0.5 1.0 1.5 2.0
Rolling friction coefficient, µr
0.00
0.02
0.04
0.06
0.08
0.10
0.12
0.14
R
is
e
 s
p
e
e
d
 (
cm
/o
sc
)
Walls
µs  = 0.8
µs  = 1.0
µs  = 1.3
Figure 2.4: Regression lines obtained by tting the data of the series of
10 simulations for PBC cases (top) and walls cases (bottom). Left panels
represent the evolution of the rise speed with s for dierent values of
r, while right panels show the evolution with r for dierent values of
s.
51
have  0.01 cm/osc. On bottom right panel, we also nd a larger condence interval
for the s = 1:3 curve (0:02 cm/osc) and 0:01 cm/osc otherwise. Interestingly,
Figure 2.4 seems to imply an important dierence in behavior of the Brazil nut
between walls and PBC cases. Even in the worst case scenarios, the larger intervals
found for Walls cases do not change the visible trends.
Let us rst focus on walls cases (bottom panels). On the left panel, we see
that the lower the value of r, the more important the in
uence of s on the rise
speed of the intruder. For example, for r = 0:2 (solid line), the least-squares linear
regression shows that the rise speed can increase from approximately 0.02 cm/osc
at s = 0:8 to more than 0.12 cm/osc at s = 1:3, while for r = 2:0 (dash-dot-dot
line), the in
uence of s on the rise speed seems to be insignicant. Moreover, this
trend is strongly correlated with the values of r. The in
uence of s decreases
rapidly with increasing r, and then stagnates at almost zero. A clear trend is also
visible in the bottom-right panel of Fig. 2.4. At low levels of static friction, the
in
uence of r on the rise speed is hardly visible (solid line), however at high s
(dash-dot line), the rise speed can decrease by approximately a factor of 10 between
r = 0:2 and r = 2:0, again with a correlation between the values of s and the
importance of its in
uence on the rise speed. These observations were validated by
the same statistical test as Test 2), in which we compared the slopes two by two.
The results of this test were in agreement with the trend visible in Fig. 2.4.
On the other hand, the in
uence of static and rolling friction levels in PBC
cases is not so clearly brought out by the top panels of Fig. 2.4. Neither the in-

uence of s on the rise speed of the intruder nor the in
uence of r seem to be
52
correlated with an increase or decrease of r and s, respectively. In fact, the slope
of each regression line in the left and right panels looks relatively similar. This was
conrmed by performing Test 2) again: no clear trend and no signicant dierence
between a majority of slopes was found. In addition, the values of the intercepts of
the regression lines are also not clearly correlated to an increase of r or s. Fig-
ure 2.4 might suggest the existence of an eciency threshold for values of r and s
respectively around 0.4 and 1.0 (which give the highest intercept), but more data
would be required to conrm or refute this hypothesis.
In light of the above-mentioned results, we claim that a dierent mechanism
is driving the BNE with walls compared to with PBC. In fact, with walls, the
eciency of the process (characterized here by the rise speed of the Brazil nut)
is highly dependent on the values of r and s: at low rolling friction levels, the
process is catalyzed by an increasing level of static friction. However, a higher level
of rolling friction completely inhibits this eect. In contrast, at low levels of static
friction, the in
uence of r on the rising process is almost non-existent, whereas for
a higher s, an increase of the rolling friction level impacts drastically and negatively
the rise speed of the intruder. It is well known that, under certain conditions (level
of friction, frequency and amplitude of the oscillations), friction between particles
and walls engenders a convective motion (e.g., Knight et al. 1993), with particles
rising through the center of the container and sinking along the walls. The results
presented above seem to adequately characterize such a mechanism. The convective
motion needs a certain level of static friction between the walls and the small grains
to take place. As this level gets higher, the convection rolls form more and more
53
easily (very likely until a certain limit, apparently not reached in our simulations).
However, all this depends on the level of rolling friction. If particles are hardly able
to roll on each other, the whole inward stream (inside the layer of particles that are
in contact with the walls) is inhibited and the intruder rises more slowly.
With periodic boundary conditions, the very dierent behavior of the Brazil
nut for the same static and rolling friction conditions strongly suggests that BNE
is driven by a mechanism dierent from convection. This mechanism is marginally
aected by the level of static friction but is, however, more dependent on the level of
rolling friction. In addition, the top panels of Fig. 2.4 show a more stochastic behav-
ior of the rise speed than for walls cases, for which the correlation between rise speed,
s and r is demonstrated. To further strengthen the suspicion that convection is
not the leading mechanism for BNE when PBC are turned on, we analyzed the 
ow
of particles for several simulations, with walls and PBC, as illustrated in Fig. 2.5.
This gure was created by considering a 5-particle-wide cross-section of the center
of the container (the selection was done along the y-direction). A grid consisting
of cells that are 1 particle wide is superimposed on the x-z plane. For every cell,
the average displacement of particles from that cell to adjacent cells is tracked and
plotted. In order to properly visualize the 
ow of particles, it is necessary to isolate
the relative displacement of a particle with respect to other particles from the over-
all shaking motion of the entire medium. This is accomplished by only considering
particle positions at the same phase of an oscillation cycle. We chose the zero phase
of the oscillation, where the bottom plate has zero displacement and is moving in
the positive z-direction. Figure 2.5 shows a snapshot of this average motion of par-
54
ticles. The right panel, representing a walls case, clearly shows a convective motion
of particles, as a thin particle-wide layer is forced downwards along the walls. This
forces the center of the bed and the Brazil nut to move upwards. This pattern of
granular convection is apparent throughout all walls simulations. The left panel
does not show a similar convective motion for the PBC cases. Rather, we nd that
most PBC cases show more-randomized and less-uniform motion. There are some
cases where the convective pattern visible on the right panel is seen intermittently
in PBC cases; however, it only lasts for a few forcing cycles before it disappears.
Hence, due to the absence of boundaries, convection cycles do not appear strongly
when PBC are turned on, and despite that, we do see the Brazil nut rising. A
further analysis of the leading mechanism for BNE in the absence of conning walls
is proposed in Section 4.5.
2.4.3 Low-Gravity Cases
A key objective was to see how our ndings might extend to a realistic asteroid
gravity environment. For this we performed a suite of simulations with a local
gravity set to 104 times Earth gravity (104g), which is equivalent to the surface
gravity on an Itokawa-sized body (Michikami et al., 2008). We re-initialized the
grains in the low-gravity environment and allowed them to settle as explained in
Section 2.3.4. This was again done for walls and PBC cases. To accurately compare
these low-gravity simulations to those done in 1g, we modied the shaking frequency
of the granular bed such that  remains the same (  9; see end of Section 2.3.3).
55
4 2 0 2 4
x (cm)
0
2
4
6
8
10
12
14
z 
(c
m
)
PBC
4 2 0 2 4
x (cm)
Walls
Figure 2.5: By tracking the change in particle positions over a single
forcing period, we can visualize the convective motion of grains in sim-
ulations with wall boundaries (right) and the lack of convective motion
in simulations with periodic boundaries (left). Both of these cases are
for grains with friction parameters s = 0.9 and r = 1.4 (relatively high
friction).
56
As a consequence, we rst adopted a shaking frequency of 0.939 rad s1, while
keeping the amplitude of the oscillation at 1 cm. This allowed us to use a timestep
of t  3.7  104 s for these low-gravity simulations. We maintained a maximum
possible fractional overlap of  1% by adjusting the value of kn to  42 kg s2. The
value of kt was then
2
7
kn  12 kg s2. This adjustment in the values of the spring
constants is done for the sake of computational eciency. If we were to maintain
the same value of kn as that used in the Earth-gravity simulations, the total time
taken to complete each of the low-gravity simulations would increase by 2 orders of
magnitude.
In order to ensure that this change in the value of kn does not have a con-
siderable impact on the outcome of our simulations, we ran a suite of simulations
in a local gravity of 102g, where the value of  is again kept constant. Half of
the simulations were done with a kn and kt that is the same as that used in the 1g
cases, and the other half were done with a kn and kt chosen to maintain a maximum
possible overlap of 1% (based on a vibrating bed with an oscillation frequency of
9.39 rad s1 and an oscillation amplitude of 1 cm). The local gravity of 102g was
chosen with consideration for time constraints, since the slower simulations (those
with spring constants similar to the 1g cases) were only 1 order of magnitude slower
than simulations with spring constants chosen to maintain 1% overlaps. We per-
formed simulations with wall boundaries and chose three dierent values of s (0.8,
1.0, 1.2) and ve dierent values of r (0.2, 0.6, 1.0, 1.4, 1.8). This small suite of
simulations allows us to eciently determine the magnitude of kn's in
uence in the
rest of our simulations. The results of this test are shown in Fig. 2.6. Here we show
57
Figure 2.6: Rise speeds of the Brazil nut for simulations with a local
gravity of 102g. The open symbols represent simulations where the
particle collisions are resolved with spring constants similar to those in
the 1 g simulations, with s = 0.8, 1.0, and 1.2 represented by diamonds,
squares, and circles, respectively. The spring constants and the time-
step used in these simulations are very conservative values. The lled
triangles represent simulations where the spring constants are chosen
based on maintaining a maximum possible fractional overlap of 1%, with
s = 0.8, 1.0, and 1.2 represented by up-, down-, and left-facing lled
triangles, respectively. The time-step varies by 1 order of magnitude
between these two cases (with the lled symbols having a greater time-
step).
58
the variation in the time taken for the intruder to rise to the top of the granular
bed. We nd that for a wide range in r, the rise speeds for various s have similar
magnitudes and variances for the two spring constant cases. For both cases, we nd
that the \brazil" nut rises with fairly similar speeds, the cases with a smaller kn and
larger time-step (lled triangles) bracketing the cases with a larger kn and smaller
time-step (open symbols). We conclude that, while there is some possible eect in
the change of the particle spring constants, this eect is quite small. Furthermore,
we nd that both cases lead to similar trends across friction parameters, e.g., under
102 g, the rise-speed increases as s increases, independently of r. Therefore, us-
ing a spring constant that allows for a faster computational time, while maintaining
a small maximum fractional overlap, is acceptable for the current simulation setup.
For the cases with a local gravity 104g, the range of considered friction levels
was similar to that described in Section 2.3.2: s ranged from 0.8 to 1.3, and r
from 0.2 to 2.0. Again, we attempted to minimize any stochasticity that might arise
from sample preparation by performing three separate simulations (with separate
initializations) for each s and r pair. In total we performed 360 simulations in the
low-gravity regime (3 initializations, 2 boundary conditions, and 60 unique friction
pairs).
Overall, for a shaking amplitude of 1 cm and both PBC and walls scenarios,
we were still able to observe the BNE occurring in almost all the simulations. Every
walls case showed the intruder rising to the top of the granular bed in the time of
the simulation, whereas for PBC simulations, it occurred in approximately 84% of
the cases (152 out of 180). The PBC cases that did not lead to a complete rise
59
0.8 0.9 1.0 1.1 1.2 1.3
Static Friction coefficient, µs
0.00
0.02
0.04
0.06
0.08
0.10
0.12
0.14
R
is
e
 s
p
e
e
d
 (
cm
/o
sc
)
PBC
µr = 0.6
µr = 0.8
µr = 1.0
µr = 1.4
µr = 2.0
0.6 0.8 1.0 1.2 1.4 1.6 1.8 2.0
Rolling Friction coefficient, µr
0.00
0.02
0.04
0.06
0.08
0.10
0.12
0.14
R
is
e
 s
p
e
e
d
 (
cm
/o
sc
)
PBC
µs = 0.8
µs = 1.0
µs = 1.3
0.8 0.9 1.0 1.1 1.2 1.3
Static Friction coefficient, µs
0.0
0.1
0.2
0.3
0.4
0.5
R
is
e
 s
p
e
e
d
 (
cm
/o
sc
)
Walls
µr = 0.2
µr = 0.4
µr = 1.0
µr = 1.4
µr = 2.0
0.5 1.0 1.5 2.0
Rolling Friction coefficient, µr
0.0
0.1
0.2
0.3
0.4
0.5
R
is
e
 s
p
e
e
d
 (
cm
/o
sc
)
Walls
µs = 0.8
µs = 1.0
µs = 1.3
Figure 2.7: Rise speeds of the Brazil nut for simulations with a local
gravity of 104g. The panels show lines that represent the linear trend
of the rise speed with s for dierent values of r (left) and vice versa
(right) with data extracted from three independent suites of simulations.
The top panels show the results for cases with PBC, while the bottom
panels show the results for walls cases.
60
tended to be cases with low values of the rolling friction parameters (r < 0.6).
In addition, these low-friction cases that did not result in a BNE tended to have
the intruder begin to rise, but eventually stall at a system height lower than the
barycenter of the system. These cases were omitted from subsequent analysis of rise
speeds. Figure 2.7 shows the average rise speeds for PBC (top panels) and walls
(bottom panels) cases as a function of s (left panels) and r (right panels). The
curves represent linear ts of 3 independent simulations (3 dierent initializations of
the granular bed). Like previous studies (Tancredi et al., 2012; Guttler et al., 2013;
Matsumura et al., 2014), we nd that, at moderate friction values, the rise speeds
in a lower-gravity environment are comparable to those under Earth gravity. Here,
we also show that this result still holds even when PBC are used. However, unlike
these previous studies, we have performed an exploration of the friction parameter
(a proxy for material property) phase-space in low gravity. This has resulted in
some surprising outcomes.
For the PBC cases, we nd that the magnitude of rise speeds in 104g are
similar to the 1 g cases. Furthermore, similar to Earth gravity, we nd that both the
walls and PBC cases show a positive correlation with s (see left panels of Fig. 2.7).
However, the walls cases seem to have a stronger correlation with s than the PBC
cases. This likely points to a dierence in the mechanism driving the BNE in each
case, which is discussed further in Section 4.5. Furthermore, contrary to the Earth-
gravity cases, we nd that the rise speed in both PBC and walls cases is positively
correlated with rolling friction, r (compare right panels of Figure 2.4 to Figure 2.7).
This positive correlation appears to be stronger for the PBC cases than the walls
61
0.8 0.9 1.0 1.1 1.2 1.3
Static Friction Coefficient µs
10-3
10-2
10-1
100
R
is
e
 S
p
e
e
d
 (
cm
/o
sc
)
PBC
Walls
Walls 1 mm
0.5 1.0 1.5 2.0
Rolling Friction Coefficient µr
10-3
10-2
10-1
100
PBC
Walls
Walls 1mm
Figure 2.8: Average rise speeds of the Brazil nut for all simulations with
a local gravity of 104g. The left and right panels show the rise speed
as a function of s and r, respectively. Mean rise speeds are shown for
simulations with an oscillation amplitude of 1 cm with PBC (solid line)
and walls (dashed line). Oscillation amplitudes of 1 mm are also shown
for walls cases (dot-dashed).
62
cases. For these PBC cases, we speculate that, in a low-gravity environment, the
grains are lofted over a much longer period than under Earth gravity. This larger
lofting period may give the grains more time to recongure, allowing for the entire
system to re-settle with more voids (i.e., a smaller lling factor), as particles with
a larger rolling friction can more easily maintain larger void spaces (in the other
extreme, frictionless particles are more easily compacted into tighter spaces). As we
will discuss in Section 4.5, larger void spaces in a granular system are important for
the BNE to occur in a PBC system.
Previous studies (e.g., Matsumura et al. 2014) showed that the rise speed is
generally constant across gravitational regimes (when scaling for the number of
oscillations). In an absolute sense, Matsumura et al. (2014) showed that for a con-
stant , the rise speed of a system scales with the square root of the gravity of
the system. In this study, we nd that for low friction values, the magnitude of
the rise speeds is consistent with this result, and that across 4 orders of magnitude
in gravity, a variation of exactly 2 orders of magnitude in the oscillation frequency
results in rise speeds that are within an order of magnitude of each other (in terms
of cm/osc). Moreover, again at low friction values, comparable to those used in
Matsumura et al. (2014), the rise speeds for the walls cases are very similar to that
found in their 104g case. However, this picture becomes challenged when varying
the friction parameters toward higher values.
In 1 g walls cases, we found that increasing s from 0.8 to 1.3 leads to a factor
of  6 increase in the rise speed (for low r). In the 104g cases, a comparable
variation in s leads to an increase in the rise speed by a factor of  20. For highly
63
frictional particles, the eciency of granular convection engendered by the presence
of walls is enhanced in low-gravity. We speculate that this may be due to dierences
in the packing of the granular bed (see Section 2.3.4) in low gravity. With a smaller
overburden pressure on the grains at the bottom of the bed due to gravity, the
upward ratcheting eect of granular convection may proceed much faster. While
we can only speculate on an enhancement mechanism at this point, it is interesting
to note that the data clearly show that the behavior of the grains diers in a low-
gravity environment. By varying the material properties of the grain, we can begin
to probe the limits of possible behavior of a granular system in micro-gravity.
To further mimic an actual asteroid environment, we deemed it necessary to
test the eect of varying the frequency and amplitude of the oscillation to levels
that might correspond to seismic shaking on a real asteroid. Several authors (e.g.,
Cheng et al. 2002 and references therein) claim that the seismic shaking that occurs
on an asteroid may have frequencies of the order of a few Hertz for millimeter-size
amplitudes. In our case, the amplitude of the oscillation is 1 cm, which may still be
an upper limit corresponding to more extreme shaking situations, but may not be
representative of the vibrations induced by micro-meteorite impacts. Garcia et al.
(2015) even nd seismic wave amplitudes of the order of a few microns. Taking
this information into account, we tested the eect of a smaller amplitude (always
keeping the same value of ) and ran simulations with an amplitude of 1 mm and
an oscillation frequency of 2.96 rad s1. This was done for a unique initialization of
the granular bed. The 1 mm value was preferred to a smaller one, because a lower
amplitude would have certainly resulted in poor results, as it would have represented
64
less than 1/10th of a grain diameter.
Our initial expectation was that none or very few of the low-amplitude cases
would actually exhibit the BNE, as previous studies (Matsumura et al., 2014) un-
derlined the sensitivity of the BNE to the amplitude and frequency of the oscillation.
According to their study, even amplitudes of half the grain size could lead to the
disappearance of the BNE, even at high oscillation frequencies. For these smaller-
amplitude cases, approximately 60% of the walls cases resulted in a full rise of the
intruder, while only 6% resulted in a BNE for PBC. Both boundary conditions
showed a strong dependence of BNE occurrence on friction properties. For the Wall
cases, only cases with s > 0.9 exhibited the BNE. For PBC, the few cases (5 in
total) that did exhibit BNE had high static and rolling friction (s > 1.0 and r >
1.6). In Fig. 2.8 we show the average rise speeds as a function of s (left) and r
(right), for the three dierent types of simulation in 104g: PBC cases with 1-cm
amplitude (solid line), walls cases with 1-cm amplitude (dashed line) and walls cases
at 1-mm amplitude (dot-dashed line). Due to a lack of rising Brazil nut obtained
with some PBC cases with 1-mm amplitude cases, the corresponding average rise
speeds were not included. We found that the number of oscillations required for the
intruder to rise to the top is much greater for the 1-mm amplitude than for the 1-cm
amplitude cases (about an order of magnitude greater). This was expected, as the
particle displacement at each oscillation is much smaller. While there is no strong
correlation of the average rise speeds with r, there is a rapid increase in the rise
time with s for cases that exhibit the BNE. This is similar to the eect seen for
the 1-cm amplitude cases.
65
2.5 Discussion
2.5.1 Driving mechanisms
In light of the results presented in Section 5.4, we investigate the mechanisms
that might be responsible for generating the BNE in a medium that is not conned
by lateral walls. To analyze the dierences in the way the BNE is driven, we consider
the evolution of the height of the Brazil nut for both PBC and walls cases, for the
same range of friction parameters as before (s = 0:8 to 1:3 in increments of 0:1 and
r = 0:2 to 2:0 in increments of 0:2).
The evolution in time of the height of the Brazil nut (solid curve) is plotted in
Fig. 2.9, for one set of friction parameters (s = 0:8 and r = 0:4). We compare this
evolution to the changes in the standard deviation of the entire bed's height (dotted
curve). This latter parameter shows oscillations driven by the periodic motion of
the bottom plate. More interestingly, large spikes in the standard deviation indicate
an event where the system signicantly dilates, creating void spaces in the medium.
If we compare the top panel of Fig. 2.9 (PBC case) with the bottom panel (walls
case), we can see that the BNE takes place in both cases. However, while the walls
case shows a gradual increase in the height of the Brazil nut, the PBC case shows a
more sporadic growth with some isolated jumps in the Brazil nut's height that are
correlated with spikes in the standard deviation of the bed's height.
The dierence in the evolution of the height of the Brazil nut in the two
cases visible here is evident in most of the simulations. It suggests a systematic
66
10
5
0
5
10
P
os
it
io
n
 (
cm
)
PBC
0 5 10 15 20
Time (s)
10
5
0
5
P
os
it
io
n
 (
cm
)
Walls
Figure 2.9: The height evolution of the Brazil nut (solid curve) diers
for the PBC case (top) compared to the walls case (bottom). The dotted
line represents the standard deviation in the height of the grains in the
system. For PBC cases, the Brazil nut experiences jumps in its height,
which are directly associated with the spikes in the standard deviation
of grain heights. These spikes in the standard deviation represent ex-
pansions in the granular medium and the creation of voids. These two
cases are for systems where the particles have friction parameters of s
= 0.8 and r = 0.4.
67
dierence in the mechanism driving the BNE depending on the presence of walls,
that was hinted at in the previous section (e.g., see Section 2.4.1). For the walls cases,
the steady increase in the height of the BNE suggests that a classical convection
mechanism takes place, as shown by previous authors (e.g., Knight et al. 1993). The
convection mechanism is driven by particle-wall interactions which force a downward
stream of particles along the walls, while interior particles are driven upwards. The
size of this downward stream is limited by the size of the small grains, therefore the
Brazil nut is not able to entrain itself into the downward 
ow and remains at the
top (see right panel of Fig. 2.5 for an example of this downward stream). For the
PBC cases, the sporadic height evolution of the Brazil nut might suggest a possible
random walk of the particle to the top. However, comparing the rise speeds of
the PBC cases to the walls cases in 1g, we nd that they rise at similar speeds (see
Fig. 2.3). This does not support the idea of a random walk, which would likely result
in longer rise speeds than that found in convection-driven walls cases. Since large
\jumps" in the height of the Brazil nut are correlated with spikes in the standard
deviation of the grains' height, we suggest that with PBC a void-lling mechanism
(Poschel & Herrmann, 1995) drives the BNE. When a vibration lofts the material
such that void space is created underneath the Brazil nut, the surrounding particles
have a larger probability of lling the void space than the Brazil nut due to their
smaller size. This geometric eect is able to drive the BNE, as subsequent vibration-
driven lofting events further cause small particles to lodge themselves underneath
the Brazil nut. In order to show that this mechanism is indeed responsible for
causing the BNE in PBC cases, we analyze the standard deviation of the grains
68
height for all simulations.
For each simulation where a full BNE cycle occurs (a Brazil nut is able to
move to the top and never sinks), we use a conventional peak-nding algorithm to
pick out each spike in the standard deviation of the bed's height and measure the
corresponding change in the height of the Brazil nut due to this void-lling event,

zvf . In each panel of Fig. 2.10, we show the sum of these changes in height due to
void-lling events, 
zvf , normalized by the bed depth htotal (each panel represents
a unique initialization of the granular bed). PBC cases are represented on the
left of each panel and walls cases on the right. We see that the PBC cases have
systematically higher values of 
zvf , indicating that the \jumps" of the entire bed
are much more signicant and frequent than with walls. This supports our claim
that the BNE, under Earth gravity, is mostly driven by void-lling events with
PBC, while with walls this mechanism has only some slight contribution. Since the
outcome of each simulation can be highly in
uenced by the starting conditions, we
measure the 
zvf in each simulation for four dierent initializations to show that the
dierence between the two cases persists regardless of initial conditions (Fig. 2.10
top-right, bottom-left and -right panels). The dierence in the mean of 
zvf for
the two populations (represented on the plot by a light dot-dashed line) is shown in
the top right corner of each panel.
In the top panels of Fig. 2.11, we compare the mean of the standard deviation
spikes (z;vf ), normalized by the grain radius, for each simulation (one initializa-
tion). This gure shows that void-lling events in PBC cases are driven by large

uctuations in the system height that are present in PBC cases. For walls cases,
69
0.0
0.2
0.4
0.6
0.8
1.0
1.2
1.4
Σ
∆
z v
f/
h
to
ta
l
∆m¯=0.37 ∆m¯=0.33
PBC WALLS
0.0
0.2
0.4
0.6
0.8
1.0
1.2
1.4
Σ
∆
z v
f/
h
to
ta
l
∆m¯=0.43
PBC WALLS
∆m¯=0.30
Figure 2.10: Sum of the changes in Brazil nut height due to void-lling
events for 4 dierent initializations of the granular bed. Each point
represents a single simulation. Each simulation results in a BNE (the
intruder rises to the top of the granular system). For each simulation
the change in the intruder height due to a void-creating event is tracked.
The thick vertical dashed lines separate simulations with PBC and simu-
lations with wall boundaries. The dot-dashed lines represent these mean
values of 
zvf for each boundary condition, and 
 m is the dierence
in the mean values. The values here represent the total net change in
the intruder's height over the course of the entire simulation. Points for
PBC cases cluster at high values, while points for walls cases cluster at
low values. Evidently, two dierent mechanisms drive the BNE for these
two dierent setups.
70
these 
uctuations are limited due to particle-boundary interactions, which restrain
the amount of possible system dilation.
We performed the same analysis for low-gravity simulations (amplitude of 1
cm and frequency of 0.939 rad s1), in order to assess whether or not the void-
lling mechanism is also responsible for the BNE eect in PBC cases at low gravity.
The ratio z;vf over grain radius obtained under these conditions for one suite of
simulations is shown in the bottom panels of Fig. 2.11. We nd a similar cluster-
ing of standard deviation spikes at high values in the PBC cases compared to the
Wall cases, which supports the idea that void-lling remains, under low gravity, the
principal mechanism driving the BNE when PBC are used.
The lower incidence of BNE in PBC cases at low gravity and low friction values
(16% of the cases showed a brazil-nut incapable of rising) might be due to an atten-
uation in the eciency of the void-lling mechanism compared to convection. We
speculate that this might come about from a diminution of the energy transmitted
to the granular bed by the bottom plate. In fact, as we keep the same value of  and
the same amplitude as in the 1g simulations for the low-gravity ones, we introduce
a factor of 104 in the kinetic energy transmitted to the grains compared to the
1g cases. With an energy four orders of magnitude lower than the one involved in
Earth gravity, we can expect smaller dilation spikes, thus smaller voids created un-
der the intruder at low friction values, and lower eciency of the void-lling process
compared to 1g. However, Fig. 2.11 shows that the mechanism is still active in low
gravity.
The void-lling mechanism has been previously studied with PBC by Schroter et al.
71
(2006). The authors showed that a system bounded by walls can have a BNE
driven by convection, while a PBC case is likely driven by void-lling events. How-
ever, these authors speculated that the vibrations that induce void-lling events are
much smaller than those that drive convection. Here, we show that for the same vi-
brational energy, two dierent mechanisms arise due to dierences in the boundary
conditions.
The results presented here marks the rst attempt at showing that a void-lling
mechanism may be responsible for size-segregation of large grains on an asteroid.
By using boundary and gravity conditions that more closely resemble an asteroid
environment, we suggest that vibrations are able to drive a boulder to the surface
of an asteroid. Furthermore, the void-lling mechanism is also enhanced by a larger
size ratio between particles, and, unlike the convection mechanism seen with walls,
it does not depend on the bed depth (Schroter et al., 2006). Previous studies (e.g.,
(Taguchi, 1992)) have shown that convection-driven mechanisms are limited to the
upper 10 or so grain heights of a system. After a certain depth, convection can no
longer eciently drive BNE. Since the void-lling mechanism occurs independently
of bed depth, a boulder entrained deep below the surface may still be able to segre-
gate upwards. Finally, the void-lling mechanism is characterized by independent
lofting events and does not necessarily require a continuity in the seismic shaking.
Thus it can be driven by independent shocks undergone by an asteroid over its life-
time. All of these facts make the void-lling mechanism particularly relevant to the
asteroid environment.
72
2.5.2 Interplay between gravity and rolling friction
Unlike previous studies, we have attempted to clarify the in
uence of material
property on the rise speed of the Brazil nut for dierent boundary conditions and
dierent gravity regimes. As Section 2.4.3 shows, a low-gravity system changes the
dynamical eects of rolling friction. In Earth gravity, PBC cases showed a negative
correlation with r; meanwhile, a positive correlation was found in low gravity.
Furthermore, an inspection of Fig. 2.6 shows that at intermediate gravity, 102 g,
there seems to be no correlation with r at all. These clues hint at a variation in
the granular 
ow of this system across dierent gravity regimes.
The PBC simulations in Fig. 2.11 are uniquely identied based on the static
friction coecient of the grains. On each panel, grains are grouped by increasing
values of s (from left to right), and each of these s groupings are sorted in order
of increasing r (again from left to right), when the data are available. We plotted,
for each s group, the regression line that best ts the dependence of z;vf with r
for that particular group. We nd that for Earth-gravity conditions, z;vf clearly
decreases with r for every value of s. At 10
4g, the trend reverses, and z;vf
increases with r for every value of s.
These trends suggest that, in Earth-gravity conditions, rolling friction dimin-
ishes dilation events in the system. Grains have a more dicult time 
owing past one
another as the system locks up. Even if voids are created, the grains are incapable
of 
owing into these voids. In low gravity, rolling friction enhances the strength of
dilation events in the system. The larger average value of z;vf (dot-dashed hori-
73
0.0
0.5
1.0
1.5
σ¯
z,
v
f/
r B
N
∆m¯=0.431 g
PBC WALLS
0.0
0.5
1.0
1.5
σ¯
z,
v
f/
r B
N
∆m¯=0.4910−4 g
µs = 0.8
µs = 0.9
µs = 1.0
µs = 1.1
µs = 1.2
µs = 1.3
Figure 2.11: Mean of the grain height standard deviation during void-
lling events. The thick vertical dashed lines separate simulations with
PBC and simulations with wall boundaries. The horizontal dot-dashed
lines represent the mean values of z;vf for each boundary condition, and

 m is the dierence in mean values. The granular systems in these
two setups behave dierently in response to the periodic forcing. For
PBC cases, the grains are more free to dilate and create voids. For walls
cases, the walls inhibit the grains from expanding too much. Hence, the
standard deviation in the height of the grains is systematically larger
for PBC cases than it is for walls cases. We plot the PBC cases with
unique scatter points based on the static friction value, s, used in each
simulation. Simulations with s of 0.8, 0.9, 1.0, 1.1, 1.2, and 1.3 are
represented by circles, x-marks, crosses, squares, left-facing triangles,
and right-facing triangles, respectively. Each of the 6 groups of s values
are sorted in order of increasing rolling friction, r. The dotted lines for
each s group represent regression lines establishing the dependency of
z;vf with r for each unique s.
74
zontal lines in Fig. 2.11) for PBC cases suggests that the system is lofted higher
and for a longer time in a lower-gravity environment. Systems with high-friction
grains will be able to create larger void spaces, as the system is able to dilate due
to stronger contact forces. This stronger force network may further feed into longer
lofting times as subsequent impulses from the bottom plate are more eciently
transmitted upwards. Thus, rather than diminishing the BNE, high rolling friction,
combined with longer lofting times, may create large enough void spaces that the
void-lling mechanism is enhanced. Therefore, while these rougher grains typically
inhibit granular 
ow, they can create larger void spaces that allow grains to move
more freely in low-gravity.
While a full physical model of this process is beyond the scope of this study,
we can begin to see its outlines. Namely, there are two important dynamical time-
scales involved in this process, whose predominancy seems to reverse with gravity;
namely, the rst being a gravitational timescale, and the second being a granular

ow timescale. The gravitational timescale is well known, being simply determined
by the size of the system and the local gravitational acceleration. The granular 
ow
timescale is more elusive to determine, as there are combined non-linearities that
come into play in a granular system. However, it is clear that what will chie
y
determine the eciency of granular 
ow, and the possibility of BNE in a semi-
boundless system, is the incidence of voids in that system.
75
2.6 Conclusion
We performed numerical simulations of vibration-driven grains. We attempted
to advance the current understanding of size segregation among asteroid regolith
grains in low gravity by performing simulations with realistic grain physical proper-
ties, boundary conditions, and gravity environment. Previous studies that used sim-
ulations of the BNE in the context of asteroids ran simulations with walls. Here, we
wanted to better mimic the asteroid environment by removing these lateral bound-
aries, which are probably not the most appropriate representation of an asteroid
interior. Instead, we implemented periodic boundary conditions (PBC). We were
able to show that for a wide range in friction properties, the BNE occurs even when
PBC are used. We also explored a broad range of friction levels, from medium to
high, which led us to the conclusion that a dierent mechanism drives the BNE
in walls compared with PBC cases. While granular convection is well known to
drive the BNE with walls, we have shown that a void-lling mechanism is actually
predominant when PBC are turned on. The void-lling mechanism occurs when a
vibration causes the granular system to dilate, creating void spaces in the medium.
In Earth gravity, a suciently 
uid granular system allows the creation of void
spaces that aords an opportunity for large-scale displacement of particles. Due to
the size dierence between the large intruder and the surrounding grains, the smaller
grains have a higher probability of lling the void space created below the intruder.
This allows the BNE to occur in a system without conning wall boundaries, which
are required for granular convection to occur.
76
For a better representation of the asteroid environment, we also studied the
BNE in low gravity (104g). We found that both convection for walls cases and the
void-lling mechanism for PBC cases are still able to occur, however with striking
dierences compared to the 1g simulations. For walls cases, we found that convec-
tion is counterintuitively enhanced for higher friction and speculate that the lower
packing in low-gravity cases may be responsible for this behavior. For PBC cases,
we found that the void-lling mechanism is still able to occur, and that the total
number of forcing cycles in the lower-gravity environment is comparable to that in
Earth gravity. However, our exploration of a broad range of rolling friction levels
highlighted major dierences in the eect of this rolling friction on the BNE. Unlike
for the 1g cases, where a higher rolling friction freezes the void-lling mechanism
(the grains are no longer able to ll the voids), low-gravity cases show a reverse
trend, where a higher rolling friction level actually enhances the BNE. We suggest
that this dierence comes from a change in the predominant time scale for the sys-
tem. In 1g, the gravity does not allow long lofting time and therefore, the higher
the rolling friction, the greater the inhibition of the grains' reorganization, which
occurs on a dierent time scale, the granular 
ow time scale. Conversely, in low
gravity, the granular 
ow time scale becomes more important as the gravity allows
much longer lofting time, enabling frictional particles to ll the voids anyway. In
addition, a higher rolling friction results in the creation of more voids to be lled
by the particles, which results in the inversion of the trend of the rise speed with
r for Earth- and low-gravity cases. We propose that the transition between the
predominancy of one regime or another may occur around 102g, as the in
uence
77
of r on the rise speed is for this case almost non-existent.
Finally, we also studied small-amplitude vibrations, where the amplitude is
much smaller than the size of the grains. We found that the BNE is less likely to
occur. This is especially true for cases with PBC. The PBC cases require suciently
large void spaces to be created during the vibration in order to drive the BNE. If
the shaking amplitude is much smaller than a grain size, then the subsequent voids
are likely to be small as well. Hence, the grains are unable to ll in voids created
underneath the intruder. Despite this, we nd that a few PBC cases are able to
exhibit the BNE. These are cases where both the static and rolling friction are high.
It is possible that in these cases the friction between grains is sucient to allow
voids to slowly grow between subsequent shakes until they become large enough for
a small particle to occupy them. Furthermore, a real asteroid will have grains that
span a wider size distribution. It is possible that grains with a size smaller than
the oscillation amplitude might be able to self-segregate eciently, in
uencing the
segregation process at larger scales.
The void-lling mechanism is a promising mechanism to explain features on
actual asteroids. Images and radar observations collected during space missions
indicate the presence of boulders at the surface of asteroids. In fact, a boulder is the
target of the Asteroid Retrieval Mission (ARM) under study at NASA, which aims
at collecting a boulder from the surface of the asteroid 2008 EV5 and placing it on a
cis-lunar orbit for later rendezvous by astronauts sent to interact with it (Abell et al.,
2015). The physical properties of boulders present on asteroid surfaces are intimately
related to their origin and history. Determining the mechanism bringing them to
78
the surface (BNE or any other mechanism) can provide valuable clues to determine
their physical properties. Future work could build on the simulations performed
in the present study and explore a larger parameter space, always with the aim of
characterizing asteroid interiors and environments more realistically.
79
Chapter 3: Low-speed granular impacts:
Modeling Spacecraft-Regolith Interactions
3.1 Chapter Preface
Some of the work presented in this chapter is in preparation as \Calibrating
Granular Mechanics Codes for a Range of Granular Materials at Very Low Gravity"
by R.-L. Ballouz, D. P. Sanchez, K. J. Walsh, H. C. Connolly Jr., D. S. Lauretta, D.
C. Richardson, and D. J. Scheeres, for publication in Icarus. This chapter presents
a suite of simulations that were performed for the upcoming NASA OSIRIS-REx
mission to the asteroid Bennu. As the work has progressed, it has been presented
in three science team meetings and one American Astronomical Society Division of
Planetary Sciences meeting, and has been written up as quarterly internal reports
for the OSIRIS-REx team. The goal of this work is to determine how the granular
surface of an asteroid might respond to the intrusion of a spacecraft. By conducting
low-speed impact experiments on to a variety of granular material types, we are
able to determine the range of possible regolith behavior that a spacecraft might
encounter as it touches down on an asteroid.
80
3.2 Background
NASA's third New Frontiers mission, OSIRIS-REx, will rendezvous in 2018
with a primitive near-Earth asteroid, (101955) Bennu. After being inserted into
orbit around the body and conducting measurements, which includes the creation of
a global surface map, a Touch-And-Go-Sample Acquisition Mechanism (TAGSAM)
will be deployed to collect material from the asteroid's surface (Berry et al., 2016).
The TAGSAM consists of a sampler head with an articulated pogo-stick-like arm.
The sampling strategy is to force nitrogen gas into the regolith thereby entraining
some material into the sampling device. However, the interaction of a sampling
mechanism with a regolith surface in the low gravity environment of a small body is
dicult to predict, even if the physical conditions on the surface of an asteroid are
well known. This is because the response of granular material to low-speed impacts
is poorly understood.
Hitherto, most studies of impacts onto asteroid surfaces have focused on de-
veloping scaling relations for predicting the outcome of cratering processes. Lab
experiments typically study impacts onto cohesive and non-cohesive grains by im-
pactors traveling at speeds of the order of 100 to 1000 m/s (Schmidt & Housen,
1987; Housen & Holsapple, 2011). In such studies, the fate of the impactor is
largely ignored. In the context of planetary science, the study of low-speed im-
pact collisions (impacts at speeds of  1 m/s or less) in low-gravity conditions were
performed to investigate collisions between planetary ring particles (e.g., Colwell
2003). However, there is a rich literature in the granular physics community of
81
studies of granular impact cratering in Earth-gravity conditions (e.g., Uehara et al.
2003; Katsuragi & Durian 2007). Uehara et al. (2003) attempted to develop scaling
of the crater depth as a function of impactor density, impact speed, and the material
type of the target. They suggested that the crater depth, d, formed by a low-velocity
impact onto granular material should scale as
d
a
=
0:175
tan()


1=2ga
v2
1=3
; (3.1)
where a is the impactor radius,  is the granular target's angle of friction,  is
the impactor density,  is the target density, g is the gravitational acceleration,
and v is the impact speed. While Eq. 3.1 does imply a scaling with gravity, there
were no experiments performed to verify if this were true. Subsequent work by
Katsuragi & Durian (2007) extended this work by conducting impact experiments
of 1-inch steel spheres onto 250-350 m glass beads. Katsuragi & Durian (2007)
proposed a unied force law for granular impact cratering, where the force on a
projectile by the granular target is composed of the sum of a velocity-dependent
inertial drag and a depth-dependent friction term. While this study did develop a
theoretical model to describe the complex response of granular materials to impact,
it did not include the eect of the local gravity on the outcome.
Driven by spacecraft missions, such as OSIRIS-REx and Hayabusa-2, that
will interact with the asteroid surface, the study of low-speed impacts on to re-
golith in low-gravity environments has only recently emerged (e.g., Nakamura et al.
2013). Nakamura et al. (2013) built upon the work of Katsuragi & Durian (2007)
82
by conducting laboratory experiments and numerical simulations of 6 mm plastic
projectiles impacting 50 m-size glass beads. They studied the deceleration pro-
les of the impactors and isolated the eects of three dierent components to the
deceleration of a projectile: a depth-dependent lithostatic pressure term, a velocity-
dependent viscous-like drag term, and a velocity-squared-dependent inertial-drag
term. Of these three components, only the rst, the lithostatic pressure term, has
a gravity dependence. By measuring the contributions of each term to the decel-
eration of the impactor, Nakamura et al. (2013) found that the gravity-dependent
lithostatic-pressure term has a very marginal eect on the deceleration compared to
the other eects. The lithostatic-pressure term was measured to have a magnitude
that was 1/40th of the inertial-drag term. This suggested that the local gravity has
a negligible in
uence on the outcome of a cratering impact. Altshuler et al. (2014)
conducted impact experiments at gravities ranging from 0.4g to 1.2g, and found a
similar lack of dependence of the penetration depth of an impactor with the local
gravity. While these studies suggest little to no gravity dependence on the outcome
of a granular impact, the range of gravities were limited to values of no more than
 0.1g ( 104  the gravity on Itokawa). It is unclear whether this description of
granular impacts is still valid at very low gravities.
Spacecraft missions that have interacted with the surface of a small solar sys-
tem body have resulted in rather unpredictable outcomes. These range from soft
landings (NEAR-Shoemaker at Eros, e.g., Dunham et al. 2002 to meter-high bounc-
ing (Hayabusa at Itokawa, see Yano et al. 2006) to spacecraft displacements of 100's
of meters (Philae at Comet 67P/Churyumov-Gerasimenko, see Biele et al. 2015).
83
The dierent outcomes of these touchdown attempts are likely due to dierence in
the material properties of the surface (
uy vs. strong material) and the gravita-
tional environment (which ranges from milli-g's on Eros to micro-g's on Itokawa).
In collaboration with the OSIRIS-REx mission, we studied granular impacts in low-
gravity environments using numerical simulations. Firstly, we calibrated our code
by conducting simulations of granular impacts in 1g and compared the outcomes to
laboratory experiments by Holsapple (2016). Then we extended our simulations to
low-gravity environments, in order to compare our simulations with theoretical pre-
dictions. Finally, we carried out simulations of the intrusion of a realistic physical
model of the sampling device (TAGSAM) into a bed of cm-size spherical particles
to estimate the extent to which the granular material resists the load applied by the
device.
The exact nature of the surface of Bennu (OSIRIS-REx's target asteroid )
is unknown. For that reason, it is necessary to study impacts and penetration
experiments over a large range of possible conditions and materials. Our goal is to
identify the governing physical parameters of that surface, and to assess how they
determine the outcome. Many important questions, especially that of the role of the
near-zero gravity, cannot easily be answered with physical experiments. Physically
simulating such a wide range of possible interaction scenarios on Earth is dicult
and expensive. Therefore, numerical simulations are an ideal tool to explore possible
sampling scenarios, in order to better understand the dynamics of the sampler and
determine any possible diculties that the mission may encounter. By generating
an \atlas" of possible outcomes, we are able to provide the mission with a tool
84
to reconstruct the actual sampling attempt. This is accomplished by producing
synthetic inertial measurement unit (IMU) data that we will be able to compare with
actual spacecraft data. Furthermore, once the actual sampling attempt happens,
our atlas may be used to infer the physical properties of the surface.
The rest of the chapter proceeds as follows. First, we introduce the code
updates that were implemented in order to model the collisional interaction of a non-
spherical body with spherical particles. Afterwards, we compare granular impact
simulations with laboratory experiments that were performed to calibrate our code.
Then, we discuss the simulations that were performed to generate an atlas of possible
outcomes that will be used to reconstruct the TAGSAM sampling attempt and that
can possibly characterize the properties of the surface grains on Bennu. Finally, we
conclude this chapter by summarizing our study and discussing the applicability of
the methods developed here for future space missions that will either interact with
a small solar system body's surface or deploy a lander.
3.3 Methodology
The study presented here uses a soft-sphere discrete element method code
to solve collisions between spherical particles (see Chapter 1 for details). The same
integrator also solves collisions between a spherical particle and a geometric primitive
(termed a wall) that either provides a boundary condition (such as a container that
holds particles) or an object that can collide with a particle. Non-inertial walls
that serve as boundaries were introduced in code developments by Richardson et al.
85
(2011). Here, we describe the developments that were implemented that allow us
to combine geometric primitives (spheres, rectangles, cylinders, etc.) to construct
arbitrary shaped inertial objects that can collide with spherical particles.
3.3.1 Inertial walls
In pkdgrav, inertial walls are geometric primitives with a nite mass that can
react to impulses by spherical particles (but not other inertial walls). These walls
can be combined to form a single assembly of walls that collectively react to impulses
as a single rigid body. An inertial wall has the following key properties that allow
it to behave as a rigid body: a mass, a center of gravity, a spin vector, a moment of
inertia tensor, and an orientation matrix.
When an inertial wall overlaps with a particle, the reaction force is calculated
using the usual Hooke's law spring force that is used for particle-particle collisions.
The resulting accelerations are solved using the leap-frog integrator and are applied
to the center of gravity of the inertial wall, similar to the method used for particle
collisions. The resulting change to the spin vector and orientation, however, requires
the rigid-body Euler equations to be solved. The changes to the spin are found by
solving the following dierential equations,
I1 _!1  !2!3(I2  I3) = N1
I2 _!2  !3!1(I3  I1) = N2
I3 _!3  !1!2(I1  I2) = N3
(3.2)
where Ik are the principal moments of inertia of the body, !k are the spin components
86
in the body frame, and Nk are the net torque components in the body frame (for k
= 1,2,3). The principal moments of inertia are user supplied information based on
the internal mass distribution of the wall assembly. The net torque components are
calculated by summing up the individual torques from each particle-wall collision.
Each individual torque, N , is given by
N = Nn +Nt Nparticle; (3.3)
where Nnormal and Ntangential are the normal and tangential torques, respec-
tively, found by
Nn;t = rpw  Fn;t; (3.4)
where rpw is the moment arm, calculated as the vector from the wall assem-
bly's center of gravity to the particle's center, and Fn;t is the normal and tangential
force by the particle on the wall assembly. Fn and Fn are found by taking the neg-
ative of the normal and tangential force on the particle by the wall assembly, which
have already been calculated. Nparticle is the torque on the particle due to the static,
rolling, and tangential friction forces that the particle feels by interacting with the
wall. The changes to the orientation, as re
ected in the principal axes components,
are found by solving the following equations,
87
_^p1 = !3p^2  !2p^3
_^p2 = !1p^3  !3p^1
_^p3 = !2p^1  !1p^2
(3.5)
where p^k denote the principal axes. Both sets of equations (Eq. 3.2 and Eq. 3.5)
are solved using a 5th-order time-adaptive Runge-Kutta integrator during the up-
date of particle positions. We use a Runge-Kutta integrator instead of a leap-frog
integrator because these equations are not easily adapted to the leapfrog integrator.
3.3.1.1 TAGSAM specications
The head will contact the surface of the Bennu asteroid with an initial approach
speed of about 10 cm/s. After a certain interval of time (< 5 seconds) of positive
deceleration of the spacecraft, nitrogen gas will 
ow through an annular nozzle
around the lower part of the collector head. This will mobilize the surrounding
regolith, allowing it to 
ow into the TAGSAM and be retained inside it. It will be
stored there until it is returned to Earth for laboratory analysis.
Fig. 3.1 is a diagram of the TAGSAM. We attempt to replicate its geometry
as close as possible by combining 6 geometric primitives (3 cylinders and 3 disks, see
Fig. 3.2). The collector head is made out of an outer hollow cylinder with a diameter
of  30 cm and an inner hollow cylinder with a diameter of  15 cm. Interior to
the main cylindrical housing is a smaller cylinder that tapers downwards. The rst
disk has a hole and is used to cover the gap between the bottom of the cylinder, the
88
Figure 3.1: Diagram of the TAGSAM. The collector head has an outer
diameter of about 30 cm, and an inner diameter of about 15 cm.
second disk covers the bottom of the tapered cylinder, and the last disk is used to
completely cover the top of the two cylinders. The total mass of the TAGSAM is
20 kg. It is attached to a spacecraft (represented by a 1,280 kg sphere) by a pogo-
stick like arm (represented by a cylinder). The TAGSAM head can move with six
degrees of freedom as a single rigid body. The spacecraft and TAGSAM dynamics
are coupled by the arm, which contains a constant-force spring that is designed to
begin compressing when the TAGSAM head feels a force that exceeds 67 N. The
spring ensures that the force felt by the spacecraft does not exceed 67 N.
In order to produce realistic acceleration proles of the spacecraft-TAGSAM
assembly, we included a spring force model that modulated the system's dynamics
as it penetrated the regolith bed. Our spring force model took the following into
account:
89
Figure 3.2: The TAGSAM used in the simulations. A disk covers the
top, but is left out of this image.
 If the ground force, Fground, as felt by the approaching TAGSAM assembly, is
below the 67 N threshold, then both the TAGSAM and the spacecraft (S/C)
feel a force,
FS=C = FTAGSAM = Fground; (3.6)
where FS=C is the force on the spacecraft and FTAGSAM is the force on the
TAGSAM.
 If the ground force exceeds the 67 N threshold, then the spring is engaged.
The S/C feels a force equal to the spring force, and the TAGSAM feels a force
equal to the dierence between the ground and spring force, Fspring.
90
FS=C = Fspring (3.7)
FTAGSAM = Fground  Fspring: (3.8)
This force conguration is kept as long as the spring is compressed (i.e. the change
in the separation between the TAGSAM and the S/C with respect to their initial
separation is  0). The spring is designed to not compress more than 15 cm.
Furthermore the spring may not extend beyond its initial equilibrium length. We
treat these two boundary conditions as a perfectly inelastic collision between the S/C
and the TAGSAM. In the event one of these hard-stops is reached, the TAGSAM's
position is set to where it should be for a 15 cm or 0 cm separation, and the TAGSAM
and S/C velocity is determined based on conservation of momentum for a perfectly
inelastic collision:
vnew =
mS=CvS=C +mTAGSAMvTAGSAM
mS=C +mTAGSAM
; (3.9)
where mS=C is the mass of the spacecraft, vS=C is the speed of the spacecraft,
mTAGSAM is the mass of the TAGSAM, and vTAGSAM is the speed of the TAGSAM.
Since the mass of the S/C is much larger than that of the TAGSAM, this new veloc-
ity will typically match the S/C's old velocity. This choice in a hard stop interaction
would cause a slow down for a downward moving S/C for the case of a 15 cm hard
stop (the TAGSAM is pushing up against the S/C) and a slow down for an upward
moving S/C for the case of a 0 cm hard stop (the TAGSAM is pulling down on the
S/C). This makes physically intuitive sense for such a conguration of a massive
object attached to a much smaller object by a rigid spring.
91
Finally, rather than being constant, in actuality, the spring force increases
the more the spring is compressed. Furthermore, the spring force discontinuously
changes depending on whether it is compressing or expanding. We model the spring
force as,
Fspring = k
z + Fspring;0; (3.10)
where k is a spring constant that changes depending on whether the spring is com-
pressing or expanding, 
z is the change in spring length (or change in separation
between the TAGSAM and S/C), and Fspring;0 is the spring force at zero compression
or expansion. For a compressing spring, the force ranges from 67 N at 0 compression
to 81 N at 15 cm compression. Therefore we can write the spring force as:
Fspring = 14
15

z + 67 (3.11)
with 
z in cm. For an expanding spring, the force ranges from 47 N at 0 compression
to 49 N at 15 cm compression. Therefore we can write the spring force as:
Fspring =  2
15

z + 47 (3.12)
3.4 Impact Experiments
In order to use our code to test the response of a granular bed on the surface
of an asteroid to a spacecraft landing, we must rst ensure that our code can accu-
rately predict the dynamics of a granular medium in Earth-gravity. In this section,
we rst present two collections of simulations that were conducted to ensure that
our grains behave similarly to realistic materials. The rst part of the calibration
92
involves measuring the bulk macroscopic properties of a granular medium as a func-
tion of its microscopic properties (elastic and friction properties). In the second
part, we compare numerical simulations of impact experiments to laboratory results
by Holsapple (2016) for grains with similar macroscopic properties. After we have
demonstrated the ability of our code to accurately predict impact outcomes in Earth
gravity, we perform simulations in low-gravity environments to check whether our re-
sults are in line with theoretical predictions of gravity-dependence (Nakamura et al.,
2013; Altshuler et al., 2014). Following the presentation of the code calibration re-
sults, we describe the simulations that were conducted to model the TAGSAM's
interaction with the granular bed.
3.4.1 Material Property
Before comparing cratering simulations with lab experiments, we rst verify
that the macroscopic properties of a granular bed can indeed be accurately mea-
sured. The materials we are attempting to simulate are glass beads and dry sand,
whose macroscopic properties are listed in Table 3.1 (Holsapple, 2016). The rele-
vant properties of the material that we nd through simulations are the Young's
modulus and the Poisson ratio. The Young's modulus measures the stiness of the
granular assembly; the Poisson ratio measures the relative amount of expansion of
the granular material in the direction perpendicular to the direction of compression.
We also measure the angle of friction of the material, through the angle of repose
landslide measurements described in Chapter 2.
93
Table 3.1: The physical properties of real granular materials.
Material Dry Ottawa Sand Coarse Solid Glass Beads
Density (g/cm3) 1.5 1.5
Grain Radius (m) 600-850 600-850
Sound Speed (m/s) 200 140
Poisson's ratio 0.25 0.25
Young's modulus (MPa) 135 70
Angle of Friction 35 22
The Young's modulus and Poisson ratio can be measured through uni-axial
compression tests on a packed granular medium. We conduct these simulations
by rst lling a 4 cm  4 cm  2 cm box with 4000 mm-size particles. Fig. 3.3
shows the setup for these simulations. The right wall of the box moves inward,
compressing the granular assemblage. This results in stresses to all the walls of the
conning box. The pressure on each wall of the box is tracked throughout the entire
length of the simulation. The Young's modulus is measured by taking the ratio of
the pressure on the forcing wall to the fractional change in the size of the granular
system. The Poisson's ratio is typically measured as the ratio of the tangential
strain to the normal strain (fractional changes in the size of the granular assembly
due to compression). Since we use non-moving walls, we instead measure the ratio
of the tangential pressure to the normal pressure. Fig. 3.4 shows an example of the
stresses measured in one of our simulations.
In our simulations, the sound speed of the material is determined by the elastic
properties of the grains, set by the normal spring constant kn. Therefore, this
parameter is tuned so that the sound speed of the simulated grains match that of
94
!"#$%&'()*++(,(-"#$%&'()*++(
./(
!0*&'1&2*+(
Figure 3.3: The setup for a simulation to perform uni-axial compres-
sion tests on grains. The uni-axial compression tests allow us to mea-
sure changes in the macroscopic properties of grains (Young's modulus
and Poisson ratio), based on microscopic properties (friction and elastic
properties of individual grains).
95
Figure 3.4: The resultant pressures on the right compressing wall
(green) and a tangential wall (blue) are shown as a function of the frac-
tional change in the length of the system (strain).
the actual material.
We ran uni-axial compression tests for walls moving at 0.2 and 2 m/s, for
coecients of normal and tangential restitution of 0.5 and 0.9, for static, s, coe-
cients of friction ranging from 0.0 to 1.0, and for rolling, r, coecients of friction
of 0, 0.1, and 0.2. Fig. 3.5 shows the result of sweeps of the two friction parameters
that control the surface roughness of the spherical grains, for the case of a conning
wall moving at 0.2 m/s, and for restitution coecients of 0.9. We see that both
the Young's Modulus and Poisson ratio are mainly controlled by the static friction.
The rolling friction has little to no eect on the outcomes of these experiments. The
Young's moduli range from 50{200 MPa, and the Poisson ratios range from  0.2
to 0.8.
96
We performed a consistency check by increasing the compression rate to 2
m/s, and found no signicant dierence in the outcome. For smaller restitution
coecients (0.5) we nd that the Young's modulus is slightly dierent: the range of
possible values is smaller (100{200 MPa).
By matching values of Young's modulus and Poisson's ratio, We nd that dry
sand is best simulated with grains that have a s of 0.8 and a r of 0.2, while glass
beads are best simulated by grains that have a s of 0.4 and a r of 0.1. The
compression tests have shown that the only time we have realistic Poisson ratios (<
0.4), is for moderately high friction coecients (s > 0.3). This highlights the point
that in order to produce realistic outcomes in SSDEM simulations, realistic friction
parameter values must be used.
3.4.2 Penetration Depth Experiments
After tuning our collision parameters to values that can accurately model
realistic materials, we used this information to conduct simulations of impact ex-
periments. By comparing with laboratory results by Holsapple (2016), we can show
the level of accuracy our simulations are able to achieve in predicting the outcomes
of impact experiments. These experiments also serve as a data set that can verify
the accuracy of the code calibration described in the previous section.
The impact experiments are set up by lling a cylinder with a 10 cm radius
and 9.2 cm height (with a covered bottom) with  400,000 spherical particles. The
spheres have a normal distribution of radii with a mean of 1 mm, and a standard
97
Figure 3.5: These plots show the measured macroscopic properties of
grains under uni-axial stress for the case of a conning wall moving into
the material at 0.2 m/s, and for restitution coecients of 0.9. The top
plot shows the range in Young's moduli that we nd as a function of
friction parameters. The bottom plot shows the corresponding Poisson's
ratio for each combination of friction parameters. The static friction
coecient has a stronger in
uence on both macroscopic properties of
the grains.
98
deviation of 0.1 mm. The distribution is cut o at 1 standard deviation. This is
done to prevent the grains from crystallizing or packing too closely together. The
spring constant of each grain is set such that the typical propagation speed of an
impact matches that listed in Table 3.1 across the grain. This is consistent with
previous impact simulations that use a soft-sphere approach (Wada et al., 2006).
We conducted low-speed impact cratering simulations onto two dierent gran-
ular materials: glass beads and dry sand. For impact experiments onto glass, we use
coecients of restitution parameters derived from previous studies (Schwartz et al.,
2013): n and t  0.9. Based on the uni-axial compression tests and matching the
friction angle of 22, we use s = 0:4, and r = 0.1. The simulated glass beads have
a Young's modulus  150 MPa and a Poisson ratio of 0.31. For impact experiments
onto dry sand, we use coecients of restitution parameters similar to that of gravel
(Yu et al., 2014): n and t  0.55. Based on the uni-axial compression tests and
matching the friction angle of 35, we use s = 0:8, and r = 0.2. The simulated dry
sand has a Young's modulus  200 Mpa and a Poisson ratio of 0.23. These values
are slightly dierent than those reported by Holsapple (2016). Any divergences of
our simulations with the lab experiments may result in part from these dierences.
The impactor has a diameter of 2.54 cm. Its density is varied between that of steel,
aluminum, and nylon (7.6, 2.7, and 1.15 g/cm3, respectively).
Fig. 3.6 shows the outcome of our impact simulations onto dry sand. The
shape of each data point represents whether it is a simulation (x's) or lab (circles)
result. The colors map the impactor's material type (density). For each case (unique
impact speed and projectile material), we performed three separate simulations with
99
dierent initializations of the granular bed in order to obtain some measure of the
uncertainty of our predictions. The error bars plotted here represent 1 standard
deviation in our measurements for each unique case. For the lab experiments, the
data provided has 5 to 6 cases for each impactor density and impact speed. For each
case, the mean value of the measured penetration depths is plotted. The error bars
show the standard deviation of the data. Overall we nd that our results match the
lab experiments fairly well. The penetration depths that we measure trend in the
same manner as that of the lab experiments. For the nylon and aluminum cases
we nd very close matches. However, for the steel cases we nd that the simulated
projectiles penetrate deeper. In general, we demonstrate that we can reliably match
the penetration depths to within a few mm for these cases at 1g.
In Fig. 3.7 we compare simulation results of impacts onto glass beads to the
corresponding lab experiments. For almost all cases we see a good match between
the simulation outcome to the lab data. The general trends in the data are similar
to that of the lab experiments. As expected, higher impactor densities and higher
impact speeds result in larger penetration depths. The one data point that deviates
wildly from the experiments is the 4.38 m/s case for nylon.
3.4.3 Gravity Dependence
Since we have shown that we can match the experimental data in 1g, we pro-
ceeded to run a suite of simulations in low gravity for the case of an aluminum
impactor on dry sand. We performed simulations for a range of impacts speeds at
100
Figure 3.6: Comparison of simulations (x's) to laboratory experiments
(circles) of a spherical impactor dropped in to a bed of dry sand. The
black, blue, and red circles represent the penetration depth of steel,
aluminum, and nylon impactors, respectively. The black, blue, and red
x's are the corresponding outcomes of numerical simulations.
101
Figure 3.7: Comparison of simulations (x's) to laboratory experiments
(circles) of a spherical impactor dropped into a bed of glass beads. The
black, blue, and red circles represent the penetration depth of steel,
aluminum, and nylon impactors, respectively. The black, blue, and red
x's are the corresponding outcomes of numerical simulations.
102
Figure 3.8: Simulations of an aluminum impactor onto dry sand at
dierent gravities. For the low-gravity cases, the impactor reaches the
very bottom of the container for the highest impact speed.
1g, 103g, and 106g. Fig. 3.8 plots the nal penetration depth of the impactor
as a function of impact speed for the three gravity cases. There is a clear indi-
cation from our simulations that penetration depth increases as the local gravity
decreases. At higher penetration depths, the impactors start interacting with the
bottom of the cylindrical container and reach the very bottom. Nevertheless, the
low-speed impacts show a clear trend with gravity that contradicts recent results by
Nakamura et al. (2013) and Altshuler et al. (2014).
103
According to Nakamura et al. (2013), the deceleration of a projectile may be
described as
F = k0Sgz  k1v  1
2
CdSv
2; (3.13)
where F is the force on the impactor by the granular medium,  is the density of
the granular material, z is the depth of the impactor, Cd is the drag coecient,
S is the cross-sectional area of the impactor, k0 is a dimensionless parameter that
encodes the physical characteristics of the grains (porosity, grain size, and surface
roughness), and k1 is a constant that represents the viscous characteristics of the
grains. The three components of the force law represent a depth-dependent litho-
static pressure term (rst term), a speed-dependent viscous-like drag term (second
term), and a speed-squared-dependent inertial-drag term (third term). Of these
three components, only the rst, the lithostatic pressure term, has a gravity depen-
dence. By measuring the contributions of each term to the deceleration of the im-
pactor, Nakamura et al. (2013) found that the lithostatic pressure term contributes
only 1/40th of the pressure compared to the inertial drag term. This suggested
that there should be no gravitational dependence for the penetration depth, since
the lithostatic pressure has a negligible eect even at 1g. Therefore, at even lower
gravities, the dierence in the penetration depth should be minimal.
However, our simulations do show that there is indeed a dierence in the pen-
etration depth with gravity. We postulate that the likely cause of this phenomenon
is a change in the granular 
ow properties of grains in a low-gravity environment.
This is likely caused by a dependence of the drag term Cd or the viscous term k1 on
104
gravity.
To test this idea, we can measure the force on the projectile as it penetrates
through the granular medium. If the resistance to penetration is dominated by an
inertial drag eect, then the magnitude of the force on the impactor is described by
the third term on the right hand side of Eq. 3.13:
F =
1
2
CdSv
2: (3.14)
Fig. 3.9 shows the force on the aluminum impactor as a function of its speed for two
dierent impact speeds (1.4 and 4.58 m/s) and gravity strengths (1g and 103g). We
also plot Eq. 3.14 for 5 dierent cases of the drag coecient, Cd. We see that both
cases have similar slopes initially. However, after some deceleration, large dierences
begin to emerge. For both cases, the 1g case is better described by a drag force with
a higher drag coecient than the 103g case. The 1g cases are best-t by a drag
coecient Cd between 1.5 and 2.0, while the 10
3g cases are described by a drag force
with a drag coecient between 0.25 and 0.5. The data for 106g cases (not shown
here for clarity) are best t with Cd between 0.1 and 0.15. These results suggest that
the drag properties of a granular system depend on the local gravity. We propose
that the drag coecient varies as log(g). This functional form suggests that changes
in granular 
ow properties of regolith only begin to manifest in the very low-gravity
environments of small bodies. This also explains the results of previous authors
(Nakamura et al., 2013; Altshuler et al., 2014) who found virtually no dependence
of the penetration depth with gravity because their studies were limited to local
gravities of at most 101g.
105
100
v [m/s]
10-2
10-1
100
101
F
o
rc
e 
(N
)
Cd =2. 00
Cd =1. 50
Cd =1. 00
Cd =0. 50
Cd =0. 25vimp =  1. 40 m/s
1g
10−3g
100
v [m/s]
10-2
10-1
100
101
102
F
or
ce
 (
N
)
vimp =  4. 58 m/s
Figure 3.9: These plots show the force on an impactor as a function of
its instantaneous speed. The top panel shows the cases for an impact
speed, vimp, of 1.4 m/s, while the bottom panel shows the cases with a
vimp of 4.58 m/s. The solid black curves represent impactors at 1g, while
the dot-dashed curve represent impactors at 103g. The blue, green,
yellow, magenta, and red lines are Eq. 3.14 for a drag coecient, Cd, of
2.0, 1.5, 1.0, 0.5, and 0.25, respectively.
106
3.5 OSIRIS-REx TAGSAM
In this section, we describe simulations of the TAGSAM penetrating a bed of
granular material. The goal of this study was to construct an atlas of simulation
results that can be used to detail the range of possible outcomes of a sampling
attempt. This atlas of results has multiple uses. First, by comparing IMU data
from the actual sampling attempt to our simulation results, we will be able to
reconstruct what occurs to the S/C and TAGSAM during the sampling attempt.
Second, we will be able to characterize the regolith since the response of the S/C
and TAGSAM depends on the physical properties of the surface grains. Finally,
observations during the orbital survey of Bennu will measure some of the surface
properties of the grains, such as their mean size, size distribution, and mineralogy. If
these characteristics are well constrained, then the sampling attempt can also act as
a useful and rare impact experiment onto grains in a very-low-gravity setting. Based
on the S/C's reaction, we may be able to make a measurement of the granular 
ow
properties of grains in low gravity.
3.5.1 Simulations
Fig. 3.10 shows snapshots of the TAGSAM penetrating a bed of grains. Our
simulations are set up by rst lling a cylindrical container with a height of 60 cm and
a radius of 60 cm with 150,000 grains and allowing them to settle in low gravity (20
m/s2 or  2 g). Then, the actual touch-down simulation is conducted by placing
the TAGSAM  4 cm from the surface of the granular bed. The TAGSAM has an
107
Figure 3.10: A model TAGSAM (visualized with some slight trans-
parency) penetrates a bed of  150,000 particles (seen in cross-section)
with an initial downward vertical speed of 10 cm/s in the microgravity
environment of Bennu ( g = 20 m / s2). The cylindrical container has
a heigh and radius of 60 cm.
initial downward speed of 10 cm/s. The total simulation time is 5 s, and the time
step is 105 s. Each grain has a density of 2.7 g/cm3, typical of basalt. The regolith
is composed of a power-law distribution of spherical particles (power-law index -2.5)
with a minimum radius, Rmin, of 0.5 cm, and a maximum radius, Rmax, of 1.5 cm.
We explored the eect of variations in the friction of the grains. However, we kept
the normal and tangential coecients of restitution constant at 0.55. These were
the values adopted for dry sand in the uni-axial compression simulations described
earlier.
An important aspect in these simulations was modeling the eect of the spring
force between the TAGSAM and the S/C. The spring modulates the force that
the S/C feels from the ground. It only engages when the TAGSAM feels a force
greater than 67 N. The equations governing the eect of the spring were described
previously in Eq. 3.6 and Eq. 3.8. Fig. 3.11 shows how the spring force is tracked
(cyan curves), and its eect on the displacement between the TAGSAM and the
108
S/C (black curve), for a case with particles that have properties similar to that of
dry sand. The numbered regions indicate key events in the dynamical interaction
between the S/C, the spring, and the TAGSAM:
1. The TAGSAM feels a force from the ground that exceeds the 67 N threshold
and the spring engages and starts compressing (time  0.5 s).
2. As the spring is further compressed, the spring force ramps up (as described by
Eq. 3.11), until the displacement reaches a minimum, then begins increasing
(time  1.7 s).
3. Afterwards, the spring expands, and the spring force is described by Eq. 3.12
instead. This continues until we reach the 0 cm hard stop (time  3.3 s).
4. At this point the spring is no longer engaged and the ground force is below
the spring threshold. The spring force remains at zero until the end of the
5-second simulation.
In order to put some more context into the behavior of this dynamical system,
we analyzed the time evolution of the ground force. Fig. 3.12 shows the magnitude
of the force imparted by the grains on the TAGSAM head (red scatter points, right
axis). Furthermore, we show the vertical distance from the surface of the S/C
(purple curve) and TAGSAM (blue curve, left axis).
Overall, the force prole begins with a sharp peak and then proceeds in a
constant-force-step prole. The magnitude of the force at each step is roughly equal
to the spring force as determined by Eq. 3.11 and Eq. 3.12. The numbered regions
109
Figure 3.11: Plot of the displacement between the S/C and TAGSAM
(black curve) as a function of time. Changes in the relative positions
of the S/C and TAGSAM are driven by the spring that connects the
two objects. The spring force (cyan curve) engages after the TAGSAM
feels a ground force that exceeds 67 N. The force by the spring depends
on whether it is undergoing compression or expansion. This behavior is
described by Eq. 3.11 and Eq. 3.12. See text for details on numbered
key events.
110
Figure 3.12: Plot of the ground force on the TAGSAM (red scatter
points, right axis). On the same plot, we show the vertical distance from
the surface of the TAGSAM (blue curve, left axis) and the S/C (purple
curve, left axis). See text for details on numbered key events.
111
on the plot indicate key events:
1. At t  0.5 s, the TAGSAM makes contact with the ground and experiences
a sharp impulse of  120 N. The spring immediately engages as the force
threshold of 67 N is exceeded.
2. The spring engagement causes a force balance between the downward moving
TAGSAM and the ground. The force by the ground on the TAGSAM is
roughly at the spring threshold (67 N). The TAGSAM has a slight downward
speed of 2 cm/s. Due to the force balance this remains roughly constant.
3. Once the S/C stops moving downwards, the spring begins to expand. This
causes the spring force to drop to the 47{49 N range (Eq. 3.12). Since the
TAGSAM is forced downward by the spring, it continues to interact with the
ground in this force range. The ground pushes back at the TAGSAM with an
equal force. The TAGSAM remains in force balance.
4. Once the S/C begins to move up, the spring expands to the 0 cm hard stop.
At this point we have a perfectly inelastic collision between the TAGSAM and
S/C (Eq. 3.9). The S/C drags the TAGSAM upwards with it.
Fig. 3.11 and Fig. 3.12 show how the spring dictates the dynamical interaction
between the TAGSAM and the S/C, as well as modulating the total force that the
TAGSAM feels during the sampling attempt. In the following section, we discuss
how we model this interaction for dierent surface mechanical properties and how
we will be able to infer these properties from the spacecraft's response.
112
3.5.2 Penetration Depth
In order to quantify the eect that dierent materials might have on the sam-
pling attempt, we conducted simulations that explore a wide parameter space of
two material friction coecients (static and rolling). These parameters combine to
provide a proxy for material roughness (a parameter that may also be in
uenced by
grain size, size distribution, shape, etc.). Fig. 3.13 shows a map of the penetration
depth of the TAGSAM as a function of the rolling and static coecients of friction.
We see that the penetration depth of the TAGSAM depends sensitively on the
friction properties of the grains. When the static friction coecient is low (< 0.4) or
when there is no rolling friction, the granular bed oers almost no resistance and the
S/C - TAGSAM assembly is able to sink into the granular bed almost unimpeded.
A nal penetration depth of  50 cm would signify a freely falling spacecraft. If the
regolith surface were made of material with friction properties similar to glass beads
(s = 0.4, and r  0.1), the spacecraft would feel very little resistance, and sink
deeper than  20 cm after 5 s. If instead the material were like dry sand (s = 0.8,
and r  0.2), the spacecraft would feel a strong response from the granular bed,
and only penetrate 2{3 particle diameters into the granular bed ( 4 cm). While
we show that there can be a wide range of responses to the penetration, the results
clearly segregate into two general categories: almost unimpeded penetration (yellow
squares) and strong resistance (light and dark red squares). Realistic materials
will almost always have friction properties that would place them in one of these
red squares, whose lowest friction value corresponds to a material with an angle of
113
Figure 3.13: A map of the maximum penetration depth that the
TAGSAM reaches for each combination of static and rolling friction
parameters. Whiter/yellower colors signify more penetration, and red-
der/blacker colors signify less penetration.
114
friction of  22. Therefore, the outcome of the actual sampling attempt will likely
result in the TAGSAM penetrating into at most 20 cm of the regolith bed.
In order to understand why the penetration depth depends so sensitively on the
friction properties of the grains, we can draw comparisons with recent studies into
shear-thickening 
uids (Mukhopadhyay et al., 2014). Mukhopadhyay et al. (2014)
conducted studies on granular material (cornstarch) suspended in 
uids. Such 
uid
suspensions have the surprising property of resisting heavy loads when the impact
speed of the projectile is suciently high. As a projectile impacts a 
uid suspension,
a shock wave travels through the suspension that causes the entire medium to dilate,
leading to the interlocking of nearby grains.
It is possible that a granular material in low gravity can behave like grains
suspended in a 
uid. Once the sampler makes contact with the surface, the material
underneath the sampler shear-thickens. As the material underneath the sampler ex-
pands, the highly frictional grains form more particle contacts. This creates a tran-
sition where the grains move from 
uid-like 
ow to a solid-like response. Fig. 3.14
shows snap-shots of the end of two simulation. The top row shows a case with high
friction (s, r = 0.8), and the bottom row shows a case with low friction (s, r
= 0.2). The gures on the right show the change in the column density of grains
from the beginning of the simulation. The large swaths of dark blue show regions
where particles have been excavated. While the case with low friction shows regions
of over-densities (red), the case with high friction sees some slight under-densities.
This suggests that particles with high friction are able to shear-thicken by dilat-
ing. This results in more grain-grain contacts, and stronger force chain networks
115
Figure 3.14: The plots on the left show the dierence in the penetration
depth of the TAGSAM for a case with high friction (top) and low friction
(bottom) 5 seconds after initial contact. The plots on the right are cor-
responding maps of the column density through the granular bed. Bluer
colors correspond to under-densities with respect to the initial density
distribution of grains. Redder colors correspond to over-densities.
116
throughout the system, leading to an overall stronger material. For the other case,
insucient frictional contact causes the material to behave like a 
uid in micrograv-
ity. The sampler plunges through, creating transient gaps and clusters.
3.5.3 TAG Reconstuct
We have demonstrated that friction properties of the grains play a crucial
role in determining the nal penetration depth of the TAGSAM. However, for most
cases that mimic realistic grains, the range in penetration depths we nd is relatively
narrow. For grains with angles of friction similar to that of typical terrestrial rocks
or sand (>  30), we expect the penetration depth to be < 10 cm after 5 s.
Therefore, in order to dierentiate the outcomes for realistic materials, we analyze
three dierent indicators: the force prole on the TAGSAM-S/C assembly, the
change in the S/C speed, and the duration of spring compression. Using these
indicators might allow us to more accurately reconstruct the sampling attempt, and
determine the properties of the surface grains.
3.5.3.1 Force Proles
Since onboard accelerometers are limited to the S/C bus (no such devices
will be on the TAGSAM head), our only indication of the type of material we
will be interacting with (and the regolith's reaction to the sampling attempt) will
be the data from S/C accelerometers. Therefore, any dierence in the measured
deceleration due to dierent materials would only be detected before the spring
117
has engaged. We analyzed four dierent proles that span the range of expected
material properties, s = 0.4{1.0 and r = 0.4{1.0, to test whether any dierence
in the force proles can be detected before the spring is engaged. Fig. 3.15 shows
two cases with moderate friction. The left panel shows a case where both s and
r = 0.4. The right panel shows a case where both s and r = 0.6. Both of these
cases also exhibit a \stepped" force prole similar to that shown in Fig. 3.12. The
\stepped" regions of the force prole indicate times when the spring is engaged.
Therefore, the only valuable information that can be extracted is from the peak of
the force prole. These two cases do show some dierences in their overall force
prole. In the left panel, the \steps" are shorter. For this case, the spring doesn't
compress very much and extends back to its original length relatively quickly. Once
it does so, the TAGSAM and S/C dynamics are coupled again, and the reaction force
on this assembly becomes much lower due to the drastic decrease in their downward
speed. Thus, they are able to penetrate further into the granular bed. By the time
the simulation ends, the S/C is close to stopping its downwards trajectory. The
case with s; r = 0.6 is slightly dierent. The spring is compressed for a longer
time, allowing the S/C to slow down suciently to begin an upward trajectory that
eventually leads to a 0 cm hard stop that drags the TAGSAM upwards. Furthermore,
the initial impulse that the TAGSAM feels upon rst contact is higher. Therefore,
the TAGSAM is able to penetrate less in this case. Once the 0 cm hard stop is
reached and the S/C-TAGSAM assembly begins moving upwards, the TAGSAM is
embedded close to 10 cm into the regolith bed.
For even higher values of the friction coecient, Fig. 3.16 shows that the force
118
Figure 3.15: Force proles for two cases with dierent friction proper-
ties. The left panel shows a case where both s and r = 0.4. The right
panel shows a case where both s and r = 0.6. The ground force on
the TAGSAM (red scatter points, right axis), the vertical distance from
the surface of the TAGSAM (blue curve, left axis), and the S/C (purple
curve, left axis) are shown.
proles still have the \three-stepped" pattern. The left panel shows a case where
both s and r = 0.8. The right panel shows a case where both s and r = 1.0.
One key feature that dierentiates these two cases is the peak of the force prole.
The case with s; r = 1.0 has a peak force (at  130 N) that is about 30 N higher
than the case with s; r = 0.8 (at  100 N). The lowest friction case presented
here (s; r = 0.4) follows this trend, as it has a peak force of  90 N. However,
the case with s; r = 0.6 (which has a peak force of 120 N) does not follow the
general trend. Overall, a measurement of the peak force could be an indicator of
the frictional properties of the granular matter. However, a calibration study with
the actual instrument would be required if this technique were to be used for a
measurement of the grain properties.
119
Figure 3.16: Force proles for two cases with dierent friction proper-
ties. The left panel shows a case where both s and r = 0.8. The right
panel shows a case where both s and r = 1.0. The ground force on
the TAGSAM (red scatter points, right axis), the vertical distance from
the surface of the TAGSAM (blue curve, left axis), and the S/C (purple
curve, left axis) are shown.
3.5.3.2 S/C vertical speed
Since we will have no inertial measurements from the TAGSAM itself, a useful
indicator for dierent regolith dynamics may be the behavior of the S/C. Fig. 3.17
compares the time evolution of the S/C's vertical speed for the range of material
parameters discussed in the last section. Here, we clearly see that there is a marked
dierence in the S/C's deceleration for dierent material properties. The black and
purple lines represent the upper end of our friction spectrum where things begin to
behave in a similar manner (similar penetration depths and TAGSAM ground force
proles). At a simulation time of about 2 seconds (or about 1.5 seconds after initial
contact), the eect of dierent materials on the deceleration prole becomes very
noticeable (in order of increasing friction, the vertical speeds are: 3:7, 1:7, 1:2,
120
Figure 3.17: The vertical speeds of the spacecraft as a function of time.
The red, cyan, magenta, and black curves represent cases where the
grains have friction coecient s and r = 0.4, 0.6, 0.8, and 1.0, respec-
tively. For dierent material types, the S/C will decelerate dierently.
The dierence in deceleration becomes very noticeable  2 seconds after
contact.
and 0:99 cm/s). Therefore, by tracking the deceleration prole of the S/C, we may
be able to determine the properties of the regolith itself. Furthermore, since each
deceleration prole is unique, we would also be able to link this to the behavior of
the TAGSAM, eectively reconstructing the sampling attempt.
121
3.5.3.3 Spring Compression Time
Finally, another possible indicator that would allow us to reconstruct the sam-
pling attempt and determine regolith properties is the total amount of time the
spring remains in a compressed state. Fig. 3.18 is a map of the total time a spring
is compressed for each combination of s and r. By comparing the total time the
spring compression time for dierent material properties, we see that there is some
correlation with the rigidity of a material. This can be a useful diagnostic tool if
the amount of compression can be measured in real time.
3.6 Conclusions & Future Work
This chapter presented a study of low-speed impacts onto granular material.
While low-speed impacts have been studied by the granular science community
(e.g., Katsuragi & Durian 2007) for some time now, the study of low-speed granular
impacts on small solar system bodies has only recently emerged due to interest
generated from spacecraft missions to asteroids and comets. The study of granular
material in a low-gravity environment is still in its infancy. It is dicult to develop
a deterministic theory for granular materials in any environment due to the fact
that granular systems are composed of a large number of elements that interact
through a non-linear combination of various forces (mechanical, gravitational, and
electrostatic, for example) leading to a high degree of stochasticity. Therefore,
numerical simulations are important for placing constraints on the range of possible
behaviors that can be expected from a granular system.
122
Figure 3.18: A map of the total time the spring between the S/C and
TAGSAM is compressed. Redder colors signify more time, and bluer
colors signify less time.
123
The OSIRIS-REx mission to Bennu will provide the community with a very
unique low-speed impact experiment. We have demonstrated that the outcome of
a sampling attempt by the TAGSAM can vary depending on the properties of the
regolith. Granular material behaves more 
uid-like in a low-gravity environment,
and does not resist penetration as much as it would in 1g. This was also demon-
strated through low-speed impact experiments of an aluminum sphere onto dry
sand. We rst calibrated our simulation to the behavior of real granular material,
and showed that we can replicate the results of laboratory experiments. Extending
our simulations to low-gravity conditions, we showed that the penetration depth
has a dependence on the local gravity. However, unlike previous results, which sug-
gested that changes in penetration depth are likely due to the lack of lithostatic
pressure in low gravity, we showed that it is more likely that it is a change in the
granular 
ow properties of the material in low gravity that causes a change in the
nal penetration depth.
A key goal in this study was to demonstrate the ability of our simulation
software to predict the outcome of a spacecraft's interaction with a regolith surface.
We were able to create an \atlas" of simulation results that the future sampling
attempt can be compared against. This \atlas" will be able to reconstruct the
sampling attempt and possibly help determine the properties of the regolith. We
studied the usefulness of three dierent indicators that might allow us to better
utilize our atlas. The S/C deceleration prole and the spring compression time are
promising metrics that might allow us to better constrain the regolith properties.
Future work will include extending the results presented here for grains with
124
non-zero cohesion. We expect the dynamical in
uence of cohesion to be similar to
that of friction. Furthermore, we intend to use the software that was developed
for this task on future space missions to small-body systems. In particular, we can
simulate the dynamics of the MASCOT lander on Hayabusa-2 and the MASCOT-2
lander on the proposed Asteroid Impact Mission (AIM) under study at the European
Space Agency (ESA).
125
Chapter 4: Catastrophic Disruption of Gravitational Aggregates
4.1 Chapter Preface
The work presented in this chapter was published in the Astrophysical Journal
(ApJ) as Ballouz et al. (2014), and in the Planetary and Space Sciences Journal as
Ballouz et al. (2015). In this chapter, we carry out a systematic exploration of the
eect of pre-impact rotation on the outcomes of low-speed collisions between plan-
etesimals modeled as gravitational aggregates. A rotating body has lower eective
surface gravity than a non-rotating one and therefore might suer more mass loss as
the result of a collision. What is less well understood, however, is whether rotation
systematically increases mass loss on average regardless of the impact trajectory.
This has important implications for the eciency of planet formation via planetes-
imal growth, and also more generally for the determination of the impact energy
threshold for catastrophic disruption (leading to the largest remnant retaining 50%
of the original mass), as this has generally only been evaluated for non-spinning
bodies. We nd that for most collision scenarios, rotation lowers the threshold en-
ergy for catastrophic dispersal. For head-on collisions, we develop a semi-analytic
description of the change in the threshold description as a function of the target's
pre-impact rotation rate, and nd that these results are consistent with the \uni-
126
versal law" of catastrophic disruption developed by Leinhardt & Stewart (2009).
Using this approach, we introduce re-scaled catastrophic disruption variables that
take into account the interacting mass fraction of the target in order to translate
oblique impacts into equivalent head-on collisions. Furthermore, we analyze how
rubble pile of dierent material types catastrophically disrupt. We nd that ma-
terial with higher dissipative and friction properties requires larger impact energy
to catastrophically disrupt. However, we nd that the increase in mass loss due to
pre-impact rotation is constant across dierent material types.
4.2 Background
Much of the evolution of small Solar System bodies (SSSBs) is dominated by
collisions, whether from the initial build-up of planetesimals (Lissauer, 1993) or the
subsequent impacts between remnant bodies that exist today (Michel et al., 2004).
Outcomes of collisions between SSSBs are divided into two regimes: those domi-
nated by material strength and those dominated by self-gravity (Holsapple, 1994).
The transition from the strength to the gravity regime may occur at body sizes
as small as a few kilometers or less for basalt (Benz & Asphaug, 1999; Jutzi et al.,
2010). After their formation, planetesimals interacted with one another in a dy-
namically cold disk (Levison et al., 2010). This allowed planet-size objects to form
through collisonal growth.
Since the dominant source of conning pressure for planetesimal-size SSSBs is
self-gravity, rather than material strength, they can be assumed to be gravitational
127
aggregates (Richardson et al., 2002). Hence, the collisions can often be treated as
impacts between rubble piles, the outcomes of which are dictated by collisional dis-
sipation parameters and gravity (Leinhardt et al., 2000; Leinhardt & Richardson,
2002). Understanding the eects that contribute to changes in the mass (accretion
or erosion) of gravitational aggregates is important for collisional evolution models
of the early solar system (Leinhardt & Richardson, 2005; Weidenschilling, 2011).
The outcomes of impacts in these models are paramaterized through a catastrophic
disruption threshold Q?D (Benz & Asphaug, 1999), which is the impact energy per
unit mass (termed the specic impact energy) required to gravitationally disperse
half the total mass of the system, such that the largest remnant retains the other
half of the system mass. However, few studies have accounted for the eect of pre-
impact rotation on the size evolution of SSSBs.
A rotating body has lower eective surface gravity than a non-rotating one
(with the dierence being greatest for surface material at the equator and decreas-
ing for material closer to the rotation axis). Therefore, a rotating body might suer
more mass loss as the result of a collision. A recent laboratory study by Morris et al.
(2012) suggests this is true for solid bodies. What is less well understood, however,
is whether rotation systematically increases mass loss on average regardless of the
impact trajectory.
In order to explain the collisional evolution of rotation rates of asteroids,
Dobrovolskis & Burns (1984) evaluated analytically the sensitivity of mass loss to
rotation for cratering impacts on rigid bodies. They found that the angle-averaged
mass loss for cases with rotation is enhanced by factors of  10{40% compared to
128
cases without rotation for rotation speeds  40{80% of the critical spin rate (see
their Fig. 2). Analytic and numerical work by Cellino et al. (1990) showed that
catastrophic disruptions, rather than cratering events, were a bigger contributor to
the rotational evolution of asteroids through an angular momentum \splash" pro-
cess; however, they and subsequent authors Love & Ahrens (1996) focused on the
eects of spin-state evolution change rather than mass loss.
Other authors have included pre-impact rotation in their numerical simulations
of planetesimal and protoplanet collisions; however, except for Takeda & Ohtsuki
(2009) none have systematically studied its contribution to mass loss in the disper-
sive regime. Using a hard-sphere model, Leinhardt et al. (2000) performed numeri-
cal simulations of collisions of equal-size bodies with pre-impact rotation; however,
their work focused on the eect of rotation on the shape of the largest remenant.
Canup (2008) studied the eect of pre-impact rotation on lunar formation; however,
the work focused on a non-dispersive collision regime. Using a soft-sphere collision
code, Takeda & Ohtsuki (2009) performed simulations of hyper-velocity impacts on
rotating  10-km size bodies. They found that mass loss is only sensitive to ro-
tation when the target has an initial spin period close to break-up; otherwise, the
collisional energy needed to disrupt a rubble-pile object is not aected by initial
rotation. They argued that, upon collision, the ejection speeds of fragments in the
hemisphere rotating away from the projectile (prograde direction) are accelerated
by the initial rotation, but this is balanced by fragments in the hemisphere rotat-
ing toward the projectile (retrograde direction) being decelerated. However, their
analysis was restricted to targets with initial rotations of 2.6 and 4.6 revolutions
129
per day (9.23 and 5.58 hours, respectively). Furthermore, their work focused on the
eciency of angular momentum transfer in catastrophic collisions.
In this chapter, we expand upon the work of previous authors by performing
a systematic study of the eect of pre-impact rotation on the energy required to
disperse material from km-size gravitational aggregates rotating with spin periods
of 3, 4.5, and 6 hours. We solve numerically the outcomes of rubble-pile collisions
using our soft-sphere discrete element method (SSDEM) collisional code and the
numerical gravity solver, which is needed to accurately model the reaccumulation
stage. SSDEM has the numerical resolution to determine the mechanics involved in
enhancing or diminishing the amount of mass loss associated with collisions onto a
rotating target. SSDEM permits realistic modeling of multi-contact and frictional
forces between discrete indestructible particles. Thus, it is well suited to study low-
speed (a few to tens of m s1) impacts, as it can robustly model collisions that
do not produce irreversible shock damage to material (as in hypervelocity, km s1
impacts). In the quasi-steady-state collisional system generally present in a pro-
toplanetary disk, impact speeds are typically on the order of the escape speed of
the largest body in the vicinity. Until the largest body becomes protoplanet sized,
impacts will be typically at speeds less than the sound speed of the assumed rocky
material. Hence, we limit our study to collisions that occur at subsonic speeds.
Most signicant collisions today occur at supersonic speeds; however, studies of su-
personic collisions require the use of shock physics codes, which include the eects
of irreversible shock deformation.
Furthermore, we attempt to revise the dependence of catastrophic disruption
130
on the impact parameter b = sin , where  is the angle between the projectile's
path and the target's center at impact (see Section 4.3.2). Previous studies (Canup,
2008; Leinhardt & Stewart, 2012) have shown that the increase in the threshold
for catastrophic disruption for oblique impacts is due to a reduction in interacting
projectile mass. These authors provide a formulation parameterized by the frac-
tion of interacting mass, . We show that this does not account adequately for the
increase in the catastrophic dispersal threshold for impacts with b 6= 0 but   1
(impacts where most of the projectile interacts with the target). We discuss a possi-
ble revision to the formulation of the catastrophic disruption variables that includes
the eective interacting target material. By only taking into account the mass of
material that interacts in the collision, oblique impacts are rescaled into equivalent
head-on collisions such that they are well described by the so-called \universal" law
for catastrophic disruption (Leinhardt & Stewart, 2012).
Our results have important implications for the eciency of planet formation
via planetesimal growth, and for the determination of the impact energy thresh-
old for catastrophic disruption, as this has generally only been evaluated for non-
spinning bodies. In Section 4.3 we explain the computational methods and outline
the parameter space that we explore. In Section 4.4 we provide our results. In
Section 4.5 we discuss these results in the context of the \universal" law for catas-
trophic disruption and formulate a semi-analytic description of the dependence of
catastrophic disruption on pre-impact rotation. We summarize and oer perspec-
tives in Section 4.6.
131
4.3 Methodology
4.3.1 Rubble-Pile Model
Our simulations consist of two bodies with a mass ratio of  1 : 10: a station-
ary target with mass Mtarg and a projectile with mass Mproj (Mproj = 0.1Mtarg for
this work) which impacts the target at a speed of vimp. Both the target and pro-
jectile are gravitational aggregates of many particles bound together by self-gravity.
The particles themselves are indestructible and have a xed mass and radius. In
most of the simulations reported here, the only friction that is modeled is static
friction, for which we assume s = 0:5, corresponding to an internal angle of friction
of tan1(s)  27. However, we do discuss the dependence on material properties
and SSDEM parameters in Section 4.5.3.
The rubble piles are created by placing equal-sized particles randomly in a
spherical cloud and allowing the cloud to collapse under its own gravity with highly
inelastic particle collisions. Randomizing the internal structure of the rubble piles
reduces articial outcomes due to the crystalline structure of hexagonal close pack-
ing (Leinhardt et al., 2000; Leinhardt & Richardson, 2002). Due to symmetry lines
and planes in crystalline packing, there is a dependency of the collision outcome
on the initial orientation of the target's principal axes. To test the dependence of
the collision outcome on initial orientation for a spherically collapsed rubble pile,
a series of simulations was performed where the simulation parameters were kept
constant except for the initial orientation of the target's equatorial principal axes,
132
which were varied by increments of 45 about its polar axis. The results of these
simulations show that, for a spherically collapsed rubble pile, the dependence of
collision outcome on initial orientation is small (mass loss deviations of less than 1%
from the mean).
For the simulations presented here, the target had an average radius of Rtarg 
1:0 km and bulk density of targ  2 g cm3. The projectile had an average radius
of Rproj  0:5 km and bulk density of proj  2 g cm3. In order to determine
accurately the physical properties (size, shape, mass, angular momentum) of the
target after the collision, the rubble piles were constructed with a relatively high
number of particles (Ntarg = 10
4; Nproj = 10
3).
The collisional properties of the constituent particles are specied prior to each
simulation. These values were xed at n = 0:8 (mostly elastic collisions with some
dissipation) and t = 1:0 (no sliding friction). Furthermore, since SSDEM models
treat particle collisions as reactions of springs due to particle overlaps, the mag-
nitude of the normal and tangential restoring forces are determined by the spring
constants kn and kt =
2
7
kn, respectively. For rubble-pile collisions, the value of kn is
given by
kn  m

vmax
xmax
2
; (4.1)
where m corresponds to the typical mass of the most energetic particles, vmax is the
maximum expected speed in the simulation, and xmax is the maximum fractional
particle overlap, which we set to be  1%. Thus, for our rubble-pile collisions with
speeds  10 m s1, kn  41011 kg s2. The initial separation of the projectile and
133
Figure 4.1: Schematic of two collision scenarios. Panel (a) shows a target
impacted in its equatorial plane by a projectile moving left to right with
speed vimp. The impact angle  is the angle between the line connecting
the centers of the two bodies and the projectile's velocity vector, at the
time of contact. Panel (b) shows a projectile that impacts a target at an
angle  between the rotation axis (z) and the projectile's velocity vector.
target, d, for all cases was  4Rtarg, far enough apart that initial tidal eects were
negligible. In order for the post-collision system to reach a steady state, the total
run-time was set to  3 the dynamical time for the system, 1=pGtarg  2 hours.
Furthermore, a time-step 
t  3 ms was chosen on the basis of the time required
to sample particle overlaps adequately, for the choice of kn and xmax given above.
134
4.3.2 Simulation Parameters and Collision Geometries
In order to probe the eect of rotation on collision outcome, simulations with
the target rubble pile having an initial spin period Pspin of 3, 4.5, and 6 hours (values
well above the spin break-up limit for a rubble pile of bulk density  2 g cm3),
were compared against runs with the target having no initial spin. For every spin
period, simulations were done with a range of impact speeds such that there was
adequate coverage of the gravitational dispersal regime (collisions that result in a
system losing 0:1{0:9 times its total mass).
Furthermore, three dierent collision geometries were explored in this work,
each of which depended on two dierent collision parameters. The rst was the
impact parameter b = sin , where  is the angle between the line connecting the
centers of two bodies and the projectile's velocity vector (see Fig. 4.1a). The second
parameter was the angle , which is the angle between the target's rotation axis
and the projectile's velocity vector (see Fig. 4.1b ). In this study, the eect of each
parameter on the collision outcome was studied seperately and compared against
the standard case of a head-on collision. In a head-on collision, the projectile's ve-
locity vector is normal to the target's rotation axis and is directed towards its center
(b = 0 and  = 90).
For oblique impacts, the impact parameter b 6= 0. The impact parame-
ter has a signicant eect on the collision outcome because the total mass of
the projectile may not completely intersect the target when the impact is oblique
(Yanagisawa & Hasegawa, 2000; Canup, 2008; Leinhardt & Stewart, 2012). Thus,
135
when the projectile is large enough compared to the target, a portion of the pro-
jectile may shear o and only the kinetic energy of the interacting fraction of the
projectile will be involved in disrupting the target. For any given collision speed,
an oblique impact will erode less mass than a head-on collision. In this study, four
values of b were used,  0:5 and  0.7. For b > 0, the projectile impacts the target
on the hemisphere that rotates towards the projectile, which we dene as the ret-
rograde hemisphere (see Fig. 4.1a). For b < 0, the projectile impacts the target on
the hemisphere that is rotating away from the projectile, the prograde hemisphere.
Hence, if the collision outcome is sensitive to initial rotation, then it is expected
that the sign of b will also aect the amount of mass that is dispersed.
For non-equatorial impacts, the polar angle  < 90 (see Fig. 4.1b). If the
collision outcome is sensitive to the target's pre-impact rotation, then material from
the target's equator may preferentially be dispersed due to its lower specic binding
energy. However, it is uncertain whether the projectile more eciently transfers its
energy to the target's equator or to its poles upon impact. Hence, this study tests
the eect of three dierent polar impact angles:  = 90 (collisions directed at the
target's equator),  = 45, and  = 0 (collisions directed at the target's pole). In
reality, most collisions will have a combination of non-zero values for both b and .
4.4 Results
The collision of two rubble-pile objects typically results in either net accre-
136
tion, where the largest remnant has a net gain in mass compared to the mass of
the target, or net erosion, where the target has lost mass. Alternatively, a collision
could result in no appreciable net accretion or erosion (Leinhardt & Stewart, 2012).
These latter types of collisions, called hit-and-run events, typically occur for grazing
impacts that have an impact parameter, b, that is greater than a critical impact
parameter bcrit (Asphaug, 2010), where bcrit = Rtarg=(Rproj+Rtarg). In this chapter,
we focus on the dispersive regime, where impact velocities, vimp, are greater than
the escape speed from the surface of the target, vesc (assuming no rotation). The
impact speeds in our simulation range from 4{30 vesc, where vesc  1 m s1 is the
escape speed from a spherical object with mass Mtot = Mproj +Mtarg and density
1 = 1 g/cm
3. At impact speeds of 4{30 vesc, the mass of the largest remnant in
each simulation, MLR, ranges between 0:2{0:8 Mtot. The amount of mass loss at the
end of a simulation is found by measuring the nal mass of the largest remnant and
all material gravitationally bound to it (material with instantaneous orbital energy
 0). Furthermore, we analyze the mechanics behind rotation-dependent mass loss
by comparing the number of escaping particles that originate from dierent regions
of the target. The result of each simulation is summarized in Table 4.1.
4.4.1 Head-on Equatorial Collisions
For the nominal case of an equatorial-plane head-on collision, b = 0 and
 = 90. Fig. 4.2 shows that the amount of mass dispersal increases monotoni-
137
Figure 4.2: Mass loss for head-on equatorial impacts (b = 0,  = 90)
is sensitive to pre-impact rotation. The amount of mass that is gravita-
tionally dispersed is proportional to the impact speed. Green triangles,
blue stars, and purple squares represent impacts where the target has
a pre-impact spin period of 6 h, 4.5 h, and 3 h, respectively. Red-lled
circles represent impacts where the target has no pre-impact spin. For
head-on equatorial collisions, we derive the dependence of the reduced
mass catastrophic disruption threshold, Q?RD, on the target's pre-impac
rotation rate in Section 4.2.
138
cally with collision speed, and that, for head-on equatorial collisions, the amount
of mass dispersal is sensitive to initial rotation, as cases with shorter spin periods
systematically result in more mass loss. Furthermore, the collision outcomes for the
target with an initial period of 6 hours approach the outcomes where the rubble pile
initially has no spin. Since the spin limit for cohesionless rubble piles of this size and
density is  2:3 hours, these results are fairly representative of all possible head-on
collisions with initial spin (for s = 0:5). Through a simple linear regression of the
mass loss as a function of impact speed, we nd that the catastrophic disruption
threshold Q?D decreases by a range of  10{30% for the cases with pre-impact spin
studied here. Since the transition from merging to catastrophic disruption may oc-
cur over a dierence in energy of  30% (Leinhardt & Stewart, 2012), our results
show that pre-impact spin can play a crucial role in the formation of planetesimals
and protoplanets.
In order to obtain a better understanding of the underlying mechanics of
rotationally enhanced mass loss, we considered the geometrical eects associated
with a collision. By tracking the provenance of escaping particles originating from
the target, we studied the likelihood of a particle's escape as a function of its initial
longitudinal and latitudinal point of origin on the target. Fig. 4.3 shows a mass-
loss map for the case of a head-on collision with specic impact energy close to
catastrophic disruption (vimp = 9 m s
1). The collision creates an extended impact
region proportional to the projectile's size (Rproj  0:5 km), shown in Fig. 4.3 as the
black region in the middle of the map. Material within this region is retained by the
largest remnant as it is enclosed between the incoming projectile and the target's
139
Figure 4.3: The cylindrical equi-distant map projection of the areas of
mass loss from the head-on equatorial impact with vimp = 9 m s
1. For
an impactor of nite size, the regions near the impact point [coordi-
nates (0; 0)] do not experience mass loss as they are conned between
the antipodal material and the incoming projectile during the collision.
Rather, a ring of material surrounding the extended impact region is
ejected (here shown as a rectangular region due to the distortion of the
map projection).
140
Figure 4.4: The mass loss distribution with impact longitude of cases
with pre-impact spin normalized by the case without spin for the near-
catastrophic impact speed, vimp = 9 m s
1. Mspin is the mass in the
largest remnant originating from a given longitude range of the target.
The value of Mspin at each longitude bin is normalized by the case with
no pre-impact rotation, Mnospin. The sensitivity of mass loss to pre-
impact rotation exists due to an enhancement in the amount of escaping
material from the prograde hemisphere, which more than compensates
for a corresponding greater retention of material from the retrograde
hemisphere.
141
antipodal region (material 180 from the impact point). The escaping material orig-
inates from a nearly symmetrical ring about the impact region. Closer inspection
of this escaping material reveals that the enhancement in mass loss is due to a pref-
erential escape of material from prograde (negative) longitudes (Fig. 4.4 ). For this
analysis, particles in the initial rubble-pile target were binned into 30 longitudinal
spherical wedges. Since the target is nearly spherically symmetric, we consider lon-
gitudinal bins of equal size. Fig. 4.5 shows the number and sign convention that is
used. The 0 longitude point is dened as the meridian of the target aligned with
the impactor's velocity vector at the beginning of the simulation. Fig. 4.4 shows
the mass distribution of the largest post-impact remnant as a function of the par-
ticle origin. When a target has pre-impact spin, material located in the retrograde
longitudes (positive values) is preferentially retained by the largest remnant, and
material located in the prograde longitudes preferentially escapes. This is due to
retrograde material having angular velocities that are anti-aligned with the impact
velocity vector, and prograde material having angular velocities that are aligned.
However, there is a slight imbalance between these two outcomes, which favors a
net increase in escaping mass. This can be seen in the increase in the amplitude
of the troughs of Fig. 4.4 compared to the peaks. This phenomenon is observed in
all cases of pre-impact spin. It is this process that eventually determines the total
net enhancement in mass loss as a function of increasing pre-impact spin seen in
Fig. 4.2.
However, head-on collisions are not the most likely collision geometry; rather,
a collsion with b = 0:7 ( = 45) is the most common on average (Love & Ahrens,
142
Figure 4.5: In order to determine the longitudinal origin of escaping par-
ticles, the target is divided into equal-size bins of longitudes. Following
the right-hand rule, negative values of the longitude correspond to neg-
ative values along the y-axis when the spin-angular momentum vector is
aligned with the positive z-axis.
1996). Furthermore, the mechanics of mass-loss enhancement due to pre-impact
spin appears to be linked directly with the manner in which loading on the target
occurs. Therefore, we extend our analysis to study the dependence of mass loss with
the combined eects of non-zero b and  by studying each in isolation.
4.4.2 Oblique Equatorial Collisions
Previous studies have shown that the catastrophic disruption criterion is highly
sensitive to the impact parameter (Leinhardt et al., 2000; Leinhardt & Stewart,
2012). For oblique impacts, the energy of the projectile may not completely in-
tersect the target. For certain values of b, a segment of the projectile may be able to
shear o, and, consequently, this material does not interact with the target, eec-
143
Figure 4.6: Mass loss distribution for oblique equatorial impacts for
vimp = 15 m s
1. Panel (a) shows the distribution for impacts onto the
prograde hemisphere (negative b), where the solid lines are for collisions
with b = 0:5, and dotted lines are for b = 0:7. Panel (b) shows
the distribution for impacts onto the retrograde hemisphere (positive b),
where the solid lines are for collisions with b = +0:5, and dotted lines
are for b = +0:7. See text for discussion.
144
tively lowering the specic impact energy. Hence, a greater impact speed is required
to reach catastrophic disruption compared to a head-on collision. Previous studies
used a simple geometric model to determine the fraction of interacting projectile
mass (Canup, 2008; Leinhardt & Stewart, 2012). The revised mass is then used in
scaling the catastrophic disruption criteria (Leinhardt & Stewart, 2012).
We nd that mass dispersal is only sensitive to pre-impact rotation when b < 0,
for the range in b-values considered. For b > 0, the projectile impacts the target
on its retrograde hemisphere; hence, the angular velocities of rotating particles are
either anti-aligned or perpendicular to the impact velocity vector. Fig. 4.6 shows
the mass distribution as a function of longitudinal origin for vimp = 15 m s
1. The
solid and dotted curves are for b  0:5 and b  0:7, respectively. Fig. 4.6a shows
that, for collisions onto the prograde hemisphere, pre-impact rotation systematically
increases mass loss. This is similar to the results found for head-on collisions in Sec-
tion 4.4.1. Conversely, Fig. 4.6b shows that collisions onto the positive retrograde
hemisphere result in a rough balance between mass retention and mass-loss enhance-
ment when pre-impact spin is introduced (i.e., the integrals under the curves are
approximately zero). Hence, the net eect is that, for prograde impacts, pre-impact
spin does not enhance mass loss very much. Furthermore, for the oblique impacts
considered here, the rubble pile is ecient at dissipating the impact energy such
that gravitational dispersal is localized to the hemisphere of the target that was
impacted, rather than being a global eect.
145
4.4.3 Head-on Collisions with  < 90
Since the eective acceleration of a particle on a rotating body is a function
of its colatitude, we would expect that the amount of mass loss a rubble pile experi-
ences is a strong function of the latitude of impact. Fig. 4.7 shows the mass dispersal
as a function of impact speed for two dierent values of the polar impact angle, .
Contrasting with the head-on equatorial case, for low collision speeds, pre-impact
rotation does not systematically increase the amount of mass loss. However, for
near- and super-catastrophic speeds, the amount of mass loss is greatly enhanced.
For such energetic collisions, a target with a 3 h pre-impact spin period experiences
 25% (for  = 45) and  50% ( = 0) more mass loss than a non-spinning target.
On average, a pole-on impact onto a spinning target requires  15% less specic
impact energy to reach catastrophic disruption compared to an equatorial head-on
impact ( = 90, b = 0).
The enhancement in mass loss for near-polar impacts is due to the escape of
equatorial surface particles. As discussed in Section 4.4.1 and shown in Fig. 4.3, the
material that escapes originates from a ring about the extended impact region. For
the case of pole-on and near-pole-on impacts, this dispersal region extends to the
equator, where the eective pre-impact acceleration of particles is greatest. Unlike
equatorial impacts ( = 90), material at the equator is not conned by impact-
ing projectile material. Instead, for polar impacts, polar material is trapped by
the merging projectile material, and the collisional wave disperses material outside
the immediate polar region. For highly energetic collisions, this mechanical wave
146
extends to the equator, where particles are unhindered and more readily escape.
Hence, for the cases of near- and super-catastrophic collisions, the vertical transfer
of impact energy (from the pole to the equator) leads to the ejection of low-latitude
particles. Experiencing a lower eective gravitational potential at higher rotation
rates, these particles escape more easily. This leads to the enhancements in mass
loss seen in Fig. 4.7.
4.5 Discussion
4.5.1 The \Universal" Law for Catastrophic Disruption
In order to account for the dependence of mass ratio on catastrophic disrup-
tion criteria, Leinhardt & Stewart (2009) introduced new variables into their for-
mulation for predicting collision outcomes: the reduced mass   MprojMtarg=Mtot,
the reduced-mass specic impact energy QR  0:5v2imp=Mtot, and the correspond-
ing reduced-mass catastrophic dispersal limit Q?RD. Through this new formulation,
Leinhardt & Stewart (2009) showed that the outcome of any head-on collision, re-
gardless of projectile-to-target-mass ratio, can be described by a single equation that
they call the \universal" law:
MLR=Mtot = 0:5(QR=Q0RD) + 0:5; (4.2)
where the prime (0) notation in Q0?RD was introduced by Leinhardt & Stewart (2012)
to denote a collision that could have a non-zero impact parameter. Leinhardt & Stewart
147
Figure 4.7: Mass of the largest remnant as a function of impact speed
for cases with  < 90. Signicant enhancements in mass loss only
occur for near-catastrophic and super-catastrophic collisions. See text
for discussion.
148
Figure 4.8: Head-on collision outcomes are described well by the \univer-
sal" law for catastrophic disruption, Eq. 4.2 (dotted-line). Most head-on
collisions (red circles) show < 1% deviations from the universal law.
Non-equatorial impacts ( < 90) onto rotating targets, show deviations
of up to  10% due to non-linear eects present when collisions are near-
and super-catastrophic (see text for discussion). Oblique impacts (lled
triangles) systematically deviate from the universal law.
(2009) veried that, for head-on impacts (b = 0), Eq. 4.2 agrees well with results
from both laboratory experiments and numerical simulations of binary collisions
with a range of mass ratios and material properties. Leinhardt & Stewart (2012)
found that their numerical simulations showed deviations inMLR=Mtot of  10% for
near-normal impacts (b = 0 and b = 0:35), and larger and more varied deviations
for more oblique impacts. In order to compare our results (summarized in
Table 4.1) to those of Leinhardt & Stewart (2012), the outcomes of each group of
collisions (MLR as a function of QR) were t with a linear function to determine
149
empirically the value of Q0?RD (collision groups are uniquely identied by single val-
ues of b, , and Pspin.) For most collision scenarios the results are well described
by the linear t. In these cases, deviations from the model were less than 1%. For
non-equatorial collisions ( < 90) onto a rotating target, we showed in Section
4:4:3 the existence of two dierent mass loss outcome regimes. For near- and super-
catastrophic collisions, enhanced mass dispersal from the equator causes a non-linear
increase in mass loss. Hence, for these cases, a single linear t to determine Q0?RD
leads to deviations of  10%.
The results of our simulations are presented in Fig. 4.8, superimposed on the
universal law for catastrophic disruption. For collisions with b = 0, we observe that
the law predicts accurately the mass of the largest remnant (deviations < 10%).
Furthermore, varying the polar-impact angle, , or the spin period, Pspin, does not
aect the mass of the largest remnant if the specic impact energy is normalized by
Q0RD.
Cases of b = 0:5;0:7 seem to be better t by a shallower slope. For oblique
impacts, Leinhardt & Stewart (2012) introduced a parameter   mint;proj=Mproj,
which accounts for the mass of the projectile that interacts with the target, such
that the appropriate reduced mass is
  MprojMtarg
Mproj +Mtarg
: (4.3)
For the case of a 1:10 mass ratio and b = 0:5, we nd   1:0. Yet, we observe in
our experiments that the specic impact energy required for catastrophic disruption
is much greater than that for a head-on collision. Therefore, there must be some
150
other mechanism that accounts for this discrepancy in required specic impact en-
ergy.
During a collision, a compressive wave travels through the projectile, and upon
encountering the projectile's edge, is re
ected as a tensile wave that travels through
the contact points and disrupts and disperses target (as well as projectile) material
(Ryan, 2000). The width of this wave depends on a number of physical parameters
such as the strain loading rate, the type of material, and the size ratio of the projec-
tile and target. We hypothesize that the size of this wave determines what fraction
of the target interacts in the collision. For cases with non-zero b, we observe that
the particles that are able to escape are mostly localized to the same hemisphere of
the collision; hence, some of the material in the opposite hemisphere of the target
shears o and is not involved in the collision. For head-on collisions, the dispersive
wave would originate near the center and propagate symmetrically through the tar-
get, maximizing the amount of material aected by the collision. The exact wave
mechanics involved in computing the fraction of the target that does interact during
the collision is dicult to determine accurately analytically. For now, we attempt
to nd an empirical determination of this interacting target mass by adjusting the
disruption criteria variables to include this eect. Hence, we introduce a further
revised reduced mass:
A  MprojAMtarg
Mproj + AMtarg
; (4.4)
where A  mint;targ=Mtarg, the fraction of the mass of the target that interacts
in the collision. Thus, the equivalent head-on specic impact energy is QAR =
151
0:5Av
2
imp=Mtot. Furthermore, the mass of the largest remnant is now normalized by
a total interacting mass,MAtot  MprojAMtarg. Using these new variables, we rescale
the oblique impacts into equivalent head-on impacts. By using a 2-minimization
routine, we determined the best-t values for  and A (Fig. 4.9) that rescale the
oblique impact results such that they t the universal law for catastrophic disrup-
tion. We consider two cases, b = 0:5 and b = 0:7. Since the eect is considered
to be purely geometric, it is independent of the sign of b. For both cases of b, we
nd single 1=2 peaked regions in the parameter space. Possible values of A are well
constrained between 0:8 and 0:9 for both cases of b. The best-t values of  diers
for the two jbj cases, with b = 0:5;0:7 peaked at   0:7; 0:75, respectively.
In order to constrain the possible values of  and A, we also determined the
interacting mass fraction of the projectile, geom, through a simple geometric model.
Following Leinhardt & Stewart (2012), we determine geom by considering the vol-
ume of the projectile whose cross section intersects with the target. The volume of
a spherical cap can be expressed as:
Vint;proj =
l2
3
(3r  l); (4.5)
where l is the height of the spherical cap, and r is the radius of the sphere. For
an oblique impact, l is the projected length of the projectile overlapping the target,
and can be expressed as
l = (1 b)(Rtarg +Rproj): (4.6)
Therefore,
geom =
mint;proj
Mproj
=
Vproj
Vint;proj
=
3rl2  l3
4r3
: (4.7)
152
Figure 4.9: Using a 2 analysis, we determined the best values of  and
A that adjust the catastrophic disruption variables such that oblique im-
pacts are well modeled by the universal law for catastrophic disruption.
The white dashed line is the geometric constraint placed by geom. For b
=  0:7, geom = 0:44. The 2 analysis shows that the geometric model
may be underestimating the fraction of interacting projectile material
(see text for discussion).
153
The value of geom gives a lower boundary to the possible value of , since it is
expected that the minimum amount of interacting material would be the mass that
overlaps the target geometrically. We determined the best values of  and A as
constrained by the geometric model, and nd that as jbj increases,  and A decrease
(Table 4.2). This trend is expected, since more projectile and target material can
shear o when the impact is close to grazing. For the cases with b = 0:7, we nd
that the best-t value of  (0.745) is much greater than geom (0.44). This suggests
that the simple geometric model underestimates the value of . A fraction of the
projectile whose cross section does not overlap with the target does not completely
shear o; rather, it is involved in the collision process, contributing to the impact
energy delivered to the target.
For the cases with b = 0:5, we nd that the best-t values for  and A are
excluded when  is imposed to be greater than or equal to geom. If the geometric
constraint is not considered, then a combinaton of  = 0:685 and A = 0:845 gives
the minimun 2. However, this would imply that the cases of b = 0:5 have a lower
value for  than the cases of b = 0:7 ( = 0:745). This is unlikely, as a larger
fraction of projectile material is expected to interact when a collision is closer to
head-on.
Hence, while our 2 analysis does a good job of tting the data to the uni-
versal law, the results are unphysical unless they are constrained by a geometric
model. Therefore, more work must be done in order to verify whether the energet-
ics of oblique impacts are aected by the interacting fractional mass of both the
target and the projectile as we suggest here. In Fig. 4.10, we show that oblique
154
impacts follow the linear universal law if the axes are changed to our rescaled catas-
trophic disruption variables. However, we nd that for super-catastrophic collisions,
the data points seem to tail-o rather than follow a linear relationship. Labora-
tory experiments (Matsui et al., 1982; Kato et al., 1995) and disruption simulations
(Korycansky & Asphaug, 2009; Leinhardt & Stewart, 2012), have shown that the
mass of the largest remnant follows a power law with QR for super-catastrophic
collisions. At these high energies, mass dispersal results in the formation of a large
number of fragments of roughly equal size, rather than a single large remnant, such
that for incrementally higher impact energies, the largest remnant remains constant.
This may explain the discrepancy that we see between the high-energy collision out-
comes and the universal law.
4.5.2 Rotation Dependence of Catastrophic Disruption
Since rotation decreases the eective gravitational binding energy of a body,
we rst describe the size-dependence of catastrophic disruption so that we may be
able to formulate an analytic description of the dependence on pre-impact rotation.
The catastrophic disruption criterion is a function of radius, with two regimes: a
strength-dominated regime and a gravity-dominated regime (Housen & Holsapple,
1990). For rocky bodies, the transition from strength to gravity occurs at a radius of
 100 m (Leinhardt et al., 2008). In the strength regime, the catastrophic disrup-
tion criterion decreases with increasing radius; this is due to multiple factors, such
as the increase in the size of the largest internal crack and the total number of 
aws
155
Figure 4.10: Adjusting our results for the correct interacting projectile
and target produces a better t for oblique impacts, for low QAR=Q
?
RD
(cf. Fig. 4.8). The corrected interacting mass fractions are summarized
in Table 4.2.
156
with target size. In the gravity regime, disruption increases as the radius increases
since disruption requires shattering and gravitational dispersal, and the gravitational
binding energy of a body, U , is proportional to the square of the body's radius. For
a binary collision, the gravitational binding energy can be approximated as,
U =
3GMtot
5RC1
=
4
5
1GR
2
C1; (4.8)
where G is the gravitational constant and RC1 is the spherical radius of the combined
projectile and target masses at a density of 1 = 1 g cm
3. Leinhardt & Stewart
(2009) introduced RC1 in order to compare collisions of dierent projectile-to-target-
mass ratios.
By determining the dependence of mass ratio on the catastrophic disruption
criterion, Leinhardt & Stewart (2012) found that, in the gravity regime, the disrup-
tion criterion for equivalent equal-mass impacts, Q?RD, of dierent materials all fall
along a single curve that scales as the radius squared. Since catastrophic disruption
and the gravitational binding energy scale similarly with radius, the authors dene
a principal disruption curve, where Q?RD is a scalar multiple of U . They dened a
dimensionless material parameter, c?, that represents this oset, such that
Q?RD = c
?U: (4.9)
For bodies with a diverse range of material properties (strengthless hydrodynamic
targets, rubble piles, ice, and strong rock targets) and with 0:5 km < RC1 < 1000
km, Leinhardt & Stewart (2012) found that the threshold for catastrophic disruption
is dened by a single principal disruption curve with c? = 5  2. They nd that
157
a similar curve can describe the disruption of planet-size bodies (RC1 > 3000 km)
with c? = 1:9 0:3.
In order to formulate a description of the dependence of catastrophic disruption
on rotation, we rst consider a rotating 
uid body that has an eective specic
gravitational binding energy, Ue , given by
Ue = U  j!rj
2
2
; (4.10)
where ! is the constant angular velocity of the body, and r is the position vector of
a particle relative to the center of the target. Similarily, we propose that in a binary
head-on equatorial collision with pre-impact rotation, the catastrophic disruption
criterion, Q?RD;rot, is a function of the angular speed of the target, !targ. For a
head-on equatorial collision with pre-impact rotation (collision geometry with no
dependence on b and ), this is represented by a substraction of a latitude-averaged
centrifugal term as in Eq. 4.10 such that
Q?RD;rot = Q
?
RD;no-rot  !2targR2targ; (4.11)
whereQ?RD;no-rot is the catastrophic disruption criterion without pre-impact spin. For
non-equatorial or oblique impacts, we expect that the change in the catastrophic
disruption criterion is a more complicated function of collisional angular momentum
and the eciency of its transfer. For the purposes of this chapter, we restrict our
analysis to the simpler equatorial head-on collisions. Dividing Eq. 4.11 by Q?RD;no-rot,
and substituting from Eq. (9), we nd
158
Q?RD;rot
Q?RD;no-rot
= 1K

!targ
!crit
2
;
K  5
3
RC1
Rtarg
Mtarg
Mtot
1
c?
;
(4.12)
where !crit  (GMtarg=R3targ)1=2 is the spin break-up limit of the target. We t
Eq. 4.12 with the empirically derived catastrophic disruption values (described in
Section 4.5.1) for the three dierent pre-impact spin cases normalized by the spinless
case (Fig. 4.11). We nd a best-t value of K = 0:3814 (maximum deviations were
 3%), corresponding to a value of c? = 4:83. Our value of c? falls within the range
for small bodies (R < 1000 km in size) that Leinhardt & Stewart (2012) nd.
Eq. 4.12 predicts that for pre-impact rotations close to break-up, the catas-
trophic disruption criterion can decrease by a factor of  40%. In reality, a spherical
target spinning close to break-up would reach a new 
uid-equilibrium shape, and the
nature of its ellipsoidal shape will likely aect the catastrophic disruption criterion.
In our simulations, the impacts occur quickly enough that the target does not reach
a 
uid equilibrium before disrupting. Future work could study the dependence of
catastrophic disruption on the pre-impact shape of the target body.
4.5.3 Dependence on Material Properties
Finally, we measure the dependence of Q?RD on material properties. Our pre-
vious simulations kept the frictional properties of the particles constant. Here, we
wish to examine the range of possible Q?RD that may be obtained by varying the
friction coecient values. We model three dierent types of materials: smooth par-
159
Figure 4.11: Head-on equatorial impacts have catastrophic disruption
thresholds that are sensitive to pre-impact spin. We develop a semi-
analytic description of Q?RD as a function of the pre-impact spin-rate
of the target, !targ (see text), and nd that our data imply a value
of c? = 4:83 for the dimensionless material parameter of the principal
disruption curve. The dashed curve is the best-t function (of the form
described by Eq. 4.12) to the catastrophic disruption thresholds with
pre-impact rotation normalized by the case with no rotation (circles).
160
ticles, glass beads, and gravel. The properties of glass beads were determined by
Schwartz et al. (2013) in impact experiments. The properties of gravel were deter-
mined by Yu et al. (2014) in landslide experiments. The dissipation and friction
parameters that correspond to each material were previously discussed in Chapter
2 and are shown in Table. 2.1. We performed simulations of binary collisions where
we varied the material property of the rubble pile and the pre-impact spin rate of
the target. In Fig. 4.12, we show that a rubble pile made up of rougher and less elas-
tic material has a higher catastrophic disruption threshold. Furthermore, we also
see that regardless of the material property type, the decrease in Q?RD due to some
non-zero pre-impact spin is constant. Therefore, we expect that the semi-analytic
prescription developed in the previous section to be valid for a rubble pile made up
of particles of any material property.
4.6 Conclusions & Future Work
We have studied the eect of initial rotation on the outcome of rubble-pile col-
lisions by analyzing the properties of the largest remnant and material that is grav-
itationally bound to it. By simulating dierent collision geometries and speeds, we
have begun to explore a parameter space that is wide enough that we can formulate
a phenomenological description of collision outcomes. Our main conclusion is that
mass dispersal is a function of initial rotation period, with faster-rotating rubble-pile
targets dispersing more mass. By analyzing the initial spatial distribution of the
gravitationally bound masses, we have shown that there is an enhancement of mass
161
Figure 4.12: The reduced-mass catastrophic disruption threshold Q?RD,
normalized by the gravitational binding energy, U , for the three material
properties studied here. Red circles represent collisions with no pre-
impact spin. Blue triangles represent collisions where the target has a
pre-impact spin of 6 hours. Spin systematically decreases Q?RD=U by 
6%.
162
loss when the impact energy is eciently transferred to the prograde hemisphere of
a rotating rubble pile. For head-on impacts onto regions near the pole, the collision
eciently disperses equatorial material for near- and super-catastrophic collisions,
and fast rotation increases mass loss by factors of up to 50%. The mass of the largest
remnant of head-on impacts is well described by the \universal law" for catastrophic
disruption rst put forward by Leinhardt & Stewart (2009), independent of initial
pre-impact rotation. Hence, for a given impact speed, pre-impact rotation decreases
MLR=Mtot, and the corresponding decrease in Q
0?
RD can be described by the uni-
versal law. However, oblique impacts follow a linear relationship with a shallower
slope than the universal law. When the interacting mass fraction of the projectile
and target are factored into the catastrophic disruption variables, the outcomes for
oblique impacts can be rescaled for a better match to the universal law.
By subtracing a centrifugal term from the catastrophic disruption criterion of
the case with no pre-impact rotation, we developed a prescription that describes
the change in the catastrophic disruption criteria of head-on equatorial impacts
onto a rotating target. We independently nd a dimensionless material parameter,
c? = 4:83, that agrees with the principal disruption curves of Leinhardt & Stewart
(2012) for small bodies. Our simplied description does not take into account the
eects of angular momentum transfer in oblique impacts or the mechanism for en-
hanced mass loss from near- and super-catastrophic polar impacts described in Sec-
tion 4.4.3. These eects will have to be further studied so that the change in the
catastrophic disruption criterion can be determined for any impact trajectory.
Lastly, we calculated the catastrophic disruption threshold for three dierent
163
material types. We found that higher friction and dissipation values increase the
catastrophic disruption threshold by about half a magnitude. Furthermore, we nd
that pre-impact rotation systematically increases mass loss on average, regardless
of the target's internal conguration.
In the future, a wider parameter space should be explored in order to strengthen
the conclusions drawn here. In particular, the eectiveness of rotation in enhancing
mass loss must be studied for dierent projectile-to-target-mass ratios. Since mass
loss is sensitive to rotation, the spin-up or spin-down of the post-impact largest
remnant plays an important role in the size evolution of a population of km-size
bodies. Furthermore, the change in spin likely has an eect on the reaccumulation
process, changing the nal mass. In some cases, spin-up may lead to a remnant
crossing the rotational disruption threshold (Pspin  2:3 h, for   2 g cm3). In
order to understand how rotation aects the long-term size, shape, and spin evo-
lution of a population of SSSBs, a semi-analytic description of the dependence of
Q0?RD on the parameters explored here must be formulated, such that MLR can be
determined for any given collision. Therefore, future work will also need to study
the sensitivity of mass loss on rotation for larger (> 1 km) bodies. It is unclear
whether the collisional dynamics explored here scales to larger bodies. This will
help inform future planet-formation studies by giving a more accurate prescription
for collision outcomes, a necessary component for models that study the collisional
growth of planetesimals.
164
Table 4.1: Summary of Collision Geometries and Mass Loss Outcomes. b is the
impact parameter,  is the polar impact angle, Pspin is the spin period of the target,
vimp is the impact speed of the projectile, and MLR=Mtot is the mass of the largest
remnant normalized by the total mass of the system.
b  () Pspin (h) vimp (m s1) MLR=Mtot
0 90 1 6.0 0.831089
0 90 1 7.0 0.764666
0 90 1 8.0 0.683734
0 90 1 9.0 0.607632
0 90 1 10.0 0.517345
0 90 1 11.0 0.428664
0 90 6 6.0 0.811708
0 90 6 7.0 0.742641
0 90 6 8.0 0.661843
0 90 6 9.0 0.579304
0 90 6 10.0 0.494084
0 90 6 11.0 0.389184
0 90 4.5 6.0 0.800477
0 90 4.5 7.0 0.728418
0 90 4.5 8.0 0.649403
0 90 4.5 9.0 0.567309
0 90 4.5 10.0 0.464814
0 90 4.5 11.0 0.351837
0 90 3 6.0 0.760627
0 90 3 7.0 0.696922
0 90 3 8.0 0.599811
0 90 3 9.0 0.514277
0 90 3 10.0 0.421599
0 90 3 11.0 0.305056
165
Table 4.1: Summary of Collision Geometries and Mass Loss Outcome, continued.
b  () Pspin (h) vimp (m s1) MLR=Mtot
0 45 1 5.0 0.876133
0 45 1 6.0 0.807526
0 45 1 7.0 0.739923
0 45 1 8.0 0.652526
0 45 1 9.0 0.575735
0 45 1 10.0 0.496087
0 45 6 5.0 0.877857
0 45 6 6.0 0.805796
0 45 6 7.0 0.725621
0 45 6 8.0 0.634082
0 45 6 9.0 0.554188
0 45 6 10.0 0.485355
0 45 3 5.0 0.846789
0 45 3 6.0 0.767449
0 45 3 7.0 0.684930
0 45 3 8.0 0.565474
0 45 3 9.0 0.441618
0 45 3 10.0 0.378795
0 0 1 5.0 0.881670
0 0 1 6.0 0.817889
0 0 1 7.0 0.734517
0 0 1 8.0 0.649638
0 0 1 9.0 0.562160
0 0 1 10.0 0.472488
166
Table 4.1: Summary of Collision Geometries and Mass Loss Outcome, continued.
b  () Pspin (h) vimp (m s1) MLR=Mtot
0 0 6 5.0 0.874406
0 0 6 6.0 0.812157
0 0 6 7.0 0.724058
0 0 6 8.0 0.627248
0 0 6 9.0 0.531327
0 0 6 10.0 0.433020
0 0 3 5.0 0.865095
0 0 3 6.0 0.777822
0 0 3 7.0 0.684112
0 0 3 8.0 0.481345
0 0 3 9.0 0.347406
0 0 3 10.0 0.206721
-0.5 90 1 5.0 0.855201
-0.5 90 1 10.0 0.694707
-0.5 90 1 15.0 0.513617
-0.5 90 1 20.0 0.299089
-0.5 90 6 5.0 0.834180
-0.5 90 6 10.0 0.671891
-0.5 90 6 15.0 0.490690
-0.5 90 6 20.0 0.312836
-0.5 90 3 5.0 0.797677
-0.5 90 3 10.0 0.638241
-0.5 90 3 15.0 0.445009
-0.5 90 3 20.0 0.289541
167
Table 4.1: Summary of Collision Geometries and Mass Loss Outcome, continued.
b  () Pspin (h) vimp (m s1) MLR=Mtot
+0.5 90 1 5.0 0.852571
+0.5 90 1 10.0 0.699535
+0.5 90 1 15.0 0.540372
+0.5 90 1 20.0 0.327492
+0.5 90 6 5.0 0.865060
+0.5 90 6 10.0 0.702961
+0.5 90 6 15.0 0.511956
+0.5 90 6 20.0 0.326386
+0.5 90 3 5.0 0.860442
+0.5 90 3 10.0 0.690634
+0.5 90 3 15.0 0.481739
+0.5 90 3 20.0 0.315096
168
Table 4.1: Summary of Collision Geometries and Mass Loss Outcome, continued.
b  () Pspin (h) vimp (m s1) MLR=Mtot
-0.7 90 1 5.0 0.876840
-0.7 90 1 10.0 0.796897
-0.7 90 1 15.0 0.715336
-0.7 90 1 20.0 0.619042
-0.7 90 1 25.0 0.532207
-0.7 90 1 30.0 0.424618
-0.7 90 6 5.0 0.863776
-0.7 90 6 10.0 0.786793
-0.7 90 6 15.0 0.695417
-0.7 90 6 20.0 0.601023
-0.7 90 6 25.0 0.504266
-0.7 90 6 30.0 0.414065
-0.7 90 3 5.0 0.840861
-0.7 90 3 10.0 0.750761
-0.7 90 3 15.0 0.656733
-0.7 90 3 20.0 0.551418
-0.7 90 3 25.0 0.462672
-0.7 90 3 30.0 0.386483
+0.7 90 1 5.0 0.874388
+0.7 90 1 10.0 0.795250
+0.7 90 1 15.0 0.712159
+0.7 90 1 20.0 0.628877
+0.7 90 1 25.0 0.534306
+0.7 90 1 30.0 0.432360
+0.7 90 6 5.0 0.880263
+0.7 90 6 10.0 0.804069
+0.7 90 6 15.0 0.720063
+0.7 90 6 20.0 0.628323
+0.7 90 6 25.0 0.525654
+0.7 90 6 30.0 0.390580
+0.7 90 3 5.0 0.884772
+0.7 90 3 10.0 0.787585
+0.7 90 3 15.0 0.702943
+0.7 90 3 20.0 0.599283
+0.7 90 3 25.0 0.477597
+0.7 90 3 30.0 0.359355
169
Table 4.2: Summary of 2 analysis: best-t interacting mass fractions. b is the
impact parameter,  and A are the interacting mass fraction of the projectile and
target, respectively, found the the chi2 minimization technique. 2 is a goodness-
of-t parameter, and geom is the interacting mass fraction of the projectile found
by comparing the cross-sections of the the projectile and the target.
b  A 2 geom
0:5 0.85 0.94 0.61 0.85
0:7 0.745 0.845 0.68 0.44
170
Chapter 5: Large-scale Structure in Saturn's B ring
5.1 Chapter Preface
The work presented in this chapter has been submitted for publication to the
Astrophysical Journal (ApJ) as Ballouz et al. (2017). In this chapter, we study
the B ring's complex optical depth structure. The source of this structure may be
the intricate dynamics of the Keplerian shear and the self-gravity of the ring par-
ticles. The outcome of these dynamic eects depends sensitively on the collisional
and physical properties of the particles. Two mechanisms can emerge that domi-
nate the macroscopic physical structure of the ring: self-gravity wakes and viscous
overstability. Here we study the interplay between these two mechanisms through
simulations of a local ring patch.
5.2 Background
Observations have shown that Saturn's rings exhibit rich and varied structure.
Images from the Cassini spacecraft show that this structure exists on many dierent
scales, from 100's of meters to 100's of kilometers (Porco et al., 2005; Colwell et al.,
2007). The B ring, in particular, exhibits a varied and complex structure as uncov-
171
ered by physical optical depth measurements (Colwell et al., 2009). This structure
likely comes about from the combined eects of ring particle collisions, ring self-
gravity, and tidal forces of Saturn. For the origin of structure of sub-km scales, two
mechanisms have been proposed.
The rst of these mechanisms is the formation of gravitational wakes in the
dense rings. Self-gravity causes the ring particles to clump, forming transient struc-
tures that are quickly disrupted by the background shear. These wakes have typical
wavelengths of a few 100's of meters (the wavelength is dened as the average radial
separation between two wakes) , and they are chie
y characterized by their pitch an-
gle of 20{30, the angle that the opaque clumps form with respect to the azimuthal
direction of motion. While they have never been imaged directly, their presence
have been inferred from Cassini Ultraviolet Imaging Spectrograph (UVIS), Visual
and Infrared Mapping Spectrometer (VIMS), and radio occultation optical depth
measurements of the A and B rings (Colwell et al., 2006; Hedman et al., 2007).
The self-gravitational interactions of the particles create a dense packing of
particles that are able to dissipate energy through inelastic collisions, leading to
an increase in the angular momentum transport in the ring. For a sucient num-
ber of dissipative collisions, the viscosity of the ring may evolve to the point that
a viscous overstability develops (Salo et al., 2001). This allows particles in over-
dense regions to 
ow towards under-dense regions. As particles 
ow radially away
to under-dense regions, the opacity of the ring patch is temporarily smoothed out
as particles don't preferentially occupy any region of space. However, this smooth-
ing process overshoots as particles moving radially inward towards Saturn from an
172
overdense structure collide with particles moving radially outwards from an adja-
cent structure. This allows new over-dense regions to develop. These pulsations in
ring density generate a wave pattern characterized by axisymmetric over-dense re-
gions (regions of high vertical coherence, see Hedman et al. 2014). Similar to gravity
wakes, these opaque clumps have sub-km wavelength, typically  100 m (Salo et al.,
2001; Daisaka et al., 2001; Rein & Latter, 2013). Recent Cassini infrared occulation
measurements of the B ring (Hedman et al., 2014) show that sub-km structure exists
and appear to be structurally stable and axisymmetric at scales of up to 3,000 km,
many orders of magnitude larger than the radial overstability wavelength, suggesting
that such overstable structures do indeed exist in the B ring.
The development of the viscous overstability depends sensitively on the colli-
sional dynamics of the ring particles. Initially, a viscous instability was hypothesized
to explain the structure in the B-Ring. In this scenario, particle collisions were as-
sumed to be highly elastic, which results in the dynamic shear viscosity of the ring
to decrease with density, and the particle 
ux will be directed towards regions of
high density. This results in an amplication of density 
uctuations. However,
laboratory experiments by Bridges et al. (1984) showed that collisions between fric-
tionless and spinless icy spheres at impact speeds typical of Saturn's dense rings
(a few mm/s) are highly inelastic. Subsequent laboratory and numerical studies
showed the rst hints that the viscous overstabiltiy could be possible in plane-
tary rings (Borderies et al., 1985; Salo et al., 2001). These laboratory experiments
(Bridges et al., 1984; Borderies et al., 1985) revealed that the restitution coecient
of the ice spheres was a function of the impact speed. This functional dependence
173
of the restitution coecient with the impact speed predicts an energy equilibrium
that allows for the viscous diusion of ring material.
The study of the self-gravity and viscous overstabiltiy mechanisms has been
advanced through numerical simulations using N -body integrators that simulate
the collisional and gravitational evolution of some 100,000 particles in a co-moving
patch (e.g., Salo 1992; Yasui et al. 2012; Rein & Latter 2013). These studies adapt
the functional dependency of the dissipative parameters on the impact speeds found
in Bridges et al. (1984) to model the collisions of hard-sphere particles. In the hard-
sphere method, particle collisions occur as pair-wise interactions and are resolved
instantaneously. These simulations have predicted that the viscous overstability
may only develop when the internal density of ring particles is low (200-300 kg/m3),
implying highly porous ring material or rubble-pile like structure. At larger parti-
cle densities, the viscosity of the ring is insucient to keep it stable against self-
gravitational perturbations, and gravity wakes develop. Studies have shown that
both viscous overstabiltiy and self-gravity wakes may co-exist for certain ring phys-
ical and collisional properties (Salo et al., 2001). However, the exact nature of the
relationship between these two mechanisms is not well known.
In this study, we perform direct numerical simulations of a patch of particles
in Saturn's dense B ring. Our N -body code allows us to study the collisional and
gravitational dynamics in high optical depth regions in Saturn's rings. The colli-
sional method we employ in this study, which allows multi-contact and frictional
forces (Schwartz et al., 2012), has never before been used in studies of Saturn's
rings. This state-of-the-art numerical modeling enables a better understanding of
174
the interplay between the two mechanisms that create the structure of Saturn's B
ring. The goal of this study is to generate a library of synthetic observations of the
optical depth structure mapped to particle physical and collisional properties. By
comparing these synthetic observations with actual data, we may be able to also
determine the physical properties of ring particles.
5.3 Methodology
We simulate the collisions of ring particles in Saturn's B ring, a distance of
100,000 km away from Saturn. The typical collision speeds of particles due to
Keplerian shear is expected to be of the order of a few times 
0rp, where 
0 is the
Keplerian orbital angular frequency for an object orbiting Saturn at 100,000 km
( 2  104 s1), and rp is the typical radius of the ring particles. Thus, for our
ring simulations, particles typically collide at relative speeds of  0:5  1 mm s1.
Therefore, in order to maintain overlaps that are  1%, we set kn  4:71 kg s2 and
the time step t  1:5 s.
5.3.1 Periodic Boundary Conditions & Patch Properties
We perform local simulations of ring particles by restricting the computational
volume to a small region (a \patch") of the entire ring. Particles are modeled in a
co-moving patch frame (Perrine & Richardson, 2012). The linearized equations of
motion, called the Hill's equations, are (Hill, 1878):
175
x = Fx + 3

2x+ 2
 _y;
y = Fy  2
 _x;
z = Fz  
2z;
where F is the acceleration due to particle self-gravity, x, y, and z are the
coordinates of the particle in the local coordinate system (whose origin is located at
the center of the patch), and the derivatives are with respect to time. We use peri-
odic boundaries at the radial and azimuthal edges of the patch. Periodic boundary
conditions (PBCs) allow us to simulate a patch of the disk with a computational
feasible number of particles. PBCs are imposed on a patch of particles by repli-
cating the patch in the radial (x) and azimuthal (y) directions. Each replicated
patch contains ghost particles that match the relative positions of the original par-
ticles. Figure 5.1 illustrates the PBC setup. Particles in the central patch feel the
gravitational eects of ghost particles, and can collide with ghost particles at the
boundaries. The azimuthal and radial extents of the patch are small compared to
the orbital distance to the planet, but large compared to the radial mean free path
of the particles inside the patch. Particles exiting one side of the azimuthal or radial
boundary reappear on the other side, preserving their general properties (mass and
spin). When a particle crosses the radial boundary, its azimuthal speed is modied
to account for the Keplerian shear. Therefore, depending on whether it crosses over
from a smaller radial position to a larger one or vice versa, its azimuthal speed will
be adjusted by 3
2

Lx, where Lx is the radial width of the patch.
176
!!
!
!
!
!
!
!
!
!
!
!
!
!
!
!
!
!
To Saturn 
x
y
Figure 5.1: A schematic diagram of a local sliding-patch model with
shearing periodic boundary conditions. The three white particles in the
center solid box are the simulated particles; grey particles are in repli-
cated patches (dashed boxes) that provide boundary conditions. The x
coordinate is the radial direction, with Saturn located in the negative x
direction; y is the azimuthal direction, and the entire patch orbits Saturn
in the positive y direction. z points out of the page, forming a right-
handed coordinate system. The simulation is carried out in the orbital
frame of the center of the patch, so particles to the left shear upward
(positive y direction); on the right, they shear downward (negative y
direction). The replicated patches similarly shear past the center box in
Keplerian fashion; each bold X marks the center of each patch, with the
bulk motion of each patch indicated by black arrows.
177
We simulate a patch of the B ring that has a dynamic optical depth of 1.4 and
a surface density of 840 kg m2. We choose this value of the optical depth as it is
typical of the B ring, and it corresponds to a number of particles that can be feasibly
simulated in a timely fashion using our code. The value of the surface density is
well within the range found through analysis of density waves (400{1400 kg m2,
Hedman & Nicholson 2016). In order to ensure that we allow sucient space for
the formation of large-scale structures, we simulate a patch that is at least 20 cr in
the radial direction and 4 cr in the azimuthal direction, where cr is the Toomre
wavelength, and is dened as:
cr =
42G


; (5.1)
where G is the gravitational constant and  is the surface density. Previous simu-
lations (Salo, 1992; Daisaka & Ida, 1999) have shown that the typical wavelengths
of wakes is close to cr.
For each simulation, we choose a single combination of a particle density and
particle radius that corresponds to a constant optical depth of 1.4 and a patch
surface density of 840 kg m2.This is done in order to keep cr constant across all
simulations, allowing us to draw more accurate comparisons of simulation outcomes
across dierent particle density cases. Therefore, a simulation can contain particles
with densities of 0.45, 0.60, 0.75, or 0.90 g/cm3 that have radii of 1.0, 0.75, 0.6, and
0.5 m, respectively. The number of particles in our simulations range from 120,000
to 500,000 particles (see Table 5.1)
178
5.3.2 Diagnostic Tools
We develop diagnostic tools to study the bulk behavior of the particles in the
patch and to compare our results to observations. Firstly, in order to measure the
dierent contributions to the angular momentum transport in the ring, we describe
how we calculate the viscosity of the patch. This is important for characterizing
the in
uence of the individual ring properties to the larger-scale dynamics, and it
allows us to directly compare our results to previous theoretical treatments of ring
dynamics. Secondly, we describe a Monte Carlo ray-tracing algorithm that measures
the physical optical depth of our simulations. This allows us to directly compare
the outcome of our simulations to Cassini results by producing synthetic stellar
occultation data.
5.3.2.1 Viscosity Calculation
In a planetary ring system, the viscosity of the disk describes how eciently
the system is able to transport angular momentum from interior parts to its exte-
rior. In its simplest sense, viscosity is a measure of a 
uid system's resistance to
stresses, either shearing or tensile. For a ring system such as Saturn's, there are
three components that contribute to the viscosity in a disk: translational viscosity,
collisional viscosity, and gravitational viscosity (Daisaka et al., 2001). The rst is
the translational component due to Keplerian shear, trans. This can be calculated
179
by measuring the radial velocity dispersion in the ring patch, and is dened as:
trans =
2
3
imi
imi _xi _yr;i; (5.2)
where imi is the sum of all particle masses, _xi is the x-component of the velocity,
and _yr;i = _yi + (3=2)
xi is the y-component of the velocity relative to the mean
shear speed at radial position xi.
The second component of the viscosity is due to the transfer of angular momen-
tum due to collisions. Previous authors (Wisdom & Tremaine, 1988; Daisaka et al.,
2001) have relied on event-driven hard-sphere collision codes. Thus, they calculate
the exchange in momentum in the azimuthal direction for each discrete collision.
The collisional component of the viscosity is given by:
coll =
2
3
MT
py(x>  x<); (5.3)
where M is the total mass in the patch, and the summation is over all pairwise
collisions that occur during the time interval T , which results in an exchange of
momentum py in the y (azimuthal) direction to the particle in the exterior orbit
with radial distance x> whose collisional partner has an interior orbit with radial
distance x<.
The last component of the viscosity is due to the self-gravitational interaction
of ring particles. While this is an important contribution to the overall viscosity of
a ring system, we do not measure this component in this study. Rather, the focus
of this study is on how collisional interactions of the disk are eected by dierent
physical properties of the particles. We refer the reader to Daisaka et al. (2001) for
a full discussion on gravitational viscosity.
180
Since our soft-sphere code treats collisions as nite duration with possible
multiple simultaneous contacts and frictional forces, we must revise this formulation.
Following the treatment in Wisdom & Tremaine (1988), we can reformulate the
exchange of momentum between two particles from the collisional force balance
calculated in our soft-sphere collision method, so the collisional viscosity becomes:
coll =
2
3
M
Fy>(x>  x<); (5.4)
where Fy> is the collisional force in the azimuthal direction on a particle due to
its colliding partner(s). Formulated in this way, we are able to take into account
multi-contact and frictional forces.
Another method that involves measuring the change in orbital elements due to
collisions has been devised to measure the total viscosity (sum of translational, colli-
sional, and gravitational viscosity) in the ring patch (Tanaka et al., 2003). Yasui et al.
(2012) recently used this method to calculate the viscosity of spinning self-gravitating
particles in planetary rings. They were also able to take the rst steps in measuring
the eect of friction on the viscosity by including the eects of a tangential resti-
tution coecient. Yasui et al. (2012) were able to demonstrate the validity of their
calculation by performing simple tests to measure the viscosity in a small patch and
compare their results to those found by Wisdom & Tremaine (1988).
In order to validate our re-formulation of the Wisdom & Tremaine (1988) pre-
scriptions for calculating the collisional viscosity, and to ensure that our code is
capable of replicating results found in recent numerical simulations of ring systems,
we performed a simple simulation test similar to that done by Yasui et al. (2012).
181
Using a square 125 m  125 m ring patch made up of 1 m particles, we ran ring
simulations with a speed-dependent normal coecient of restitution, n, based on
laboratory impact experiments (Bridges et al., 1984), given by:
n(v) = minf0:32(v=vc)0:234; 1g; (5.5)
where v is the relative speed between colliding particles at the moment of rst
overlap, and vc is a normalization equal to 1 cm s
1. We performed these simulations
for ring patches with optical depths ranging from 0.4 to 1.6, and for values of t of
0.5 and 1.0. These runs had the self-gravity between particles turned o since we are
only interested in the equilibrium viscosity achieved through ring particle collisions.
This also speeds up the calculation. The simulations were run until the particles
reached collisional equilibrium (a few 10's of orbits), determined when the patch
reaches a steady radial velocity dispersion. We nd excellent agreement between
the viscosities we measure and those found in Yasui et al. (2012) (see Figure 5.2).
5.3.2.2 Physical Optical Depth Calculation
The dynamical optical depth dyn of a ring is dened as
dyn =
Z
s2n(s)ds; (5.6)
where s is the particle radius and n(s) is the surface number density of ring particles
whose sizes are between s and s+ ds.
The physical normal optical depth phys is dened as (e.g., Colwell et al. 2007)
phys =  sin jBj log( I
I0
); (5.7)
182
Figure 5.2: We plot the total viscosity, , normalized by 
0r
2
p as a funci-
ton of the dynamical optical depth of the ring. Our total viscosity calcu-
lations (translational and collisional) match those found by Yasui et al.
(2012) (see their Fig. 1b). The simulations are for a velocity-dependent
normal restitution coecient with t = 0:5.
183
where B is the elevation angle of the observer, I is the measured intensity, and I0
is the unocculted intensity of the light source.
If particles are randomly distributed in the horizontal directions and the mu-
tual separation is large enough, phys becomes independent of observational geometry
and coincides with dyn (the limit of classic radiative transfer). If the mutual sep-
aration is comparable to the particle radius, phys becomes larger than dyn. This
eect is called the volume lling factor eect (e.g., Salo & Karjalainen 2003). If the
ring structure is not horizontally uniform due to wakes or waves, phys depends on
the azimuthal angle of the observer. Usually, phys takes its minima when the line of
sight is aligned to the wake direction (e.g., Colwell et al. 2006; Hedman et al. 2007).
Fig 5.3 illustrates this principle.
We developed a ray-tracing code to derive phys using outputs of N -body sim-
ulations. The basic schemes employed in the code are similar to those given in Salo
and Karjalainen (2003). Brie
y, the algorithm is as follows:
 The simulated region is partitioned into two-dimensional sub-cells, with par-
ticles allocated by sub-cell.
 A ray is generated by randomly choosing x and y coordinates in the ring
midplane, where z=0, and forming the unit vector n directed towards the
light source. In the Hill coordinate system, n is dened by its components:
nx = cos s sinBs;
ny = sin s sinBs;
nz = cosBs;
184
Figure 5.3: Schematic representation of azimuthal brightness asymme-
try in the ring system due to the presence of wake structure. At low
elevation angles, B, the wakes, which trail here by about 21 with re-
spect to the observer's tangential direction, are seen along their long axis
at ring longitudes, , of 69 and 249. Therefore, the rareed regions are
visible, reducing the total re
ecting surface area. This corresponds to
minimum brightness (top plot) or minimum optical depth (bottom plot).
At longitudes 90 away, the rareed regions are hidden, maximizing the
total re
ecting surface area. This corresponds to maximum brightness
levels and optical depths. Source: Schmidt et al. (2009)
185
where s and Bs are the azimuthal and elevation angle of the light source.
 We rst check whether the ray coming from the light source intersects any
sub-cells, and then check whether the ray intersects any particles within an
intersected sub-cell. If multiple sub-cells can be intersected by a single ray,
then these sub-cells are checked in the order of which is intersected rst by
the ray's trajectory.
 If the ray passes through the simulated region without intersecting any parti-
cles, it may still intersect with a ghost particle. Therefore, if an intersection is
still possible, the ray is regenerated with a new position given by the periodic
boundary conditions of the shearing system. Intersections are still possible if
pz < max(zi + ri) (5.8)
or
pz > min(zi  ri) (5.9)
where pz is the z coordinate of the ray, zi is the z coordinate of any particle
in the simulation, and ri is the corresponding particle radius. Therefore, the
algorithm compares the z coordinate of the ray to the maximum and minimum
vertical extent of the simulated particles.
 We continue generating randomly positioned rays until we numerically con-
verge on a value for the physical optical depth. The simulated physical optical
depth is found by setting I0, in Eq. 5.7, to the number of rays produced, and
186
I to the number of rays that pass through the ring without intersecting any
particles.
5.4 Results
We have performed simulations that explore the eects of particle density
and friction on the collisional dynamics of the ring particles and the subsequent
large-scale structure evolution. Each simulation was run for at least 100 orbits
in order for the patch to safely reach collisional equilibrium. Since the surface
properties of actual ring particles are not well constrained, we explored a wide range
of friction properties. This allowed us to gauge the range of possible behaviors of
large-scale structure formation. Furthermore, by matching collisional outcomes in
our simulations to observations, we may be able to constrain the frictional (surface)
properties of real ring particles.
Table 5.1{5.1 lists the simulations that were performed. Here we also sum-
marize the outcome of each simulation by providing the equilibrium collisional and
translational viscosities. Furthermore, we also give the Toomre Q of the ring, where
Q is dened as:
Q =
cr

G
; (5.10)
where cr is the radial velocity dispersion, and  is the surface density. Q attains a
time-averaged equilibrium value of the order of 2 in the case of strong wake structure.
Previous results have shown the emergence of strong wake structure when the radial
187
velocity dispersion due to collisions and gravitational encounters does not exceed
that corresponding to Q  2 (Salo, 1995; Salo et al., 2001; Ohtsuki & Emori, 2000).
Table 5.1: Summary of ring patch parameters and simulations results.  is the ring
particle internal density, rp is the ring particle radius, r is the coecient of static
friction, r is the coecient of rolling friction, coll is the collisional viscosity (in
units of 
0r
2
p) , trans is the translational viscosity (in units of 
0r
2
p), and Q is the
Toomre parameter.
 (g cm3) rp (m) s r coll trans Q
0.45 1.0 0.0 0.0 5.70 8.29 1.34
0.45 1.0 0.0 0.2 5.64 8.23 1.31
0.45 1.0 0.2 0.0 8.88 14.12 1.86
0.45 1.0 0.2 0.2 9.91 16.41 2.09
0.45 1.0 0.4 0.0 9.24 16.32 2.07
0.45 1.0 0.4 0.2 18.89 30.81 2.65
0.45 1.0 0.6 0.0 13.98 20.98 2.30
0.45 1.0 0.6 0.2 17.81 28.21 2.58
0.45 1.0 0.8 0.0 13.12 21.06 2.43
0.45 1.0 0.8 0.2 20.75 40.96 3.44
0.45 1.0 1.0 0.0 12.97 20.08 2.04
0.45 1.0 1.0 0.2 20.35 40.86 3.22
0.60 0.75 0.0 0.0 6.81 14.41 2.34
0.60 0.75 0.0 0.2 6.98 15.65 2.25
0.60 0.75 0.2 0.0 10.97 22.43 2.53
0.60 0.75 0.2 0.2 10.71 20.04 2.71
0.60 0.75 0.4 0.0 10.96 21.10 2.48
0.60 0.75 0.4 0.2 13.34 24.81 2.96
0.60 0.75 0.6 0.0 7.80 15.07 3.05
0.60 0.75 0.6 0.2 16.67 34.48 3.09
0.60 0.75 0.8 0.0 10.51 20.48 2.67
0.60 0.75 0.8 0.2 14.57 25.74 3.13
0.60 0.75 1.0 0.2 11.82 21.10 3.32
188
Figure 5.4: Snapshots of some of our simulations after 100 orbits of col-
lisional and gravitational evolution. The larger images show the position
of ring particles in radial (x) and azimuthal (y) space. The smaller im-
age shows the positions of particles in radial and vertical (z) space. The
left and right columns show ring particles with internal densities of 0.45
and 0.9 g/cm3, respectively. When we allow the particles to interact
with higher frictional forces (lower rows in the gure have higher fric-
tion), we nd that the ring particles are better at maintaining azimuthal
(y-direction) axisymmetric features, brought about by the viscous over-
stability. For lower friction, the patch is not viscous enough to be stable
against self-gravity wakes (inclined features most prominent in top pan-
els).
189
Table 5.1: Summary of ring patch parameters and simulations results, continued.
 (g cm3) rp (m) s r coll trans Q
0.75 0.60 0.0 0.0 7.45 22.19 2.60
0.75 0.60 0.0 0.2 7.45 20.62 2.73
0.75 0.60 0.2 0.0 12.01 26.76 2.95
0.75 0.60 0.2 0.2 14.03 28.32 3.14
0.75 0.60 0.4 0.0 13.43 30.97 2.92
0.75 0.60 0.4 0.2 16.52 32.17 3.22
0.75 0.60 0.6 0.0 13.47 27.77 2.85
0.75 0.60 0.6 0.2 17.02 32.20 2.99
0.75 0.60 0.8 0.0 13.87 30.43 2.99
0.75 0.60 0.8 0.2 17.68 32.26 2.93
0.75 0.60 1.0 0.0 14.29 31.33 2.70
0.75 0.60 1.0 0.2 16.06 28.80 2.91
0.90 0.50 0.0 0.0 7.45 22.19 1.816
0.90 0.50 0.0 0.2 7.45 20.62 2.155
0.90 0.50 0.2 0.0 12.01 26.76 2.49
0.90 0.50 0.2 0.2 14.03 28.32 2.59
0.90 0.50 0.4 0.0 13.43 30.97 2.94
0.90 0.50 0.4 0.2 16.52 32.17 3.27
0.90 0.50 0.6 0.0 13.47 27.77 3.86
0.90 0.50 0.6 0.2 17.02 32.20 3.74
0.90 0.50 0.8 0.0 13.87 30.43 3.73
0.90 0.50 0.8 0.2 17.68 32.26 4.01
0.90 0.50 1.0 0.0 14.29 31.33 4.24
0.90 0.50 1.0 0.2 16.06 28.80 4.33
5.4.1 Correlation between Friction and Viscosity
Figure 5.4 shows snapshots of our simulations after 100 orbits of collisional
and gravitational evolution. We nd that when particles have higher surface fric-
tion, axisymmetric features become more prevalent. Visual inspection of multiple
frames allows us to ascertain that these features are likely brought about by viscous
overstability wakes. At lower friction values, the ring particles form self-gravity wake
190
structure that are inclined with respect to the orbital direction (inclined features
most prominent in top panels of Figure 5.4).
The viscous overstability comes about by the combined eects of inelastic
collisions and Keplerian shear, which lead to an energetic steady state. In this
state, inelastic collisions cause overdense regions to develop. Due to energy input
from the Keplerian shear, the particle 
ux is directed away from regions of high
density. The mass in over-dense regions 
ows to under-dense regions; however,
this mass 
ow overshoots due to the momentum of the particle 
ow, leading to
new overdense regions. Studies show that this overstability develops for regions in
the ring where the viscosity increases with optical depth beyond some criticial rate
(Salo et al., 2001; Daisaka et al., 2001; Rein & Latter, 2013):
 =
d ln 
d ln 
> cr  1
where  is the viscosity, and  is the optical depth.
For moderately high-friction cases, we see that the patch also reaches an equi-
librium state with characteristic viscous overstability pulsations. Unlike previous
simulations that studied the overstability (Salo et al., 2001), we see very little sign
of gravitational wakes in this case even for a moderately high internal particle den-
sity of 0.45 g cm3. Salo et al. (2001) showed that the viscous overstability is present
when ring particles have low internal density (0.20{0.30 g cm3), and may be slightly
featured (less prominent axisymmetric features) at intermediate internal densities
(0.45 g cm3), and likely non-existent for particles with the density of actual water
ice 0.90 g cm3. We nd that, for suciently high surface friction, overstability
191
features may exist even for non-porous water ice with densities of 0.90 g cm3 (see
bottom-right panel of Figure 5.4).
Quantitatively, we can compare the change in the collisional viscosity of the
patch with the Toomre Q parameter, as a function of particle static friction s
(Figure 5.5). The Q parameter is a ratio of the particle's radial velocity dispersion
to the ring patch's surface density. Self gravity wakes develop for Q < 2. A high
friction value increases the viscosity of the patch, and inhibits the formation of
self-gravity wakes.
5.4.2 Correlation with Density
We studied the dependence of the collisional viscosity with density. We nd
that for suciently large friction values (here s = 1), the viscosity stabilizes the
disk against self-gravity wakes even for very dense particles. Previous results that
did not include explicit treatment of friction and multi-contact physics (Salo et al.,
2001; Hedman et al., 2014) have shown that the internal density of particles needs
to be < 0.30 g cm3 for viscous overstability to develop. Here we nd that even
particles with 0.90 g cm3 can exhibit the overstability.
To further test the in
uence of particle density on large-scale structure for-
mation, we plotted the relationship between , the friction parameters, and Q in
Figure 5.6. We show the value of Q found in simulations with three dierent par-
ticle densities, for values of s ranging from 0.0{1.0, and for a rolling friction value
of r = 0.2. Here, we see that Q rises most steeply with s for the lowest particle
192
Figure 5.5: Quantitatively, we can compare the change in the collisional
viscosity (left axis, solid black curve) of the patch with the Toomre Q
parameter (right axis, dashed blue curve), as a function of particle static
friction s. A high friction value increases the viscosity of the patch, and
inhibits the formation of self-gravity wakes.
193
Figure 5.6: The relationship between particle density (denoted by dier-
ent linestyles and colors), the static friction component, and Toomre Q.
These are cases where r = 0.2. The solid cyan line represents cases with
 = 0.45 g cm3 and rp = 1.0 m. The dashed magenta line represents
cases with  = 0.60 g cm3 and rp = 0.75 m. The dotted black line
represents cases with  = 0.75 g cm3 and rp = 0.60 m.
density. The friction parameter here acts to increase the viscosity of the patch, and
subsequently drives the formation of large-scale structure and overstability wakes.
At suciently high friction values, the value of Q seems to saturate at  3.
5.4.3 Physical Optical Depth Measurements
Using the ray-tracing algorithm described in Section 2.3.2, we can generate
synthetic optical depth proles of our library of simulations, in order to compare
with Cassini observations.
194
By varying the azimuthal observing angle, , we are able to investigate the
nature of the large-scale structure in our simulations. Structures generated by the
viscous overstability are axisymmetric; therefore, the physical optical depth, phys,
is at a minimum when  is aligned parallel to the azimuthal direction and at a
maximum when  is aligned parallel to the radial direction.
In Fig. 5.7, we show our generated phys proles for two dierent particle den-
sity cases, 0.45 and 0.60 g cm3. In these gures,  = 0 corresponds to an observing
angle aligned with the radial direction facing inwards towards Saturn.  increases
in a counterclockwise fashion from this initial geometry. For each simulation, phys
is calculated for an elevation angle of B = 30, and for  ranging from 0 to 2
sampled at 4 intervals. This value of B is chosen to match that of the Cassini
spacecraft. Furthermore, intermediate values of B results in variations phys with .
As B approaches 90, the viewing geometry will be less in
uenced by changes in ,
resulting in a near-constant phys. For both particle density cases, we see that the
amplitude of the phys variation decreases as a function of the friction parameters.
This suggests that friction acts to create a larger dispersion in the physical location
of particles in the patch. This is consistent with the idea that friction acts against
self-gravitational clumping. For most of the cases seen here for both particle densi-
ties, we see that the minima are located at  values between  90 and 110. This is
consistent with observations of  90 to 120 tilts in the structure found in Saturn's
dense rings (e.g. Colwell et al. 2007). The typical value is 90 for the outer B ring,
indicating the existence of overstablity features. Although the typical pitch angle
seen in all of our simulations is consistent with observations, the average values of
195
phys we measure ( 1) are much lower than those for the outer B ring, which are
typically a few factors larger than 1 and can be greater than 5. It is possible that
inter-wake gaps are lled by small particles that are ignored in our simulations, but
contribute most of the observed opacity.
We also nd that for higher friction values, the minima are shifted towards 90,
or axisymmetry. This is true for both density cases. In order to obtain an accurate
estimate of the azimuthal angle that coincides with a given minimum phys, a sinu-
soidal curve that best ts the data was generated using a least-squares minimization
technique. For the case of 0.45 g cm3, we nd minima at values of 100.1, 93.6,
and 88.5 for s = 0.6, 0.8, and 1.0, respectively. For the case of 0.60 g cm3, we nd
minima at values of 109.5, 103.8, and 95.8 for s = 0.6, 0.8, and 1.0, respectively.
The shift is more apparent in the top panel of Fig. 5.7, the lower-density case. This
is expected, as larger particle densities tend to lead to stronger self-gravitational
wakes. Hence, for similar friction values, lower-density particles are more strongly
driven to axisymmetry and viscous overstability.
5.4.4 Comparison with Cassini UVIS occultations
Finally, we attempt to compare the results of our simulations to Cassini UVIS
stellar occultations (Colwell et al., 2006). Observations of stellar occultations by the
UVIS instrument reveal the transparency of the ring as a function of observation
geometry. These observations vary both the azimuthal observing angle, , and the
spacecraft elevation angle B. Therefore, in order to eectively compare our simu-
196
Azimuthal Observing Angle − θ
0.92
0.94
0.96
0.98
1.00
1.02
N
o
rm
a
li
ze
d
 A
v
er
a
ge
 τ
p
h
y
s
0.45 g/cm3
0 50 100 150 200 250 300 350
Azimuthal Observing Angle − θ
0.85
0.90
0.95
1.00
N
or
m
al
iz
ed
 A
v
er
ag
e 
τ p
h
y
s
0.60 g/cm3
µs =0.6
µs =0.8
µs =1.0
Figure 5.7: The normalized average physical optical depth as a function
of observing angle () at an elevation B = 30 for simulations with a
particle density of 0.45 g cm3 (top panel) and 0.60 g cm3 (bottom
panel). For both panels, the solid, dash-dotted, and dotted curves rep-
resent simulations where particles have static friction values of 0.6, 0.8,
and 1.0, respectively. Individual data points are plotted on the solid line
in the top panel to show the 4 sampling frequency that was used for each
curve. An increase in the particle surface roughness evidently causes an
observational shift from tilted self-gravity wake features to azimuthally
coherent viscous overstability features.
197
lation data to observations, we calculate the average transparency of our simulated
patch for unique combinations of  and B. The transparency, T , of the ring patch
is dened as
T = ephys=; (5.11)
where  is the projection factor,  = j sin(B)j. In Fig. 5.8, we plot T as a function
of  for two particle densities and two dierent values of static friction s. For
each value of , we plot results for a range of the azimuthal observing angles, . As
expected, for larger elevation angles, B, the range of transparencies becomes smaller
across the range of . This is due to the fact that at high elevation angles, variation
in  would result in very minimal changes to the observing geometry. Overall, we
nd that a simulation where particles have lower densities and higher friction ( =
0.45 g/cm3 and s = 1.0) is the best t to the observations. Furthermore, in general,
a higher value of static friction results in lower transparencies.
5.5 Conclusions & Future Work
We presented ring simulations that reveal the relationship between particle
friction and large-scale structure in Saturn's B ring. A comprehensive explanation
of the interplay between the mechanisms of self-gravity and viscous overstability has
remained elusive. It is understood that both play an important role in generating
large-scale structure in Saturn's rings. Here we have shown that it is possible that
the surface properties of particles, here modeled through two friction parameters
(static and rolling friction), may control whether a patch of ring particles exhibits
198
0.1
0.2
0.3
0.4
0.5
T
ra
n
sp
a
re
n
cy
 (
e
−τ
/µ
)
0. 45 g/cm3Cassini UVIS data
µs =1. 0
µs =0. 6
0.2 0.4 0.6 0.8 1.0
µ=sin(B)
0.0
0.1
0.2
0.3
0.4
0.5
T
ra
n
sp
ar
en
cy
 (
e
−τ
/µ
)
0. 60 g/cm3
Figure 5.8: The transparency of the simulated data (circles and squares)
is compared to Cassini UVIS stellar occultation data (crosses). The top
and bottom panels plot transparency of two friction cases (cyan circles
for s = 1:0 and magenta squares for s = 0.6) as a function of the
projection factor for simulations with particle densities of 0.45 g/cm3
(top panel) and 0.60 g/cm3 (bottom panel), respectively. We nd that
higher values of friction lead to lower ring patch transparency.
199
self-gravitational wakes or axisymmetric overstable features, or both.
We nd that higher interparticle friction increases the viscosity (collisional and
translational), such that the ring particles become stable against self-gravitational
perturbations. Furthermore, we have shown that the presence of axisymmetric fea-
tures via viscous overstability does not necessarily require highly porous particles. In
fact, non-porous particles may generate viscous overstability features if their surface
roughness is suciently high.
This has important implications for the possible mass of the B ring itself.
Previous studies (Robbins et al., 2010) have relied on similar simulations to estimate
the masses of the A and B rings. Robbins et al. (2010) suggested an upper limit for
ring particle density of 0.45 g/cm3, based on comparisons of simulation data to ring
opacity measurements. However, they did not include the eects of inter-particle
friction. If highly frictional particles exist in the B ring, then axisymmetric features
in the B ring may be produced by particles with densities as high as 0.90 g/cm3,
and current values for the mass of the B ring may be underestimated by at least a
factor of 2.
We have performed comparisons to Cassini UVIS stellar occultation data and
found that particles with lower densities and higher friction produce ring transparen-
cies that best match the observations. This suggests that axisymmetric features are
likely prevalent in the B ring, as both conditions of low densities and high fric-
tion are conducive to the generation of viscous overstability. However, the average
physical optical depths measured for these simulations ( 1) are much lower than
those observed in the B ring (> 4 in some locations). If we were to simulate a ring
200
patch with a larger size distribution of particles (and with more small particles in
particular), the transparencies would likely be lower for all our cases. Simulations
by Robbins et al. (2010) showed that for high particle densities (> 0.45 g/cm3), a
particle size distribution does lead to more opaque rings. However, this eect was
not as pronounced for simulations with lower particle densities. Nevertheless, this
suggests a possibility that particles with densities greater than 0.45 g/cm3 may still
be able to produce features that are suciently opaque to match observations, with
the caveat that they would need to have high surface roughness.
Since friction plays an important role in controlling the macroscopic properties
of the ring, our simulations have shown that understanding the surface roughness of
the ring particles is essential if observations of large-scale structures are to be used
to interpret their physical properties (size, densities, etc.). Here, we have chosen
to explore a wide range of the friction parameter space in order to understand the
range of possible behaviors. The actual surface roughness of ice particles in a vacuum
are unknown. So far, simulations have relied on experimental work (Bridges et al.,
1984) that was performed on perfectly smooth and non-rotating ice spheres. There
have been some experiments that have measured the friction properties of ice on
ice. Sukhorukov & Lsetl (2013) conducted eld tests on the friction of sea ice on
sea ice in the Bahrents sea in temperatures of 253 to 271 K, and found that the
coecient of static friction can range from  0.6 to 1.26 depending on the time the
ice is held together before sliding is allowed to begin (higher hold times result in
higher friction). Schulson & Fortt (2012) measure the coecient of kinetic friction
of freshwater ice sliding slowly (with speeds up to 1 mm/s) at temperatures of 98 to
201
263 K. They found that the kinetic friction can range from 0.15 to 0.76 depending
on the temperature and sliding speed. While these results can provide some ground
truths to numerical studies, the environmental conditions they were conducted in
do not match that found around Saturn (vacuum at temperatures that range from
20 to 40 K (Flasar et al., 2005)) . Therefore, future experimental and theoretical
studies should explore the collisional interaction of realistic ice particles in order to
place better constraints on their surface properties.
Future work should consider a size distribution of particles (within a single
simulation) and a wider range of initial geometric optical depths in order to better
quantify the eect that particle surface roughness has on observable properties of the
rings. Furthermore, recent studies (Hedman et al., 2014) have shown that axisym-
metric features in the B ring are azimuthally coherent in length scales much larger
than those simulated here. If inter-particle friction does indeed allow high-density
particles to form overstabilities, then simulation with much larger azimuthal extents
will have to be performed to show that such developed structures are stable in the
long-term. Finally, we intend to study the eect of inter-particle cohesive forces on
the formation of km-scale structure. Tremaine (2003) showed that cohesion can have
an important role in creating structure in Saturn's ring and predicted the formation
of shear-free ringlets. We will study the extent to which cohesive forces can enhance
or diminish the generation of self-gravitational wakes and axisymmetric features, or
if they generate separate unique large-scale features.
202
Chapter 6: Conclusions & Future Work
6.1 Conclusions
Our study of granular physics in the Solar System has shown that granular
matter behaves quite dierently in low-gravity environments. In particular, we nd
surprising that the manner in which grains 
ow past each other is dierent for small
solar system bodies compared to Earth. In Chapter 2, we conducted a study of
vibration-driven size-segregation of grains to study the eectiveness of the Brazil-
nut eect (BNE) in dierent boundary and gravity conditions that closely resemble
those found on asteroids. While previous studies had shown that the classical BNE
eect can occur in low-gravity environments, the studies were limited to containers
with lateral boundaries. We were able to better mimic the asteroid environment
by removing these lateral boundaries, using periodic boundary conditions (PBC).
We were able to show that for a wide range of friction properties, the BNE occurs
even when PBC are used. However, rather than granular convection, a void-lling
mechanism drives the BNE when PBC are present. By analyzing the 
ow patterns of
grains during the oscillatory driven motion, we were able to show that the vibration
causes the granular system to dilate, creating void spaces in the medium. Due to the
size dierence between the large intruder and the surrounding grains, the smaller
203
grains have a higher probability of lling the void space created below the intruder.
This allows the BNE to occur in a system without conning wall boundaries, which
are otherwise required for granular convection to occur. Furthermore, we nd hints
that the void-lling mechanism may be more ecient in lower-gravity environments.
Our study suggests that the BNE driven by a void-lling mechanism can occur on
an asteroid, and it may be able to explain the presence of large boulders on their
surfaces.
In Chapter 3, we studied low-speed impacts onto granular material in low-
gravity environments in the context of the OSIRIS-REx mission to the asteroid
Bennu. We demonstrated that the outcome of a sampling attempt by the TAGSAM
varies depending on the properties of the regolith. By generating an atlas of simula-
tion results, we will be able to develop maps of the sampleability of dierent regions
of the surface of Bennu. When touchdown does occur in 2020, we will be able to
reconstruct the response of the OSIRIS-REx spacecraft as it touches down on the
surface. Furthermore, by comparing our synthetic accelerometer data to the actual
behavior of the spacecraft, we may be able to derive the properties of the surface
regolith.
We showed indications that granular material behaves in a more 
uid-like
manner in a low-gravity environment, and does not resist penetration as much as
it would in 1g. This was also demonstrated through low-speed impact experiments
of an aluminum sphere onto dry sand. Our analysis shows that, compared to low-
speed impacts on Earth, an impactor experiences less drag as it penetrates through
a granular medium on an asteroid. The simulations suggest that the drag coecient
204
varies as log

g


. Further simulations and experiments are required to fully develop
a coherent theory of granular impacts in low-gravity environments. Calibrating our
code to low-gravity experiments will be an important next step to further elucidating
such a theory. An important breakthrough in this process will be the development
of an experimental apparatus that is able to replicate the gravity conditions on an
asteroid. Some low-gravity experiments have been performed in Atwood-type ma-
chines (Altshuler et al., 2014), but the realistically achievable gravity level in these
experiments are still orders of magnitude greater than that found on an asteroid.
Developing a coherent theory of granular impacts in low gravity enables us to better
predict the outcome of a robotic or human interaction on the surface of an asteroid.
In Chapter 4, we moved from discussions of cm-size particles on the surface or
near-surface of small bodies to exploring the collisions of rubble piles. We studied
the eect of initial rotation on the outcome of rubble-pile collisions by analyzing the
properties of the largest remnant and material that is gravitationally bound to it. We
found that the collisional mass loss is sensitive to the initial rotation period, with
faster-rotating rubble-pile targets dispersing more mass. The mass of the largest
remnant of head-on impacts is well described by the \universal law" for catastrophic
disruption rst put forward by Leinhardt & Stewart (2009), independent of initial
pre-impact rotation. However, the actual interacting mass fraction of a rubble pile
needs to be considered for oblique impacts to match the universal law. Furthermore,
we developed a prescription that describes the change in the catastrophic disruption
criteria of head-on equatorial impacts onto a rotating target. Finally, we found that
dissipation and frictional eects combine to increase the catastrophic disruption
205
threshold in rubble-pile collisions. A logical next step for this study is to apply our
revised mass-loss prescription to a larger study of proto-planetary disk evolution,
where we can test whether the eects of pre-impact rotation drastically in
uence
the growth of large planetesimals.
In Chapter 5, we studied large-scale structure in Saturn's B ring. The physical
properties of the icy particles that make up Saturn's rings are poorly constrained,
since even the highest-resolution Cassini images of the rings are unable to resolve
individual ring particles. These properties must be inferred from the observable
structures in the rings that are generated from the collective collisional, gravita-
tional, and possibly electrostatic interactions of ring particles. Two types of km-
scale macroscopic structures have been observed in the B ring: self-gravity wakes,
which are tilted with respect to the azimuthal direction of the ring system, and
axisymmetric viscous overstability structures. A comprehensive explanation of the
interplay between the two mechanisms of self-gravity and viscous overstability has
remained elusive. It is understood that both play an important role in generating
large-scale structure in Saturn's rings. We showed that the surface properties of
particles controls whether a patch of ring particles exhibits self-gravitational wakes
or axisymmetric overstable features, or both. We nd that higher inter-particle
friction increases the viscosity (collisional and translational), such that the ring par-
ticles become stable against self-gravitational perturbations. This allows particles
with relatively high densities to form overstable features. This has important impli-
cations for the possible mass of the B ring, as previous studies have suggested that
the presence of overstability features precludes the existence of non-porous particles
206
in the B ring. By comparing synthetic observational data generated from analyzing
simulation output with Cassini UVIS stellar occultation data, we nd that parti-
cles with lower densities and higher friction produce ring transparencies that best
match the observations. However, the average physical optical depths measured for
these simulations ( 1) are only found in the innermost B1 portion of the B ring.
Therefore, we expect that simulations with a larger size distribution of particles, or
a larger optical depth, are necessary if we wish to study the higher-opacity regions
of the B ring.
6.2 Future Work
In this section, we outline work that will build upon the results presented in
this thesis. First, we discuss on-going work to study the interaction of the cube-
shaped Mobile Asteroid Surface Scout (MASCOT) that will be deployed by the
Hayabusa 2 mission to the asteroid Ryugu. Then, we discuss work that we have
proposed to do on studying the plateaus in Saturn's C ring. These two applications
further show the scope of possible problems that can be approached with numerical
simulations of granular material. With each application, code improvements are
required to enable the study of that application and to ensure that the software's
performance and capabilities does not stagnate. Future possible improvements in-
clude a full treatment of polyhedral particles to model irregularly shaped granular
materials. Such an implementation would enable even more applications of the code,
such as studies of situations where grains may fracture or cases where the shape of
207
the particles are well constrained.
6.2.1 MASCOT lander
The Hayabusa 2 mission, a direct successor of the Hayabusa mission, is the
second asteroid-sample-return mission of the Japan Aerospace Exploration Agency
(JAXA) (Tsuda et al., 2013). It is expected to reach the near-Earth asteroid Ryugu
in July 2018. Hayabusa 2 will deploy a cube-shaped lander, MASCOT, that is
expected to bounce and eventually settle on the asteroid surface. The purpose of
the lander is to characterize the surface features of Ryugu using four instruments:
a camera, a hyperspectral microscopic imager, a magnetometer and a radiometer.
By studying the lander's mechanical response to the surface, we may also be able
to determine the mechanical properties of the asteroid's regolith.
In collaboration with instrument designers at the German aerospace center
(DLR) and the French space agency (CNES), we have been performing simulations
in order to characterize the bouncing of the lander as a function of assumed re-
golith properties, in the gravitational environment of the target asteroid Ryugu (see
Fig. 6.1). A key goal of this project is to determine the eective coecient of resti-
tution (eCOR) of a lander impacting the surface of an asteroid, as a function of
both regolith properties and initial impact conditions (angle of approach, rotation
velocity, impact speed, etc.). The eective coecient of restitution is a measure of
the ratio of the nal linear and rotational impact to their initial values. Laboratory
experiments by MASCOT designers have shown that the value of the eCOR can vary
208
quite wildly. For very special circumstances it can even be greater than unity. Our
simulations have attempted to constrain the range of possible values, and attempt to
nd correlations with the initial conditions of the lander. Using simulations has the
added benet of testing whether a very-low-gravity environment aects the value of
eCOR, compared to an Earth-gravity experiment. Preliminary results have shown
that the eCOR is sensitive to the initial orientation of the lander.
Understanding how the eCOR can vary is important for ensuring spacecraft
safety. By estimating the eCOR for an individual bounce, we can then predict
subsequent impact conditions and formulate the spacecraft's future trajectory. Fur-
thermore, we can measure the response of the regolith to the bouncing lander. The
mass and speed of grains ejected by the impact depends on the impact conditions.
If we are able to know the pre-impact state of MASCOT, we can measure the prop-
erties of the ejecta created in our simulations. The mass and speed of the ejecta
also depends on the properties of the regolith. By generating an atlas of simulation
results, and comparing with spacecraft observations of an impact event, we can then
make some determination of the regolith properties.
6.2.2 The Plateaus in Saturn's C ring
We have shown that numerical simulations are a crucial tool for the study
of large-scale structure in Saturn's rings. By recreating the conditions necessary
for the formation of these macroscopic structures, the microscopic properties of
the constituent ring particles can be inferred. Saturn's C ring exhibits peculiar
209
Figure 6.1: Snapshots at three dierent times (ordered by column) of
four dierent simulations (ordered by row). For each row, MASCOT has
a unique pre-impact orientation. From top to bottom these are: a 
at
orientation, where the lander impacts the surface on a 
at face; edge-
on-parallel, where the lander impacts the surface along an edge that
is oriented parallel to the direction of motion; edge-on-perpendicular,
where the edge is oriented perpendicular to the direction of motion; and
a case where the corner of the lander touches the surface rst. The post-
collisional properties of MASCOT are sensitive to its initial orientation.
210
Figure 6.2: Radial prole of Saturn's C ring from UVIS stellar oc-
cultation data. High optical depth regions (plateaus) are interspersed
between lower optical depth regions. Source: Colwell et al. (2009).
structures, whose presence has eluded explanation. Figure 6.2 Colwell et al. (2009)
shows the optical depth prole of the C ring, which displays broad regions of optical
depth less than 0.1, interspersed with narrow \plateaus" (labeled P1 through P11 in
the gure) with optical depths that are 3 to 5 times the optical depth of most of the
C ring. Hedman & Nicholson (2014) have found that the surface mass densities of
the plateaus are very similar to that of the surrounding regions, in spite of the very
dierent optical depths. This implies that either the average particle size is much
smaller in the plateaus or that the interior density of the particles is substantially
lower in the plateaus.
211
Furthermore, high resolution-Cassini images of plateaus Spitale & Porco (2007)
found narrow 10-km-long striations aligned with the orbital direction, which are sim-
ilar to axisymmetric viscous overstability structures. It is possible that these features
are a result of jamming transitions where the particles are temporarily locked up
into narrow, low-shear zones where the thickness of the rings can locally increase
due to particles moving vertically out of the ring plane. We intend to study this
phenomenon by rst studying the jamming transition more intently. Published work
on granular 
ows shows that the jamming transition only occurs for specic combi-
nations of particle density, friction and particle lling factors (Otsuki & Hayakawa,
2011; Grob et al., 2014). By rst determining these properties in laboratory-scale
simulations, we will then be able to see if the jamming transition can occur in the
context of a ring particle simulation. Reproducing the jamming transition and the
associated gaps and ridges in the C ring plateaus, as well as the observed average
optical depth and surface mass density, would allow us to place tight constraints on
the particle properties in this part of Saturn's rings.
212
Bibliography
Abell, P. A., Barbee, B. W., Chodas, P. W., et al. 2015, in Asteroids IV, ed.
P. Michel, F. E. DeMeo, & W. F. Bottke (University of Arizona Press), 855{
880
Altshuler, E., Torres, H., Gonzalez-Pita, A., et al. 2014, Geophysical Research Let-
ters, 41, 3032
Asphaug, E. 2007, in Lunar and Planetary Science Conference, Vol. 38, Lunar and
Planetary Science Conference, 2432
Asphaug, E. 2010, Chemie der Erde / Geochemistry, 70, 199
Asphaug, E., King, P. J., Swift, M. R., & Merrield, M. R. 2001, in Lunar and
Planetary Science Conference, Vol. 32, Lunar and Planetary Science Conference
Barabasi, A.-L., Albert, R., & Schier, P. 1999, Physica A Statistical Mechanics
and its Applications, 266, 366
Barnouin-Jha, O. S., Cheng, A. F., Mukai, T., et al. 2008, Icarus, 198, 108
Benz, W., & Asphaug, E. 1999, Icarus, 142, 5
Berry, K., Sutter, B., May, A., et al. 2016, in Annual AAS Guidance and Control
Conference, Vol. 36, AAS/Division for Planetary Sciences Meeting Abstracts
Biele, J., Ulamec, S., Maibaum, M., et al. 2015, Science, 349
Binzel, R. P., Morbidelli, A., Merouane, S., et al. 2010, Nature, 463, 331
Borderies, N., Goldreich, P., & Tremaine, S. 1985, Icarus, 63, 406
Bottke, W. F., Broz, M., O'Brien, D. P., et al. 2015, in Asteroids IV, ed. P. Michel,
F. E. DeMeo, & W. F. Bottke (University of Arizona Press), 701{724
Bottke, W. F., Durda, D. D., Nesvorny, D., et al. 2005, Icarus, 179, 63
Bottke, Jr., W. F., Cellino, A., Paolicchi, P., & Binzel, R. P., eds. 2002, Asteroids
III (University of Arizona Press)
Bottke, Jr., W. F., Vokrouhlicky, D., Rubincam, D. P., & Nesvorny, D. 2006, Annual
Review of Earth and Planetary Sciences, 34, 157
Bridges, F. G., Hatzes, A., & Lin, D. N. C. 1984, Nature, 309, 333
Canup, R. M. 2008, Icarus, 196, 518
Cellino, A., Zappala, V., Davis, D. R., Farinella, P., & Paolicchi, P. 1990, Icarus,
87, 391
Chapman, C. R. 1976, Geochim. Cosmochim. Acta, 40, 701
Cheng, A. F., Izenberg, N., Chapman, C. R., & Zuber, M. T. 2002, Meteoritics and
Planetary Science, 37, 1095
Chesley, S. R., Ostro, S. J., Vokrouhlicky, D., et al. 2003, Science, 302, 1739
Colwell, J. E. 2003, Icarus, 164, 188
Colwell, J. E., Cooney, J. H., Esposito, L. W., & Sremcevic, M. 2009, Icarus, 200,
213
574
Colwell, J. E., Esposito, L. W., & Sremcevic, M. 2006, Geophysical Research Letters,
33, L07201
Colwell, J. E., Esposito, L. W., Sremcevic, M., Stewart, G. R., & McClintock, W. E.
2007, Icarus, 190, 127
Cotto-Figueroa, D., Statler, T. S., Richardson, D. C., & Tanga, P. 2015, The As-
trophysical Journal, 803, 25
Daisaka, H., & Ida, S. 1999, Earth, Planets, and Space, 51, 1195
Daisaka, H., Tanaka, H., & Ida, S. 2001, Icarus, 154, 296
Delbo', M., dell'Oro, A., Harris, A. W., Mottola, S., & Mueller, M. 2007, Icarus,
190, 236
Delbo, M., Mueller, M., Emery, J. P., Rozitis, B., & Capria, M. T. 2015, in Asteroids
IV, ed. P. Michel, F. E. DeMeo, & W. F. Bottke (University of Arizona Press),
107{128
Delbo, M., Libourel, G., Wilkerson, J., et al. 2014, Nature, 508, 233
Dobrovolskis, A. R., & Burns, J. A. 1984, Icarus, 57, 464
Dollfus, A., Geake, J. E., Mandeville, J. C., & Zellner, B. 1977, in IAU Colloq. 39:
Comets, Asteroids, Meteorites: Interrelations, Evolution and Origins, ed. A. H.
Delsemme, 243{251
Dunham, D. W., Farquhar, R. W., McAdams, J. V., et al. 2002, Icarus, 159, 433
Flasar, F. M., Achterberg, R. K., Conrath, B. J., et al. 2005, Science, 308, 975
Forterre, Y., & Pouliquen, O. 2008, Annual Review of Fluid Mechanics, 40, 1
Fujiwara, A., Kawaguchi, J., Yeomans, D. K., et al. 2006, Science, 312, 1330
Garcia, R. F., Murdoch, N., & Mimoun, D. 2015, Icarus, 253, 159
Gray, J. M. N. T., & Thornton, A. R. 2005, Proceedings of the Royal Society of
London Series A, 461, 1447
Grob, M., Heussinger, C., & Zippelius, A. 2014, Physical Review E, 89, 050201
Guttler, C., von Borstel, I., Schrapler, R., & Blum, J. 2013, Physical Review E, 87,
044201
Hedman, M. M., & Nicholson, P. D. 2014, Monthly Notices of the Royal Astronom-
ical Society, 444, 1369
|. 2016, Icarus, 279, 109
Hedman, M. M., Nicholson, P. D., & Salo, H. 2014, The Astronomical Journal, 148,
15
Hedman, M. M., Nicholson, P. D., Salo, H., et al. 2007, The Astronomical Journal,
133, 2624
Hill, G. W. 1878, American Journal of Mathematics, 1, 5
Holsapple, K. A. 1994, Planetary & Space Sciences, 42, 1067
|. 2016, private communication
Housen, K. R., & Holsapple, K. A. 1990, Icarus, 84, 226
|. 2011, Icarus, 211, 856
Housen, K. R., Wilkening, L. L., Chapman, C. R., & Greenberg, R. 1979, Icarus,
39, 317
Jaeger, H. M., Nagel, S. R., & Behringer, R. P. 1996, Reviews of Modern Physics,
68, 1259
214
Jutzi, M., Michel, P., Benz, W., & Richardson, D. C. 2010, Icarus, 207, 54
Kato, M., Iijima, Y.-I., Arakawa, M., et al. 1995, Icarus, 113, 423
Katsuragi, H., & Durian, D. J. 2007, Nature Physics, 3, 420
Knight, J. B., Jaeger, H. M., & Nagel, S. R. 1993, Physical Review Letters, 70, 3728
Korycansky, D. G., & Asphaug, E. 2009, Icarus, 204, 316
Kudrolli, A. 2004, Reports on Progress in Physics, 67, 209
Lambe, T. W., & Whitman, R. V. 1979, Soil Mechanics, SI edition (Wiley, New
York)
Lay, T., & Wallace, T. C. 1995, Modern Global Seismology (Academic Press, San
Diego)
Leinhardt, Z. M., & Richardson, D. C. 2002, Icarus, 159, 306
|. 2005, The Astrophysical Journal, 625, 427
Leinhardt, Z. M., Richardson, D. C., & Quinn, T. 2000, Icarus, 146, 133
Leinhardt, Z. M., & Stewart, S. T. 2009, Icarus, 199, 542
|. 2012, The Astrophysical Journal, 745, 79
Leinhardt, Z. M., Stewart, S. T., & Schultz, P. H. 2008, in The Solar System Beyond
Neptune, ed. M. A. Barucci, H. Boehnhardt, D. P. Cruikshank, A. Morbidelli, &
R. Dotson (University of Arizona Press), 195{211
Levison, H. F., Thommes, E., & Duncan, M. J. 2010, The Astronomical Journal,
139, 1297
Lissauer, J. J. 1993, ARAA, 31, 129
Love, S. G., & Ahrens, T. J. 1996, Icarus, 124, 141
Matsui, T., Waza, T., Kani, K., & Suzuki, S. 1982, Journal of Geophysical Research,
87, 10968
Matsumura, S., Richardson, D. C., Michel, P., Schwartz, S. R., & Ballouz, R.-L.
2014, Monthly Notices of the Royal Astronomical Society, 443, 3368
Merline, W. J., & Chapman, C. R. 2001, Meteoritics and Planetary Science Supple-
ment, 36, A132
Michel, P., Benz, W., & Richardson, D. C. 2004, Planetary & Space Sciences, 52,
1109
Michel, P., Benz, W., Tanga, P., & Richardson, D. C. 2001, Science, 294, 1696
Michel, P., DeMeo, F. E., & Bottke, W. F., eds. 2015a, Asteroids IV (University of
Arizona Press)
Michel, P., O'Brien, D. P., Abe, S., & Hirata, N. 2009, Icarus, 200, 503
Michel, P., & Richardson, D. C. 2013, Astronomy & Astrophysics, 554, L1
Michel, P., Richardson, D. C., Durda, D. D., Jutzi, M., & Asphaug, E. 2015b, in
Asteroids IV, ed. P. Michel, F. E. DeMeo, & W. F. Bottke (University of Arizona
Press), 341{354
Michikami, T., Nakamura, A. M., Hirata, N., et al. 2008, Earth, Planets, and Space,
60, 13
Miyamoto, H., Yano, H., Scheeres, D. J., et al. 2007, Science, 316, 1011
Morris, A. J. W., Burchell, M. J., Price, M. C., Wozniakiewicz, P. J., & Cole, M. J.
2012, in European Planetary Science Congress 2012, EPSC2012
Mukhopadhyay, S., Allen, B., & Brown, E. 2014, ArXiv e-prints, arXiv:1407.0719
Murdoch, N., Rozitis, B., Green, S. F., et al. 2013, Monthly Notices of the Royal
215
Astronomical Society, 433, 506
Murdoch, N., Sanchez, P., Schwartz, S. R., & Miyamoto, H. 2015, in Asteroids
IV, ed. P. Michel, F. E. DeMeo, & W. F. Bottke (University of Arizona Press),
767{792
Nakamura, A. M., Setoh, M., Wada, K., Yamashita, Y., & Sangen, K. 2013, Icarus,
223, 222
Nesvorny, D., Bottke, W. F., Vokrouhlicky, D., Chapman, C. R., & Rafkin, S. 2010,
Icarus, 209, 510
Ohtsuki, K., & Emori, H. 2000, The Astronomical Journal, 119, 403
Orlova, N. S. 1962, in IAU Symposium, Vol. 14, The Moon, ed. Z. Kopal & Z. K.
Mikhailov, 411{414
Otsuki, M., & Hayakawa, H. 2011, Physical Review E, 83, 051301
Perera, V., Jackson, A. P., Asphaug, E., & Ballouz, R.-L. 2016, Icarus, 278, 194
Perrine, R. P., & Richardson, D. C. 2012, Icarus, 219, 515
Porco, C. C., Baker, E., Barbara, J., et al. 2005, Science, 307, 1226
Poschel, T., & Herrmann, H. J. 1995, EPL (Europhysics Letters), 29, 123
Pravec, P., Harris, A. W., Scheirich, P., et al. 2005, Icarus, 173, 108
Rein, H., & Latter, H. N. 2013, Monthly Notices of the Royal Astronomical Society,
431, 145
Richardson, D. C., Blum, J., Weinhart, T., et al. 2012, in AAS/Division for Plan-
etary Sciences Meeting Abstracts, Vol. 44, AAS/Division for Planetary Sciences
Meeting Abstracts, 105.06
Richardson, D. C., Leinhardt, Z. M., Melosh, H. J., Bottke, Jr., W. F., & Asphaug,
E. 2002, in Asteroids III, ed. W. F. Bottke, Jr., A. Cellino, P. Paolicchi, & R. P.
Binzel (University of Arizona Press), 501{515
Richardson, D. C., Michel, P., Walsh, K. J., & Flynn, K. W. 2009, Planetary &
Space Sciences, 57, 183
Richardson, D. C., Quinn, T., Stadel, J., & Lake, G. 2000, Icarus, 143, 45
Richardson, D. C., Walsh, K. J., Murdoch, N., & Michel, P. 2011, Icarus, 212, 427
Richardson, J. E., Melosh, H. J., & Greenberg, R. 2004, Science, 306, 1526
Richardson, J. E., Melosh, H. J., Greenberg, R. J., & O'Brien, D. P. 2005, Icarus,
179, 325
Robbins, S. J., Stewart, G. R., Lewis, M. C., Colwell, J. E., & Sremcevic, M. 2010,
Icarus, 206, 431
Robinson, M. S., Thomas, P. C., Veverka, J., Murchie, S. L., & Wilcox, B. B. 2002,
Meteoritics and Planetary Science, 37, 1651
Rosato, A., Standburg, K., Prinz, F., & Swendsen, R. 1987, Phys. Rev. Lett., 1038,
L58
Ryan, E. V. 2000, Annual Review of Earth and Planetary Sciences, 28, 367
Salo, H. 1992, Nature, 359, 619
|. 1995, Icarus, 117, 287
Salo, H., & Karjalainen, R. 2003, Icarus, 164, 428
Salo, H., Schmidt, J., & Spahn, F. 2001, Icarus, 153, 295
Sanchez, D. P., & Scheeres, D. J. 2012, Icarus, 218, 876
Schmidt, J., Ohtsuki, K., Rappaport, N., Salo, H., & Spahn, F. 2009, in Saturn
216
from Cassini-Huygens, ed. M. K. Dougherty, L. W. Esposito, & S. M. Krimigis
(Springer), 413
Schmidt, R. M., & Housen, K. R. 1987, International Journal of Impact Engineering,
5, 543
Schroter, M., Ulrich, S., Kreft, J., Swift, J. B., & Swinney, H. L. 2006, Physical
Review E, 74, 011307
Schulson, E. M., & Fortt, A. L. 2012, Journal of Geophysical Research (Solid Earth),
117, B12204
Schwartz, S., Richardson, D., & Michel, P. 2012, Granul. Matter, 14, 363
Schwartz, S. R., Michel, P., & Richardson, D. C. 2013, Icarus, 226, 67
Schwartz, S. R., Michel, P., Richardson, D. C., & Yano, H. 2014, Planetary & Space
Sciences, 103, 174
Spitale, J., & Porco, C. C. 2007, in Bulletin of the American Astronomical Society,
Vol. 39, AAS/Division for Planetary Sciences Meeting Abstracts #39, 460
Stadel, J. G. 2001, PhD thesis, University of Washington
Sukhorukov, S., & Lsetl, S. 2013, Cold Regions Science and Technology, 94, 1
Sullivan, R. J., Thomas, P. C., Murchie, S. L., & Robinson, M. S. 2002, in Asteroids
III, ed. W. F. Bottke, Jr., A. Cellino, P. Paolicchi, & R. P. Binzel (University of
Arizona Press), 331{350
Taguchi, Y.-H. 1992, Physical Review Letters, 69, 1367
Takeda, T., & Ohtsuki, K. 2009, Icarus, 202, 514
Tanaka, H., Ohtsuki, K., & Daisaka, H. 2003, Icarus, 161, 144
Tancredi, G., Maciel, A., Heredia, L., Richeri, P., & Nesmachnow, S. 2012, Monthly
Notices of the Royal Astronomical Society, 420, 3368
Tiscareno, M. S., Perrine, R. P., Richardson, D. C., et al. 2010, The Astronomical
Journal, 139, 492
Tremaine, S. 2003, The Astronomical Journal, 125, 894
Tsuda, Y., Yoshikawa, M., Abe, M., Minamino, H., & Nakazawa, S. 2013, Acta
Astronautica, 91, 356
Uehara, J. S., Ambroso, M. A., Ojha, R. P., & Durian, D. J. 2003, Physical Review
Letters, 90, 194301
Vincent, J.-B., Besse, S., Marchi, S., et al. 2012, Planetary & Space Sciences, 66, 79
Wada, K., Senshu, H., & Matsui, T. 2006, Icarus, 180, 528
Walsh, K. J., Richardson, D. C., & Michel, P. 2008, Nature, 454, 188
Weidenschilling, S. J. 2011, Icarus, 214, 671
Wisdom, J., & Tremaine, S. 1988, The Astronomical Journal, 95, 925
Yanagisawa, M., & Hasegawa, S. 2000, Icarus, 146, 270
Yano, H., Kubota, T., Miyamoto, H., et al. 2006, Science, 312, 1350
Yasui, Y., Ohtsuki, K., & Daisaka, H. 2012, The Astronomical Journal, 143, 110
Yeomans, D. K., Konopliv, A. S., & Barriot, J.-P. 1997, Journal of Geophysical
Research, 102, 23775
Yoshikawa, M., Kawaguchi, J., Fujiwara, A., & Tsuchiyama, A. 2015, in Asteroids
IV, ed. P. Michel, F. E. DeMeo, & W. F. Bottke (University of Arizona Press),
397{418
Yu, Y., Richardson, D. C., Michel, P., Schwartz, S. R., & Ballouz, R.-L. 2014, Icarus,
217
242, 82
Zhang, Y., Richardson, D., Barnouin, O., et al. 2016, submitted
218