Transcript N-body Models of Aggregation and Disruption

N-body Models of Aggregation
and Disruption
Derek C. Richardson
University of Maryland
QuickTime™ and a
TIFF (LZW) decompressor
are needed to see this picture.
Overview





Introduction/the N-body problem.
Numerical method (pkdgrav).
Application: binary asteroids.
Non-idealized & strength models.
First results: “YORP” spinup of rubble
piles & spin limits with strength.
Introduction
 Many dynamical processes in the solar
system can be modeled by gravity and
collisions alone. E.g.,
 Reaccumulation after catastrophic
disruption (collisional or rotational).
 Planetary ring dynamics.
 Planet formation.
 Problems well suited to N-body code.
The N-body problem
The orbit of any one planet depends on the
combined motion of all the planets, not to
mention the actions of all these on each
other. To consider simultaneously all these
causes of motion and to define these motions
by exact laws allowing of convenient
calculation exceeds, unless I am mistaken,
the forces of the entire human intellect.
— Isaac Newton, 1687.
The N-body problem
Ýi   
rÝ
j i
Gm j (r i  r j )
ri  r j
3
Cost = N (N – 1) / 2 = O(N2)

Tree codes
 Reduce computational cost by treating
particles in groups.
Tree codes
Replace many summations with
single multipole expansion
around center of mass.
Tree codes
 Reduce computational cost by treating
particles in groups.
 Error controlled by opening angle
criterion and order of expansion.
Tree codes
Use multipole
expansion if opening
angle  < crit.

crit
Tree codes
 Reduce computational cost by treating
particles in groups.
 Error controlled by opening angle
criterion and order of expansion.
 Particles organized into systematic
hierarchical structure.
 Ideally suited for recursive algorithms.
Tree codes
E.g. Barnes & Hut (1986)
two-dimensional tree.
Cost = O(N log N)
Reducing cost further
 Parallel methods:
 Distribute work among Np processors.
 N-body problem difficult—exploit tree.
 Adaptive/hierarchical timestepping:
 Focus work on most active particles.
 Good object-oriented code structure.
 Hard-core optimizations.
Integrating the equations of
motion
 Many techniques for solving coupled
linear ordinary differential equations.
 Most popular:
 Runge-Kutta (explicit forward).
 Bulirsch-Stoer (complex/expensive).
 Leapfrog/symplectic methods.
 Preserve phase space volume.
 Timestep adaptability issues.
Collision detection

Particles collide when separation
distance equals sum of radii.
R1
R2
Collision detection


Particles collide when separation
distance equals sum of radii.
Two approaches:
1. Predict collisions before they occur.

Need neighbour-finding algorithm (tree!).
2. Detect collisions after they occur.


Detected by mutual overlap.
Adaptive timestepping essential.
Numerical method
 Our group uses pkdgrav:
 Parallel k-D tree code.
 k-D: split along longest dimension.
 Expand to hexadecapole order.
 Second-order leapfrog integrator.
 Hierarchical timestepping.
 Collisions predicted before they occur.
 Includes bouncing and sliding friction.
Parallelism in pkdgrav
master
• controls overall flow
“mdl”
“pst”
• loops over processors
“pkd”
• loops over particles on
one processor
• interface between
pkdgrav and parallel
primitives (e.g. mpi)
Application: binary asteroids
 Use N-body code to simulate:
 Capture of collisional ejecta in Main Belt.
 Michel et al., Durda et al.: collisions that make
families also make satellites.
Application: binary asteroids
QuickTime™ and a
mpeg4 decompressor
are needed to see this picture.
Michel et al. 2001
Application: binary asteroids
 Use N-body code to simulate:
 Capture of collisional ejecta in Main Belt.
 Michel et al., Durda et al.: collisions that make
families also make satellites.
 Rotational disruption of gravitational
aggregates in near-Earth population.
 Tidal disruption.
 “YORP” thermal spin-up.
Application: binary asteroids
QuickTime™ and a
YUV420 codec decompressor
are needed to see this picture.
Tidal disruption vs. YORP
 Tidal disruption makes binaries, but also
destroys them quickly.
 Binary NEA mean lifetime only ~ 1 Myr.
 YORP thermal effect may form binaries
through rotational disruption.
 But, some internal strength/cohesion may
be necessary to prevent material from just
“dribbling” away (but that may be OK too!).
Forming binaries with YORP
 Preliminary
investigation:
 Slowly spin up
various rubble
piles.
 Find particles
leak away from
equator (no
fission).
QuickTime™ and a
YUV420 codec decompressor
are needed to see this picture.
Forming binaries with YORP
 Preliminary
investigation:
 Slowly spin up
various rubble
piles.
 Find particles
leak away from
equator (no
fission).
Recoil: new mobility
mechanism?
QuickTime™ and a
YUV420 codec decompressor
are needed to see this picture.
Forming binaries with YORP
Forming binaries with YORP
 May need
strength and/or
irregular body
shape to form
binaries.
 E.g., contact
binary can
separate.
QuickTime™ and a
YUV420 codec decompressor
are needed to see this picture.
Non-idealized models
 Treating particles as idealized, rigid,
independent spheres is convenient.
 Components with different shapes may
provide more realism. E.g.,
 Ellisoidal particles (Roig et al.)
 Polyhedral (Korycansky & Asphaug).
 We combine best of both worlds: allow
spheres to “fuse” together…
Non-idealized models
QuickTime™ and a
YUV420 codec decompressor
are needed to see this picture.
QuickTime™ and a
YUV420 codec decompressor
are needed to see this picture.
Strength model
 Colliding particles/aggregates can:
 Stick on contact;
 Bounce;
 Liberate particle(s) from aggregate(s).
 Outcome currently parameterized by
