Viktor Sergeyevich Safronov 11/XI/1917 - 18/IX/1999

1. Planetesimal Formation: Numerical Modeling of Particle Growth, Settling, and Collective Gas-Grain Interactions S. J. Weidenschilling, Planetary Science Institute Cambridge, UK, Sept 2009

2. Viktor Sergeyevich Safronov11/XI/1917 - 18/IX/1999 “S”Si Monumentum Requiris, Circumspice”“”

3. Outline: Numerical Models • Structure of a particle layer in the midplane of a laminar nebula, where settling is in equilibrium with shear-generated turbulence • Settling with coagulation; conditions for growth • Accretion of larger bodies: effects of initial planetesimal size

5. Solid bodies are not supported by the pressure gradient • Any solid body must move relative to the gas • Particles move toward higher pressure, i.e., generally inward • Radial and transverse velocities are size-dependent; relative velocity between any two bodies can be easily calculated - if they are isolated

7. Collective Effects • Particles settle to form a dense layer in the nebular midplane • Gravitational instability and collisional coagulation both require (for different reasons) this layer to be much denser than the gas • In such a layer, gas is dragged by the particles, and no longer moves at the pressure-supported velocity, but closer to the Kepler velocity • Relative velocities depend not just on sizes, but on the entire population of particles

8. Shear between the layer and the surrounding gas causes turbulence • The thickness of the layer and its vertical structure are determined by the balance between settling and turbulent diffusion • Shear-induced turbulence may halt settling at densities too low for gravitational instability

9. Response of the Gas to Particle Loading • Nakagawa et al. (1986) derived coupled equations: particles move inward as gas moves outward • Angular momentum is conserved: equal and opposite radial mass fluxes of particles and gas • The solution assumes laminar flow, with local balance of momentum at any given level • However, turbulent viscosity causes significant vertical transport of momentum, affecting the radial and transverse velocities of the gas • The momentum flux must balance overall, but does not hold locally

10. Ekman Length A characteristic length scale for the thickness of a turbulent boundary layer of a disk rotating in a fluid is the Ekman length LE, defined as LE = (t / K)1/2 where t is the turbulent viscosity, and  = K is the rotation frequency. After Cuzzi et al. (1993) we take t ~ (V/Re*)2/K, where Re*~102is a critical Reynolds number. This implies LE = V/(Re*K) where V is the velocity difference between the midplane and large Z. Turbulence decays over a distance LE.

11. Richardson Number The Richardson number (Ri) is a measure of the stability of a stratified shear flow. If a fluid element is displaced vertically, work is done against gravity and buoyancy, while kinetic energy is extracted from the flow due to the mismatch of velocity due to the shear. Ri is dimensionless, defined as The flow becomes turbulent if Ri < 0.25.

12. Rossby and Stokes numbers • In weak turbulence (Ri ~ 1/4), the eddy frequency is imposed by the system’s rotation frequency,  • In strong turbulence (Ri << 1/4), eddies have higher frequency, Ro K, where Ro ~ 10-102 is the Rossby number • The Stokes number St = te / where teis the response time to the drag force • Particle random velocity = Vturb /(1+St)

13. Vertical Transport of Momentum by Turbulent Viscosity (Youdin & Chiang 2004) • Particle concentration is greater nearer the central plane, so rotation is faster • If the layer is turbulent, viscosity is significant • The vertical velocity gradient causes upward flow of angular momentum • Gas in the midplane flows inward, while that near the surface of the layer flows outward • Inward and outward mass fluxes are equal, but higher particle concentration in midplane yields net inflow of particles, if particles and gas are perfectly coupled

14. Numerical Modeling of Particle Layer • Divide layer into a series of levels, with assumed particle abundance at t = 0 • Compute radial and transverse velocities of gas due to particle-gas momentum exchange, using Nakagawa model; assume additive for more than one size • Average velocities over Ekman length • Assume turbulent velocity is proportional to the velocity difference between local gas and particle-free gas, and a function of Ri • Turbulence propagates between levels, decaying exponentially over Ekman length • Assume eddy timescale is a function of Ri such that varies from K to 2 Ro K as Ri approaches zero

15. Compute stress tensor and gas radial velocity due to vertical shear • Add gas radial velocities due to particle-gas momentum exchange • Solve for particle radial and transverse velocities, generalized to include radial motion of gas • Compute net radial mass fluxes of particles (inward) and gas (outward); should have equal magnitudes • If particle flux is too large (small), increase (decrease) gas velocity due to particle-gas momentum exchange and solve again until fluxes balance • Distribute particles vertically by settling and diffusion • Iterate until a steady state is reached

