Low Speed Crack Propagation via Kink Formation and Advance on the Silicon ( 110 ) Cleavage Plane

James R. Kermode, Anna Gleizer, Guy Kovel, Lars Pastewka, Gábor Csányi, Dov Sherman, and Alessandro De Vita Warwick Centre for Predictive Modelling, School of Engineering, University of Warwick, Coventry CV4 7AL, United Kingdom Department of Materials Science and Engineering, Technion-Israel Institute of Technology, Haifa 32000, Israel Institute for Applied Materials IAM, Karlsruhe Institute of Technology, 76131 Karlsruhe, Germany Engineering Laboratory, University of Cambridge, Trumpington Street, Cambridge CB2 1PZ, United Kingdom School of Mechanical Engineering, Faculty of Engineering, Tel-Aviv University, Tel-Aviv, Israel Department of Physics, King’s College London, Strand, London WC2R 2LS, United Kingdom CENMAT-UTS, Via A. Valerio 2, 34127 Trieste, Italy (Received 18 June 2015; published 23 September 2015)

It has long been known that fracture occurs in brittle materials when the chemical cost 2γ of creating two new surfaces is exceeded by the elastic energy G released in the same process [1].Extending this purely thermodynamic picture, Thomson et al. [2] proposed that discrete atomic sites give rise to a periodic "lattice trapping" barrier.Within these models, the energy barriers opposing forward and reverse crack motion can be obtained analytically as a function of G [3][4][5][6].As G is increased above the Griffith critical value 2γ, the forward lattice barrier decreases, falling to zero at a higher load G þ .Similarly, the barrier to reverse motion decreases as G is reduced below 2γ, falling to zero at a lower load G − , so that at sufficiently low temperatures, cracks will remain stationary for loads in the lattice trapping range Fast dynamic fracture will typically follow imposing energy release rates above G þ once all barriers are zero.Below this value, theoretical lattice models in twodimensional brittle crystals at zero temperature indicate that no propagation will occur [7].A forbidden range of velocities for steady-state cracks is thus predicted, known as the "velocity gap" [8,9].It arises as a dynamical consequence of lattice trapping: the snapping of each bond has to generate and pass along to the next bond sufficient kinetic energy to induce its breaking, preventing the crack from falling into the trapped state.A number of molecular dynamics simulations [10][11][12] and some single-crystal experiments [10,13] in silicon have suggested that a velocity gap may remain at low temperatures, with the crack speed jumping from zero to ∼2000 ms −1 during loading.
Calculations of the lattice trapping range suggest G þ ≃ 1.2 × 2γ, the exact value depending on the system size and force model assumed [6,14,15].Lattice trapping on this scale would imply that perfect brittle crystals are much stronger than the Griffith strength 2γ and that cracks, once they move, always move very fast [9].These two conjectures have not been universally confirmed by experiments, some of which have detected cracks moving at less than 1% of the Rayleigh speed at both room temperature and 77 K [16].We note, however, that these experiments were carried out in air, and the measured speed is low enough to enable propagation by stress corrosion [17].
Most brittle fracture modeling work to date has been restricted to narrow, quasi-2D model systems, in which all bonds located along a perfectly straight crack front break simultaneously.An alternative 3D propagation mechanism analogous to the kink-based motion of dislocations [18] was first proposed by Sinclair [4] and later developed by Marder [19].A more recent application to calculate the MEP for Si(111) crack advance was carried out by Zhu et al. [20], who found that the barrier for the formation of a kink pair is much lower than the lattice trapping barrier for 2D sequential advance of the entire front.This produced a possible explanation for the directional crack propagation anisotropy on the Si(111) plane, attributed to differences in the energy required to create kink pairs in different orientations [20,21].A 3D molecular-dynamics-based study of fracture in a complex metallic alloy also observed the formation of kink pairs [22].Possibly due to prohibitively high computational cost, none of these studies were based on first-principles methods.A similar mechanism of propagation by kink-pair formation has been observed in the wetting of periodic structures, where the motion of the contact line is analogous to that of dislocations or brittle crack fronts [23].
Here, we combined density functional theory (DFT) calculations and classical molecular dynamics (MD) simulations to investigate the issue of brittle crack propagation in the low speed, finite temperature limit.We used the NEB [24,25] approach to identify transition states and to calculate potential energy barriers for various kink propagation mechanisms.Classical MD trajectories were generated using the Stillinger-Weber (SW) potential [26], with the three-body term strength increased to ensure a brittle model material [27,28].
Specifically, we studied the unit advance of a crack on the Sið110Þ½1 10 (where round brackets denote the cleavage plane and square brackets the propagation direction) cleavage system for an N sites ¼ 10 unit-cell-deep model system with dimensions 300 × 100 × 54.3 Å 3 .The system is periodic along the crack front direction and contains 81 120 atoms, with a "fixed grips" load applied by displacing and fixing the upper and lower free surfaces [29], initially at the Griffith load 2γ so that the total potential energy of the system is exactly the same before and after each crack advance step.For this model, the energy activation barrier to the quasi-2D advance mode is 0.3 eV per unit cell [Fig.1(a)].
No 2D mode can apply to activated fracture propagation in real systems, as the arbitrary number of unit cells along the crack front would imply a diverging total energy barrier [30].It is, however, interesting to investigate the load dependence of the barrier E simul a (green dashed line in Fig. 2).This decreases as the load is increased, eventually falling to zero at G simul þ ¼ 5.45 J=m 2 .In a realistic 3D picture, the bonds are broken one after another by kink-pair formation and advance.We computed the full MEP for this sequential crack advance mode [Fig.1(b)] to evaluate the energy barrier E form a to kink formation, i.e., to break the first bond on an initially straight front [31], and the barrier E adv a to kink advance, which we model by breaking the fifth bond out of ten.Both values were obtained for a number of different applied loads above 2γ.As the load was increased, the barriers to kink advance and formation reduced, eventually falling to zero at respectively (Fig. 2).We note that the closeness of the computed values is consistent with the quasi-2D picture of crack advance where kink pairs form simultaneously at all points along the crack front.
At loads G > 2γ, a net total energy decrease follows kink advance, which can be evaluated as ΔE adv ðGÞ ¼ ðG − 2γÞΔA where ΔA ¼ 14.7 Å 2 is the area of new surfaces exposed by kink advance.The barrier for kink retreat, traversing the same MEP backwards and, thus, uphill (magenta curve in Fig. 2), is larger than the advance one (red curve) for all G > 2γ.Similarly, the barrier to kink healing to recover a straight front (cyan curve) can be determined by traversing the kink formation pathway backwards.At low loads, where G=2γ only slightly exceeds unity, these barriers are low, suggesting that transient kink healing competes with kink advance in the crack propagation kinetics in this load regime.
Although recalibration of these results using DFT calculations is needed before attempting direct comparison with experiments, as discussed further below, three propagation regimes can be identified from the kinetic barrier behavior represented in Fig. 2. Namely, (I) for low loads 2γ < G < G adv þ , isolated kinks are expected, since the rate of kink formation is much lower than the rate of kink kink growth is barrierless and essentially instantaneous, but formation is still relatively slow, leading to large numbers of interacting kinks, and (III) for loads G > G form þ , kink formation also becomes barrierless, so that many kinks can nucleate.The existence of regime (I) suggests a very slow "unzipping" mode is possible for crack advance.In this regime, kink formation is very slow, and we speculate that for loads just above 2γ, it may occur only at material defects.To check that these predictions remain consistent with a dynamical picture, we performed molecular dynamics simulations in a skewed model crack system where kinks can propagate without interacting [32].We observe the unzipping mode over a range of energy release rates below G form þ .The dependence of the kink migration rate on temperature fits an Arrhenius expression for temperatures from 50 to 300 K, with rates compatible with those predicted from static barrier calculations.To check the reliability of the SW potential and enable direct comparison with experiments, we next performed barrier calculations based on a quantum mechanical (QM) approach.A full QM calculation is not practical for the ∼10 5 atom model system sizes needed to adequately capture stress concentration in the three-dimensional strip geometry.Efficient calibration of classical results can, however, be achieved by carrying out DFT calculations for the two-dimensional propagation mode (cf.green dashed line in Fig. 2).For this, we used a hierarchical multiscale approach based on the methods of Spence et al. [14] and Perez and Gumbsch [15].Namely, we computed lattice trapping barriers directly with DFT [33] and the modified SW potential, using the NEB method for the two-dimensional propagation mode in small crack tip clusters with boundary conditions set by the Irwin nearfield displacements from the continuum anisotropic linear elastic fracture mechanics solution [36].The elastic moduli of the SW potential and of the DFT approach were used to compute the energy release rates G and the corresponding stress intensity factors K I for a range of applied loads (Table I).
Remarkably, we find that barriers computed with DFT and SW collapse into a universal curve when the G axis is rescaled to align G simul þ and the E simul a axis is rescaled to align the barrier heights at 2γ [37].The data presented in Table I indicate that the DFT lattice trapping barriers are slightly lower than those predicted by the modified SW model, 499 indicates that the range of lattice trapping predicted by DFT is much more limited.This is in agreement with previous results for the range of lattice trapping which did not include direct calculations of energy barriers [15].Although three-dimensional DFT NEB calculations for the kink formation and advance processes described above are currently not viable, the observed universal form of the E a ðGÞ function (see the Supplemental Material [37]) suggests that it is reasonable to assume that the other crack advance processes can be rescaled according to the DFT G simul þ and E simul a ð2γÞ values in the same way.This leads to the estimated DFT kink barriers shown in Fig. 3.
To investigate dynamical effects, we computed crack speeds using classical molecular dynamics with the modified SW potential over a wide range of energy release rates and system sizes.NVE simulations were carried out at 300 K in a 524 × 171 × 263 Å 3 740 000 atom model system.Cracks were propagated from rest at a G ¼ 6.1 J=m 2 load for 15 ps to allow the system to reach its equilibrium state.The applied strain was then rescaled to enforce a smaller load in the 4.25 < G < 5.45 J=m 2 range, and for each load, the simulation was continued for a further 20 ps.Kink-based propagation was observed in all these low-G simulations with no evidence of a velocity gap at this temperature.Crack velocities measured from MD simulations are shown in Fig. 3, where G=2γ has been rescaled to align it with the DFT G þ as described above.The error bars on the speeds measured in the classical MD were estimated from a sensitivity analysis assuming an  Arrhenius dependence of speed on barrier height and that the MD NEB barriers, even after calibration to match the DFT ones as described above, are accurate to no better than 0.02 eV.Using larger model systems with sizes up to 1231 × 246 × 785 Å 3 containing up to 12 × 10 6 atoms leads to a very similar energy-speed relation [38].
These calculations predict that thermally activated kink motion in silicon can give rise to crack propagation at speeds of hundreds of m=s under moderate applied loads, although propagation in the few hundred m=s speed regime has not been so far reported.To test this prediction, we, thus, carried out a novel set of experiments with silicon crystal on the same cleavage system.Namely, the energy-speed relationship was examined for the low energy regime on the ð110Þ½1 10 low energy cleavage system of silicon.The experiments were carried out in an inert argon environment to prevent chemically activated ("stress corrosion") cracking [17].We evaluated the crack speed from the fracture surface features using the Wallner lines method, calculating the energy release rates G at each point of the specimen by finite element analyses [39].Our results are shown as black points and error bars in Fig. 3.
To compare our calculations and experiments, we next used two simple analytical models for vðGÞ relevant to the lower and upper speed range limits considered here.For low speeds, Lawn [5] has shown that in the steady state where kink formation and annihilation via collisions are in equilibrium, thermally activated kink motion predicts the speed dependence: where a x is the lattice spacing in the propagation direction [1.35 Å for the Sið110Þ½1 10 fracture system] and are the kink advance and kink formation rates, respectively, where the Arrhenius factors following the ν 0 attempt rate were calculated from the MEP results reported in Fig. 3.For higher speeds corresponding to barrier-free propagation, we used the continuum elastodynamic crack motion equation proposed by Freund [40], which takes the form where c R is the Rayleigh shear wave speed computed from the model elastic constants (here, 5440 m=s), and 2γ is the fracture resistance, here assumed constant for all the low speeds considered.The velocities predicted by Eqs. ( 1) and ( 4) for T ¼ 300 K are plotted as purple and blue dashed lines, respectively, in Fig. 3. Choosing ν 0 ¼ 6.8 ps −1 aligns the two expressions at G form þ and is consistent with attempt rates measured in MD simulations.The fully thermally activated model of Eqs. ( 1)-( 3) and our MD simulation provide a lower bound for propagation in the low-load interval considered, while the continuum dynamics solution provides an upper bound for the same, with opposite curvature.
Our new experimental results fall between these two bounds, exemplifying low speed propagation in Si crystals and are consistent with a crossover from thermally activated to dynamic fracture.However, a precise quantitative comparison between our models and experiments remains challenging due to the large (≥ 0.2 J=m 2 ) uncertainty of experimental G values (horizontal error bars in Fig. 3) and our use of idealized pure-Si models.Low speed fracture is observed over a range of low speeds below G simul þ where our predicted kink formation and migration barriers are sufficiently small to allow thermally activated motion.The predicted speed of this motion (magenta curve in Fig. 3) is a lower bound since these kink barriers could still be too high as additional low-barrier mechanisms leading to higher crack propagation speeds cannot be ruled out.In particular, real specimens are likely to contain defects that may act as additional low-barrier kink-pair sources.Since, once formed, kinks can propagate any distance along the crack front, even a moderate density of defects would suffice to boost crack speed in the load regions where intrinsic kink formation is the limiting step in our simulations.If these sources were associated with very low or vanishing kink formation barriers, kink advance would become the rate limiting step between G adv þ and G form þ .We expect the picture of thermally activated crack motion to apply generally to other crystalline materials, although the atomic details of kink structure will vary from one material to another.This is at variance with the common intuition informed by the everyday experience of fast catastrophic failure of brittle materials driven by stress intensification yielding barrierless crack propagation.The possibility of a velocity gap for speeds below ∼2000 ms −1 occurring at finite temperatures, which could be inferred on the basis of the results of previous simulations, is not confirmed here, once kink propagation and low-barrier kink-pair generation are enabled by a fully 3D model system geometry.

FIG. 1 (
FIG. 1 (color online).Minimum energy paths (MEPs) obtained with the nudged elastic band (NEB) method at the Griffith critical load 2γ on the Si(110) cleavage plane for (a) crack advance by simultaneous advance of all N sites ¼ 10 bonds along the crack front and (b) crack advance by sequential bond rupture.Solid blue circles are minima and black crosses the positions of replicas.Schematic insets show broken (white) and intact (black) bonds along the crack front.

FIG. 2 (
FIG. 2 (color online).Dependence of kink-pair energy barriers calculated with the modified SW potential on G=2γ.Vertical dotted lines correspond to the loads G adv þ and G form þ at which the kink advance and formation barriers fall to zero.

FIG. 3 (
FIG. 3 (color online).Estimated DFT energy barriers for crack growth processes (left vertical axis) overlaid with experimental Si(110) and classical MD crack speed data (right vertical axis).The G axis has been rescaled to align with DFT G simul þ , and the dashed purple and blue lines show simple analytical models relevant in the low and high speed limits (see text).

TABLE I .
Comparison of elastic and lattice trapping properties of Si ð110Þ½1 10 cracks using the modified SW potential and DFT the generalized gradient approximation (GGA-PW91).