Supersonic Dislocation Bursts in Silicon

Dislocations are the primary agents of permanent deformation in crystalline solids. Since the theoretical prediction of supersonic dislocations over half a century ago, there is a dearth of experimental evidence supporting their existence. Here we use non-equilibrium molecular dynamics simulations of shocked silicon to reveal transient supersonic partial dislocation motion at approximately 15 km/s, faster than any previous in-silico observation. Homogeneous dislocation nucleation occurs near the shock front and supersonic dislocation motion lasts just fractions of picoseconds before the dislocations catch the shock front and decelerate back to the elastic wave speed. Applying a modified analytical equation for dislocation evolution we successfully predict a dislocation density of 1.5 × 1012 cm−2 within the shocked volume, in agreement with the present simulations and realistic in regards to prior and on-going recovery experiments in silicon.


Dislocation Density Derivation
In 1958, Smith 1 proposed a shock front interface composed of supersonic dislocations. However, it did not predict a dislocation density increase due to shock which was resolved by sequential homogenous nucleation of dislocations at the shock front. This was calculated analytically by Meyers et al. [2][3][4] for Cu.
Here we adapt the analytical description developed by Meyers and leave in a strain-rate dependent HEL in addition to an increased separation of nucleation separation due to supersonic dislocation motion.
Dislocation density can be defined as the inverse separation of dislocations in two dimensions, lateral to the shock front, d, and aligned with the shock front, hs. Accounting for two dislocations per stacking fault, the following principal relationship for dislocation density can be written as: (1) Taking a standard expression for lateral separation of dislocations (d) in order to account for a given strain, ε, the following expressions can be written in terms of original and shocked lattice parameter (a0 and as respectively): The separation of dislocations can then be expressed in terms of lattice parameter, The inverse of lateral separation can be manipulated into the following form: Taking a typical relationship between Burgers vector, b, and lattice parameter, in addition to a relationship between instantaneous volume, V, and equilibrium volume, V0, as related to the shock and equilibrium lattice parameter, 3 3 00 We can now write an expression for d -1 in terms of only one variable, V: Following formalism previously developed 2 , the contribution of the stress field from each dislocation only produces stress in the σ12 component of the stress tensor. The sum of these contributions for a spacing of 1/n 2 dislocations extending infinitely along a planar shock front is equal to π 4 /90 and thus the expression for stress in terms of subsequent emission distances, h, can be written as: The critical shear stress, σ12 c , required to nucleate a dislocation can be defined through the HEL: Setting σ12 = σ12 c , the separation of subsequent partial dislocation nucleation sites is: The effective separation, hs, is then increased when accounting for mobile dislocations at the shock front which aid in relaxation: Here, k is an orientation factor equal to 1 for <001> shocks in fcc or diamond cubic materials. The shock speed, Us, is well defined as: While the minimum supersonic dislocation velocity relative to the shock front, vd, is defined as: The effective shock separation can now be written as:

Partial Dislocation Identification
In order to identify partial dislocations in the system, the primary identification scheme utilized a coordination based evaluation. We employ a cutoff radius of 3.0 Å which captures the addition of "new neighbors" during a partial displacement shuffle. This cutoff value must be below the expected neighbor distances in many polymorphs of silicon (3.2 Å) as well as the equilibrium distance of the second neighbor shell (~3.8 Å) which may significantly decrease during compression. Notably, this distance also corresponds to the bond distance between silicon atoms in many amorphous and 5 coordinated structures.
We show below a comparison of analysis techniques ( Fig S1) and a snapshot of partial identification near the shock front ( Fig S2).  Figure S2. Identification of partial dislocation near shock front by potential energy criterion. The partial dislocation tip is identified by the furthest dark blue atom corresponding to a shift in the lattice across the stacking fault.

Dislocation Density Calculation from Simulation
In order to extract an accurate dislocation density from the simulation, OVITO's surface identification tool 7 was used to relate surface area and volume of defective atoms. When defects have the thickness of only a single stacking fault they can be considered as thin disks of thickness, t. The relationship between surface area and volume, Rs/V, provides an expression for the circumference, and thus the partial separation, sp, and the dislocation line length, ld : (21) Relationship between uniaxial and hydrostatic stress. Figure S3. Relationship between "zz" component of the stress tensor and hydrostatic pressure compared during a steady state shock (created through a ramp) and Cij calculations for a system undergoing uniaxial compression.