25. Relative velocities in midplane for a mixture of mm- and m-sized bodies, vs. mass fraction in m-sized bodies

27. Collective Effects and Impact Velocities • Increased Vrel between small particles due to shear-induced turbulence • Impact speeds of small particles onto m-sized bodies decreased due to smaller V (~ 10 m/s instead of ~ 50 m/s) • Continued growth to > m-size results in decoupling, increased V, higher impact speeds again for small onto large • Growth may be limited by erosion, unless > m-sized bodies can accrete each other?

28. Inward flux of m-sized bodies causes net outward flux of small grains carried with gas; fast outflow in dense sublayer exceeds radial drift in the thicker layer

29. Streaming instability and pileups? Effective drift velocity varies inversely with surface density only for sizes > decimeter or surface density a few times nominal

31. Results of Equilibrium Layer Models • Relative velocities are not simple functions of particle sizes, but depend on the ensemble size distribution and abundances • Inward and outward flows of gas and particles at different Z, even for a single size • Multiple sizes result in strong stratification, size sorting • A small fraction of mass in large bodies can cause net outward flow of small particles, with radial mixing

32. Models with Coagulation and Settling • Vertical layering up to 2 scale heights, finer resolution near central plane • Logarithmic size bins from d = 10-4 cm; fractal density to 1 cm • Relative velocities: thermal, differential settling, turbulence (alpha, also shear-generated in midplane) • Gravitational stirring (Safronov number) and cross-section for large bodies

33. Impact Outcomes • Impact strength S (erg/g): if energy density > S, target is disrupted; shattering velocity Vshat = (2S)1/2 • If projectile much smaller than target: its mass added to target; mass of escaping ejecta proportional to impact energy • Excavation parameter Cex such that mex = 0.5 Cex mpV2 • Transition from net gain to erosion at critical velocity Vc = (2/Cex)1/2

34. For perfect sticking, large bodies coagulate and settle to central plane in a few thousand orbital periods

35. V = 52 m/s, Vturb = 0 S = 105erg/g, Vshat = 4.5 m/s, Vc = 10 m/s

36. V = 52 m/s, Vturb = 0 S = 106 erg/g, Vshat = 14 m/s, Vc = 32 m/s

37. V = 52 m/s,  = 10-5, Vturb = 5.3 m/s S = 105 erg/g, Vshat = 4.5 m/s, Vc = 32 m/s

38. V = 52 m/s, Vturb = 0 S = 104 erg/g, Vshat = 1.4 m/s, No Erosion

39. Results of Coagulation/Settling Modeling • Settling/growth timescale for large bodies : ~ few x 103 orbital periods • If impacts of small particles result in net erosion above a critical velocity Vc < V, growth can be halted • If Vturb is small ( ~ < 10-5), growth is possible even for very low impact strength (104 erg/g), if erosion is limited (experimental data?) • Critical parameter: Velocity threshold for net gain/loss when a small particle hits a much larger one?

40. Is There an Observational Constraint on Sizes of Original Planetesimals? • Canonical “km-sized” planetesimals from Goldreich and Ward (1972) model for gravitational instability of a dust layer; no reason to prefer that size • Asteroid belt experienced accretion of large embryos, dynamical depletion, and 4+ Gy of collisional evolution, but may retain some trace of its primordial size distribution

41. Morbidelli et al. (2009): • Excitation and depletion of early asteroid belt requires accretion of large (~ 104 km) embryos within a few My • Present-day size distribution shows excess of ~ 100 km bodies relative to a power law of equilibrium slope; had to have formed early • Survival of Vesta’s crust implies early belt was deficient in bodies ~ 10 - 100 km • Can accretion produce these features from some initial characteristic size?

43. Gravitational Accretion Code • Multiple zones of semimajor axis; 15 zones from 2 to 3.5 AU • Collisions for bodies in overlapping orbits; impact rates and velocities • Gravitational stirring of eccentricities and inclinations • Accretion, cratering, disruption depend on impact energy • Logarithmic diameter bins; fragments below minimum size are lost

44. d0 = 1 km: “bump is at ~ 10 km

45. d0 = 100 km: embryos too small

46. d0 = 10 km: too many asteroids ~ 10 km?

47. d0 = 0.1 km

48. Still To Do: • Vary initial size from 0.1 km; add dispersion about the mean • Add planetesimals over some interval, instead of instantaneously • Test behavior if mass ground down into small (< 10 m) fragments is recycled into new planetesimals instead of lost Planetesimals may have started small!