impact speed.
Strength model
 In addition, bonded aggregates can
have a size-dependent bulk tensile
and/or shear strength.
 Particles experiencing stress (relative to
center of mass) in excess of strength
are liberated.
 Global model (no fractures/cracks).
Strength model
For a demo of the new strength model
in action, see Patrick’s presentation!
Testing strength: spin limits
 One way to test the strength model is to
compare with analytical predictions of
global failure (e.g. Holsapple).
 Found good match for cohesionless
models (Richardson et al. 2005).
 Science motivation: spin-up past critical
limit could make binaries (e.g. YORP).
Spin limits: preliminary results
Work in progress!
QuickTime™ and a
TIFF (LZW) decompressor
are needed to see this picture.
QuickTime™ and a
TIFF (LZW) decompressor
are needed to see this picture.
QuickTime™ and a
TIFF (LZW) decompressor
are needed to see this picture.
QuickTime™ and a
TIFF (LZW) decompressor
are needed to see this picture.
Summary
 N-body methods allow modeling of complex
phenomena involving gravity & collisions.
 Examples include post-disruption
gravitational reaccumulation to form binaries
& families.
 Binaries: more work needed to assess YORP
(including survivability against BYORP!).
 New pkdgrav strength model provides added
realism/complexity, but needs fracture model.
Extra Slides…
What is YORP?
 Yarkovsky-O'Keefe-Radzievskii-Paddack effect.
 Irregular bodies reflect/re-radiate solar photons in
different directions: net torque  spin-up/down.
Results: Many Binaries
 High rates of
production for:




Low q.
Low v∞.
Rapid spin.
Large
elongation.
Close
approach
distance q
Spin
period P
Encounter
speed v∞
Elongation ε
Orbital Properties
Eccentricity e
 High
eccentricity.
 Range of semimajor axis.
 Binary orbit
aligned more
with approach
orbit than
progenitor spin.
 Retrograde
orbits possible.
(97% > 0.1)
Inclination I
Retrograde
Semimajor
axis a
(50% > 10 Rp)
Spinorbit
angle
Physical Properties
 Size ratio peaks at
0.1–0.2 (10–5:1).
 Obliquities:
Size
ratio
Obliquities
 Primary spin aligned
with binary orbit.
 Wide range of
secondary spin axes.
 Spin Periods:
 Primary has narrow
range (3.5  6.0 h).
 Secondary has wide
range (4.0  20+ h).
Spin
periods
Evolutionary Effects




Mutual tides damp
eccentricity in ~ 1–10 My.
Repeated encounters
may strip binary.
NEA population refreshed
by MBAs (some of which
may be binary).
Thermal effects (YORP)
important?
Steady-state (Monte Carlo)
Model
 We know…
 Binary production efficiency from tidal disruption (Walsh
& Richardson 2006);
 Planetary encounter circumstances (Bottke et al. 1994);
 Distribution of NEA lifetimes (Gladman et al. 2000);
 Shape and spin of source bodies (Harris et al. 2005);
 Tidal evolution effects (Weidenschilling et al. 1989);
 Effects of binary encounters with planets (Bottke &
Melosh 1996; this work);
 Small binary MBAs formed in collisional simulations
(Durda et al. 2004).
Steady-state (Monte Carlo)
Model
 In one timestep…
 All asteroids in the simulation are tested for:
 End of lifetime (median ~ 10 Myr);
 Close planetary encounter < 3REarth (one every ~3 Myr).
 All binaries are tested for:
 End of lifetime;
 Close planetary encounter < 24REarth: explicit 3-body integration.
 If neither happen, the binary is tidally evolved.
 Removed NEAs/binaries are immediately replaced.
 “Fresh” asteroids take spin/shape characteristics of
MBAs, with a variable percentage being binaries.
 MBA binaries have characteristics determined from
the Durda et al. 2004 simulations.
Steady-state Results
For 2000 asteroids:
 Find ~2% binary
fraction.
 Binary NEA mean
lifetime ~ 1 Myr.
 93% of removed
binaries destroyed by
planetary encounters.
 MBA initial binary
percentage has little
effect (mean lifetime
~0.32 Myr).
Steady-state Results
 The resultant
steady-state
binaries…
Observed
 Have slightly larger
semi-major axes
than observed;
Steady-state
Steady-state Results
 The resultant
steady-state
binaries…
 Have slightly larger
semi-major axes
than observed;
 Mostly have low
eccentricities (< 0.2),
consistent with
observations.
Eccentricity
A Word About Rubble Piles
 Rubble piles are low-tensile-strength,
medium-porosity gravitational aggregates.
 In simulations, rubble piles consist of perfectly
smooth spheres; some dissipation.
 Used in a variety of contexts: planetesimal
collisions, tidal disruption, spin-up.
 How do they differ from perfect fluids?
Rubble Pile Equilibrium Shapes
Mass loss: 0% < 10% > 10%
X = initial condition
Rubble Pile Equilibrium
Shapes
Mass loss: 0% < 10% > 10%
X = initial condition
YORP Spinup of Rubble Piles
QuickTime™ and a
YUV420 codec decompressor
are needed to see this picture.
QuickTime™ and a
YUV420 codec decompressor
are needed to see this picture.
Resolution Effects
QuickTime™ and a
YUV420 codec decompressor
are needed to see this picture.
Classifications
Richardson et al. 2003
Stress response may be predicted by plotting tensile
strength (resistance to stretching) vs. porosity.
Strength vs. Gravity
Asphaug
et al. 2003
Aggregates Resist Disruption
 Once shattered, impact energy is more
readily absorbed at impact site.
Coherent
Damaged
Asphaug et al. 1998
QuickTime™ and a
TIFF (LZW) decompressor
are needed to see this picture.