## Abstract

We present a detailed account of an atomistic study of three-dimensional lattice trapping barriers to brittle fracture in Si. By means of a prototypical interatomic potential model, we map out the molecular details of the evolution of atomically sharp cracks in the (111) cleavage plane with straight crack fronts along the and directions, respectively. The thermally activated processes of bond rupturing along the crack front are quantitatively characterized using a reaction pathway sampling scheme. The calculated minimum energy paths reveal a mechanism of kink-pair formation and migration in facilitating the crack front advancement. We show that the physical origin of directional anisotropy in cleavage crack propagation can be attributed to a difference in the kink-pair formation energy for different crack orientations. The effects of interatomic potentials are delineated by comparing the Stillinger–Weber model with an environment-dependent model.

## 1. Introduction

Ample experimental evidence (e.g. Wiederhorn 1967; Lawn 1993) has indicated that cracks in many brittle solids are able to undergo quasi-static extension under sustained, constant loading, at least for a range of applied loads. One of the key microscopic steps controlling this time-dependent crack growth is the stress-mediated, thermally activated bond ruptures at the crack tip, either with or without the aid of reactive foreign species from the environment. The microscopic energy barriers that limit the kinetic rate of bond rupture are the so-called ‘lattice trapping’ barriers (Thomson *et al*. 1971). Specifically, the total energy of a cracked body under stress is not a smooth function of crack length on the microscopic scale. Owing to the periodicity of the discrete lattice, the energy landscape has an atomic-scale corrugation. As a result, a crack tip could be trapped in a local energy minimum on the energy surface, just as a dislocation is trapped by the Peierls barrier (Peierls 1940). The mobility of this crack is then controlled by the stress-mediated activation energies, measured by the barrier heights of the saddle points on the edge of the energy well. It is generally recognized that quasi-static crack extension is inherently a three-dimensional process involving a series of localized bond breaking events at the crack front (Lawn 1975; Thomson *et al*. 1987; Marder 1998). For a realistic estimate on the kinetic rate of crack growth, it is then essential to characterize the lattice trapping effect quantitatively in a three-dimensional setting.

Despite recent advances in experimental techniques for probing atomic structures near the crack tip (e.g. Celarie *et al*. 2003; Wiederhorn *et al*. 2003), computer simulation remains the most effective way of studying the atomic process of crack-tip bond breaking (e.g. Sinclair & Lawn 1972; Sinclair 1975; Spence *et al*. 1993; Pérez & Gumbsch 2000*a*,*b*; Bernstein & Hess 2003; Zhu *et al*. 2004*b*). As a continuing effort along this line, we report in this paper in full details the atomistic calculation of three-dimensional lattice trapping barriers to brittle fracture in Si. The study is carried out by systematically probing the potential energy surface of a cracked system calculated using the Stillinger–Weber (SW) interatomic potential (Stillinger & Weber 1985). We determine quantitatively the transition pathways of localized bond breaking at the stressed crack front using an effective reaction pathway sampling scheme, the nudged elastic band (NEB) method (Jónsson *et al*. 1998). The atomic geometries and energy variations along the transition pathways reveal the mechanistic role of kink-pair formation and migration along the crack front. We find that the origin of the directional dependence in cleavage crack propagation can be explained in terms of the differences in the energetics of kink-pair formation for two different crack orientations. Using the geometries and energetics of crack front kink pairs presented in this paper as input, it is conceivable that a coarse-grained study of the time evolution of crack front morphology becomes now feasible along the lines of Lawn (1975) and Cai *et al*. (2000).

## 2. Model and method

### (a) Geometry

We consider an atomically sharp crack on the (111) plane with a straight crack front along the and directions, denoted as and crack, respectively (see figure 1*a*). Figure 1*b*,*c* shows the relaxed atomic configurations near the crack tip for the two crack orientations at respective Griffith loads, determined later in §§3*a* and 4. For both orientations, the atomistic simulation cell is a cylinder cut from the crack tip, with a radius *R*=80 Å. Atoms within 5 Å of the outer surface are assigned a prescribed displacement field given by the anisotropic linear elastic Stroh solution (Stroh 1958; Suo 1990), and all remaining atoms are free to move. To bring out the three-dimensional nature of crack front propagation, the simulation cell length along the cylinder axis is taken to be suitably long, 20 unit cells, with periodic boundary condition (PBC) imposed. Because the 12-atom unit cell for Si is orthorhombic, the length *l* along the crack front and the number of atoms *N* in the cell are different for the two crack orientations. For the crack, *l*=76.8 Å and *N*=77 200; for the crack, *l*=133.0 Å and *N*=133 760. In both cases, the coordinate system is such that the *x*_{3}-axis is along the crack front, the *x*_{2}-axis is [111] and the cracks extend in the direction along the *x*_{1}-axis.

### (b) Interatomic potentials

The main results presented in this paper are calculated using the three-body SW potential (Stillinger & Weber 1985). The material properties, including lattice parameter, elastic constants and relaxed {111} shuffle plane1 surface energy, are listed in table 1. For comparison, these properties are also calculated using the environment-dependent interatomic potential (EDIP) (Bazant *et al*. 1997; Justo *et al*. 1998). While the SW potential was fitted mainly to bulk diamond cubic and liquid Si in an empirical fashion, the EDIP potential was more rationally constructed with sp^{2}- and sp^{3}-hybridized bonding as targets at the very beginning, and fitted to a much larger database. Both potentials, however, describe the more local σ-bonding comparatively better in practice than the more delocalized π-bonding. Owing to the short-range nature of the two potentials (Holland & Marder 1999; Bernstein & Hess 2001), the cleaved {111} surfaces are the ideal 1×1 type without reconstruction. Figure 2 shows the cohesive responses of uniformly cleaving a perfect crystal of Si along the {111} shuffle plane into two blocks. Here, surface traction is plotted as a function of opening displacement normalized by the equilibrium interplanar spacing . We note that in order to integrate these curves to obtain the correct surface energy, the tractions calculated from the short-range SW and EDIP models are generally larger than the real value needed for the same interplanar separation (Holland & Marder 1999; Bernstein & Hess 2001). Comparing the two potentials, besides a spurious shoulder on the rising branch of the EDIP potential, significant differences are seen in the peak stress as well as in the slope on the falling side of the cohesive response. Beyond the interplanar separation of respective traction peaks, the traction calculated from the SW potential gradually decreases to zero, while the response from the EDIP potential exhibits a sharp drop. The influence of the interatomic force law on the lattice trapping effect has been investigated by Sinclair (1975), Thomson *et al*. (1987), Curtin (1990) and Gumbsch & Cannon (2000). Using simplified interatomic potentials, they found the slope of the falling side of the force law critically affects the loading range within which the lattice trapping effect exists (Thomson *et al*. 1987). The effects of the SW and EDIP potentials on lattice trapping barriers will be discussed in §5.

### (c) Finding transition pathways of crack extension at different loads

Each crack is subjected to a pure mode I load as given by the stress intensity factor *K*_{I}. Finding the transition pathways and associated activation energies for crack-tip bond rupture requires probing the potential energy surface of the system in a multi-dimensional configuration space under *K*_{I} loading. We apply the NEB method (Jónsson *et al*. 1998; Henkelman & Jónsson 2000) to sample the transition pathways. The NEB calculation requires a knowledge of both the initial and final states; one needs to identify two different local energy minima on the same potential energy surface mediated by *K*_{I}. For the present simulation setup, both the initial and final states are well defined. They correspond respectively to the equilibrium atomic configurations before and after crack extension in the *x*_{1} direction by one atomic spacing, denoted as Δ*a*. After the initial and final states are identified using the energy minimization scheme, a discrete elastic band consisting of a finite number of replicas (or images) of the system is constructed by linear interpolation to connect the two end-states (Jónsson *et al*. 1998). Note that the *K*-field load is applied via the displacement control method, which requires the positions of atoms at the outer boundary of the simulation cell to be the same for both the initial and final states, as well as for those intermediate replicas along the pathway. Then a spring interaction between adjacent replicas is added to ensure continuity of the path, thus mimicking an elastic band. An optimization of the band, involving the minimization of the forces acting on the replicas, brings the band to the so-called minimum energy path (MEP). MEP is defined as a continuous path in a 3*N*_{f}-dimensional configuration space (*N*_{f} is the number of free atoms), with the property that at any point along the path the atomic forces are zero in the 3*N*_{f}−1-dimensional hyperplane perpendicular to the path (Sorensen *et al*. 2000). The energy maximum along the MEP is the saddle-point energy, which gives the activation energy barrier. The calculation is considered to be converged when the potential force on each replica vertical to the path is less than 0.002 eV Å^{−1}.

Of particular note is that in NEB calculation, the reaction coordinate is unknown *a priori* and emerges as a natural result from the converged pathway. Specifically, we define the hyperspace arc length along the MEP, i.e. , as the reaction coordinate, where denotes a state in the hyperspace. Empirically, the bond length of the Si–Si bond undergoing rupture has been taken as an approximate reaction coordinate (e.g. Bernstein & Hess 2003). This approach of calculating the energy barrier essentially corresponds to finding the maximum point along a path which is constructed by connecting all the energy minima within a series of hyperplanes vertical to a specific direction, i.e. the direction that corresponds to the degree of freedom represented by the Si–Si bond distance in a multi-dimensional hyperspace. Clearly, it is only in some cases (e.g. is always positive) that this path agrees with the MEP defined above and thus gives an exact result for the activation energy barrier. Our calculations confirm that due to the localized nature of Si bond breaking, the path obtained from the calculation using the empirical reaction coordinate of the Si–Si bond length is close to the actual MEP in that they give the same transition states and barrier heights. Thus, the corresponding more abstract, though physically rigorous, reaction coordinate obtained from the NEB calculation can be approximately understood as the bond length of the breaking Si–Si bond. In any event, the NEB method represents a robust way of calculating the barrier height; it can be directly applied to study more complicated transition pathways when the process of bond rupture involves an interaction with foreign chemical species (e.g. Zhu *et al*. 2005).

## 3. Formation and migration of a kink pair at the crack front

In this section, we report on the study of a (111) crack with a straight crack front along the direction using the SW potential. This crack orientation has been used in most two-dimensional simulations of Si cleavage fracture in the literature.

### (a) The Griffith load *K*_{IG}

We first determine the Griffith load *K*_{IG}, at which the change in the system total energy Δ*E* is zero when the crack extends in the direction by one atomic spacing Δ*a*, with the in-plane borders of the simulation cell fixed. For the present setup with 20 unit cells along the direction (PBC), 20 bonds need to be broken along the crack front. The value of *K*_{IG} can be estimated from the relaxed surface energy *γ*_{s} using the Griffith relation in linear elastic fracture mechanics (LEFM) (e.g. Rice 1978),(3.1)where _{c} denotes the critical strain energy release rate and Δ*A*=*l*Δ*a* is the area of newly created surface. For a linear elastic, anisotropic solid, _{c}, is related to *K*_{IG} by(3.2)where *H*_{22} is the effective compliance of the cracked system on the *x*_{1}−*x*_{2} plane; it is the 22 component of the matrix * H* defined by eqn (2.13) in the paper of Suo (1990). Note that

*is defined in terms of elastic constants in the global coordinate system (*

**H***x*

_{1},

*x*

_{2},

*x*

_{3}), where the crack front is along the

*x*

_{3}-axis. For the crack, the computed

*H*

_{22}is 2.692×10

^{−11}Pa

^{−1}using the elastic constants given by the SW potential as listed in table 1. Then substitution of

*H*

_{22}and surface energy

*γ*

_{s}into equations (3.1) and (3.2) gives an estimate of the Griffith load . On the other hand, direct atomistic calculation detailed in §3

*b*gives a numerical result of . The good agreement can be taken as a measure of the accuracy in our atomistic simulations.

### (b) Two-dimensional versus three-dimensional crack extension at *K*_{IG}

Two competing pathways are studied that advance the crack along the direction by one atomic spacing. We focus on comparing the overall features of the two pathways in this subsection, further discussion of the three-dimensional pathway will be given in §3*c*.

The first pathway is characterized by *simultaneously* breaking every bond along the crack front. This is essentially a two-dimensional fracture mode, since what happens within a unit cell is repeated along the crack front direction. The corresponding MEP can be simply calculated using a two-dimensional crack configuration. That is, only one unit cell is needed along the crack front with PBC imposed. Obviously, the activation energy will grow linearly with the thickness of the simulation cell in the crack front direction. For later comparisons, figure 3*a* shows the energy variation along the MEP of breaking 20 bonds simultaneously. Here, symbols represent energies of the replica configurations, while the continuous curve is constructed by cubic-polynomial interpolation. The reaction coordinate is the normalized hyperspace arc length along the MEP. It is seen that in advancing the crack by one atomic spacing along the direction, the net change in total energy between the initial and final states is zero at . This Griffith load is obtained by trial-and-error search. Specifically, we determine *K*_{IG} by varying the value of *K*_{I} until an approximate equality between initial and final energy states is obtained; the value of *K*_{I} that accomplishes this being equal to *K*_{IG}. Zero energy change at *K*_{IG} is due to the cancellation between the elastic energy decrease and the surface energy increase. Note that the MEP shown in figure 3*a* has only a single peak; the activation energy of about 12.85 eV corresponds to an average of 0.64 eV per crack front bond.

In contrast to the first pathway which is two-dimensional in nature, the second pathway involves breaking 20 bonds *sequentially* along the crack front, which is inherently a three-dimensional process. The specific pathway of sequential bond ruptures that we study involves breaking a single bond at the centre of the 20-bond cell, and then breaking, one by one, bonds to either side of the centre bond; other competing three-dimensional pathways will be discussed in §5*b*. Figure 3*b* shows energy variation along the MEP of sequential bond breaking, where the normalized reaction coordinate *s* is defined such that each integer number *s* labels a locally equilibrated state with *s* broken bonds on the crack front. Between *s* and *s*+1, the coordinate denotes a normalized hyperspace arc length along the MEP of further breaking the (*s*+1)th bond. Correspondingly, each circle in figure 3*b* represents the energy, relative to that of the equilibrium state with a straight crack front, of a metastable state of local equilibrium with *s* broken bonds on the crack front. The curve with a single peak connecting two adjacent circles gives the energy variation along the MEP of further breaking one more bond. For clarity, we plot only the interpolated curve to indicate continuous energy variation along the MEP, and leave out the discrete data points corresponding to the calculated energies of 15 relaxed replica configurations between two adjacent circles. Note that the overall kinetic rate for the second mechanism can be determined by the maximum barrier height which corresponds to the saddle point on the MEP of breaking the 11th bond. The overall kinetic barrier is therefore about 2.02 eV for the second pathway.

Comparing energy variations along the two transition pathways at *K*_{IG}, which have the same initial and final states, we quantitatively confirm that even for the present simulation cell with nanoscale thickness (76.8 Å) along the crack front, the three-dimensional extension mode is far more favourable kinetically than the two-dimensional mode for cleavage crack propagation. In the rest of the paper, we will therefore focus on the energetics of three-dimensional extension under different load levels and for different crack orientations.

### (c) Geometry and energetics of forming a crack front kink pair at *K*_{IG}

We now consider the three-dimensional cleavage mode in detail. Figure 4 shows the profile of the crack front for a representative state of local energy minimum with 10 broken bonds (*s*=10). A continuous field of crack opening displacement across the two adjacent (111) cleaving planes is rendered by cubic-spline interpolation of the openings at discrete lattice sites. Readers should be aware that this is just a useful way to visualize completely discrete data from atomistic calculations. It can be seen that two kinks of opposite signs with sharp features are developed locally on each side of the advancing crack front. This highly localized mode of crack opening distribution is in distinct contrast with our previous result of crack-tip nucleation of an embryonic dislocation loop in an FCC Cu crystal, which exhibits a significant spread of shear fault distribution across the glide plane (see fig. 3*b* in Zhu *et al*. 2004*a*). The difference in the range of opening versus shearing distribution at the crack front can be correlated to that in the bonding characteristics: Si has directional, localized covalent bonding, while Cu has delocalized metallic bonding.

Given the crack front with *s* broken bonds shows the characteristics of a double kink, we next describe the energetics of crack front evolution in terms of kink-pair formation and migration energies. Consider first the kink-pair energies as given by the circles in figure 3*b*. Now we extract these data points and replot them in figure 5*a*; the envelope of figure 3*b* is shown in figure 5*a*. Note that the new abscissa is taken as the kink-pair separation *d*_{K}≡*s*×*l*_{0}, where *l*_{0} denotes the unit-cell length in the crack front direction, being equal to , and the integer number *s* is the reaction coordinate defined in §3*b*, which labels a local equilibrium state with *s* broken bonds on the front. We also show in figure 5*a* the kink-pair energies for a larger calculation consisting of 30 unit cells along the crack front (*l*=115.2 Å), but with otherwise identical setup. It is seen that in both cases, to break the first bond at the straight crack front requires an input of extra energy from the thermal reservoir. This trend of energy increase continues with subsequent bond breaking. However, the energy increment associated with further bond breaking will gradually decrease as the separation of kink pair *d*_{K} increases. Consequently, the kink-pair energy approaches an asymptotic value as *d*_{K} becomes large enough. To model the effect of kink-pair separation, the system total energy (in reference to a straight crack front) can be partitioned into *d*_{K}-dependent and *d*_{K}-independent parts, the former is the elastic interaction energy *E*_{el} between two kinks of opposite signs separated by *d*_{K}, while the latter is twice the self-energy *E*_{K} of an individual kink. We write(3.3)for the ideal situation of a kink pair embedded on an infinitely long, and otherwise perfectly straight, crack front. In equation (3.3), *E*_{el}(*d*_{K}) will vanish asymptotically as *d*_{K} increases. *E*_{K} then contains all the remaining atomistic energetic information, and can be interpreted as the *formation energy* of an isolated kink on the crack front.

From figure 5*a*, the kink formation energy *E*_{K} estimated from the plateau of the solid line, which comes from the simulation cell of *l*=76.8 Å, is about 0.80 eV. This value is very close to an estimate of 0.81 eV from the thicker simulation cell (*l*=115.2 Å). Furthermore, figure 5*a* shows that a further increase of *d*_{K} beyond *l*/2 will result in a decrease in the kink-pair energy. This is due to image interactions, as we are not actually in the ideal situation of equation (3.3) of an infinitely long crack front. Specifically, when *d*_{K}>*l*/2, the attractive interaction from the image of the other kink in the neighboring simulation cell will begin to dominate, leading to a decrease of the total energy. Because, by definition, the initial and final states have the same energy at *K*_{IG} in advancing the crack by one atomic spacing in the *x*_{1} direction, i.e. the crack front at *x*_{1}=0×Δ*a* and 1×Δ*a* before and after advancement, it is straightforward to show that after the image interactions have been taken into account, the kink-pair energy in the PBC cell should be roughly symmetric with respect to *d*_{K}=*l*/2, as numerically shown in figure 5*a*. Note that our calculations also show that *E*_{el}(*d*_{K}) is a rapidly decreasing function with *d*_{K}, so even in the presence of image interactions we can extract *E*_{K} to be around 0.8 eV at *K*_{IG}.

Consider next the kinetic barriers for kink migration. Recall that each curved segment connecting two adjacent circles in figure 3*b* represents the MEP of a kink moving laterally by one atomic spacing. To facilitate comparison of kink migration barriers at different values of *d*_{K}, we show in figure 5*b* the MEP of breaking the 1st, 2nd, 5th, 8th and 10th bond, respectively. These curves in figure 5*b* are essentially the folded and magnified profiles of the ‘cusps’ in figure 3*b*. The energy reference for each curve is that of the initial state of breaking an individual bond, so that all curves start from zero energy. Also shown in figure 5*b* are the data points of the energies of replica configurations along each pathway. It is seen that breaking each individual bond from the 1st to the 10th is thermodynamically unfavourable in that the final state has a higher energy than the initial state. Moreover, the kinetic barrier for a forward transition will decrease, from 0.86 eV of breaking the first bond to 0.44 eV of breaking the 10th bond. As the number of broken bonds, i.e. the separation between two kinks *d*_{K}, increases, the migration barrier will approach an asymptotic value corresponding to the activation energy for the migration of an isolated kink, denoted by *W*_{K}. This value is given approximately by the barrier for breaking the 10th bond, 0.44 eV. In addition, we compare the present *d*_{K}-dependent kinetic barrier of breaking a bond with the barrier height per bond for the two-dimensional pathway in §3*b*, in which all bonds at the crack front break simultaneously. It is of interest to observe that the kinetic barrier of 0.64 eV per bond from the two-dimensional calculation is close to the average of the upper and lower limits from the three-dimensional calculation.

We conclude this subsection by noting a useful analogy between the crack front kinks and dislocation kinks, the latter having been extensively studied (e.g. Lothe & Hirth 1959; Caillard & Martin 2003). For crystals with significant secondary Peierls barrier, such as Si, the fact that the kink mechanism should also play a central role in crack front mobility is understandable given the structural similarity between the crack front, as the core of a crack tip, and the core of a dislocation line. Indeed, it is widely recognized that in Si, dislocations move by the nucleation and migration of kink pairs. Exploiting the analogy between the crack front and the dislocation core in an explicit and operational fashion, as suggested by our results, should be considered in future studies of deformation physics of both types of defects.

### (d) Loading effect and lattice trapping range

In the literature, the Griffith load *K*_{IG} is also regarded as the neutral load (e.g. Sinclair 1975), because the initial and final states have the same energy in advancing the crack by one atomic spacing. As the applied load increases beyond *K*_{IG}, the potential energy landscape of the system will be tilted, such that a forward transition becomes favourable both thermodynamically and kinetically. This loading effect can be clearly seen from the MEP of sequentially breaking 20 bonds at a higher load of , as shown in figure 6. For comparison, the same MEP shown in figure 3*b* at the Griffith load is replotted in figure 6, but adjusted to the new energy scale. It is seen from the MEP at that there exists a critical distance between a pair of kinks, denoted by , beyond which further kink-pair separation will cause the system energy to decrease, which is thermodynamically favourable. Clearly, is finite if and only if *K*_{I}>*K*_{IG}. At , the kink pair with a separation of corresponds to the state with one broken bond on the front as shown in figure 6. In addition, it is also seen that the energy barrier for kink migration (moving in the direction of separating the kink pair) decreases compared to the corresponding barrier height with the same *d*_{K}, but at a lower load of *K*_{IG}.

To characterize *K*_{I}-dependent energetics of kink-pair formation and migration, two more critical loads, beyond the Griffith load *K*_{IG}, can be considered. The first one, denoted by , is related to the thermodynamic energy balance of kink formation. At , the energy change associated with the formation of the smallest possible kink pair, i.e. breaking the first bond beyond the straight front, vanishes. For the present simulation setup, . The MEP of breaking the first bond at is shown in figure 7*a*. Evidently, the energy change is zero between the initial and final states. But there is still an activation barrier of about 0.43 eV at this load level.

The second critical load, denoted by , is the *athermal* load which is related to the kinetics of kink-pair formation. That is, at , the activation barrier for breaking the first bond at the straight crack front vanishes. More importantly, our calculations of MEPs at different load levels have indicated that, compared to the barrier heights of later breaking other bonds, the kinetic barrier of breaking the first bond is the slowest decreasing one as *K*_{I} increases. Hence, above , the first kink-pair nucleation as well as subsequent kink migration across the crack front can take place spontaneously without the aid of thermal fluctuations. This is the scenario that leads to spontaneous cleavage fracture.

To determine , we increase the load incrementally above *K*_{IG}. At each load level, the corresponding MEP of breaking the first bond is obtained via the NEB calculation. Figure 7*a* shows energy variations along MEPs at three typical loads. Here, the symbols denote the calculated energies of replica configurations along the pathway and the continuous curves are constructed via interpolations. We find that when the load reaches , it is no longer possible to obtain a relaxed initial state geometrically similar to the configurations at the lower loads; the system tends to relax to another local energy minimum corresponding to a different deformation mode, amorphization by forming new crack-tip ring structures. Since the focus of this work is to study the transition pathway for cleavage fracture along the (111) plane, we leave the issue of competition among different deformation modes to a future investigation. Nevertheless, we can estimate by extrapolation from the activation energies Δ*E*_{act} obtained at the lower loads, see figure 7*b* where the circles are the calculated data points and the solid line is the polynomial fitting curve. It is seen that Δ*E*_{act} for breaking the first bond shows a slight nonlinear dependence on *K*_{I}. The value of , obtained by extrapolating the curve to zero activation energy, is about . This athermal load is the upper limit of the lattice trapping effect. Hence, the lattice trapping range given by the SW potential, defined here as the ratio , is about 1.55. In particular, we note that since the three-dimensional pathway is more favourable kinetically than the corresponding two-dimensional pathway at the same load level, the present estimate of lattice trapping range should be a lower bound to the two-dimensional estimate. While the result of lattice trapping range is dependent on the interatomic interaction model, the calculations and analyses presented here serve to demonstrate how three-dimensional load-mediated lattice trapping barriers to brittle fracture can be accurately characterized atomistically.

## 4. Propagation anisotropy

Although the crack (propagating in the direction) is the most frequently studied crack orientation, fractography observations of cleavage surface indicate that, instead, a {111} cleavage crack prefers to propagate in the 〈110〉 direction (e.g. George & Michot 1993). The preference of propagation direction *within* a certain cleavage plane, called *propagation anisotropy*, has been explained by comparing orientation-dependent lattice trapping ranges of two-dimensional crack configurations (Pérez & Gumbsch 2000*a*,*b*). Here, we provide a more detailed explanation of propagation anisotropy in terms of the orientation-dependent energetics of crack front kink-pair formation and migration.

We first note that it is fundamentally difficult to rationalize propagation anisotropy using only thermodynamic arguments. The Griffith criterion may suggest that preferential cleavage plane should be the one with the lowest surface energy *γ*_{s}. But, on the same cleavage plane, two cracks of different orientations would experience the same resistance *γ*_{s} to crack propagation and thus have the same critical energy release rate _{c} as can be clearly seen from equation (3.1). Indeed, there may exist a difference in the critical stress intensity factor *K*_{IG}, which arises from a variation of the effective compliance within respective *x*_{1}−*x*_{2} planes due to elastic anisotropy for different crack orientations. However, as we now show, this difference is so small that *K*_{IG} cannot be regarded as a robust measure of propagation anisotropy.

On the other hand, consider a practical situation where the material contains an assortment of flaws distributed in size and aspect ratio. It seems unlikely that a special pre-crack orientation is selected *every time* to first trigger unstable propagation, just due to the condition of the initial flaws. A more plausible scenario is that with increasing macroscopic load, an ensemble of crack orientations exceeds the Griffith limit at approximately the same time and they all start to grow, aided by thermal activations. Since the load *K*_{I} (for simplicity, we only consider mode I load) is still less than the corresponding , there will be a special orientation whose front grows the fastest due to smaller kink-pair formation and migration energies. If the system is under load control, *K*_{I} will increase with crack length in a positive feedback, until this particular crack front reaches the athermal threshold of first. At that point, runaway dynamical fracture occurs. Also, this orientational competition may occur not only at physically separate crack fronts, but also at different places in one contiguous crack front, leading to morphological evolution that favours faster growth of one orientation until dynamical fracture is triggered first in that orientation. Therefore, we see that although the lattice trapping ranges could be quite narrow in brittle materials, it is still necessary to do a careful study, because a physically important question of which orientation triggers dynamical fracture first could depend on the precise value of .

To continue this line of reasoning, we next consider the crack extending in the close-packed direction. The crack-tip atomic structure is shown in figure 1*c*. Similar to §3*a*, we estimate the Griffith load using the analytic LEFM Stroh solution. The effective compliance calculated for the new crack orientation is *H*_{22}=2.71×10^{−11} Pa^{−1}. Substitution of *H*_{22} and the (111) surface energy *γ*_{s} into equations (3.1) and (3.2) gives an estimate of the Griffith load . On the other hand, the value of *K*_{IG} determined from direct atomistic calculations is . Evidently, the value of *K*_{IG} for the crack is very close to that of for the crack, which proves numerically that the Griffith criterion cannot be applied to explain propagation anisotropy.

Although the difference in *K*_{IG} for different crack orientations is small, the difference in may be large. We begin by studying the energetics of crack front kink-pair formation and migration for the crack at its Griffith load . Figure 8*a* shows the MEP of sequentially breaking 20 bonds along the crack front. For comparison, the corresponding MEP for the crack at its Griffith load is also shown in figure 8*a*. Here, symbols represent the kink-pair formation energies, which are the sums of kink-pair interaction energy and twice of the kink self-energy, as given by equation (3.3), while the saddle points between symbols give the kink-pair migration barriers. Comparing the two MEPs, we find significantly different kink-pair formation energies, while the kink migration barriers are similar. Specifically, for the crack, the kink self-energy *E*_{K}, estimated from the asymptotic value of the lower envelope curve connecting the squares, has a much lower value of about 0.22 eV. In contrast, for the crack, the kink self-energy is much higher, at about 0.8 eV.

As the applied loads increase by the same ratio beyond *K*_{IG}, the activation energy required to extend the crack is consistently lower than that for the crack, as demonstrated selectively in figure 8*b* using two MEPs at . To estimate the lattice trapping range for the crack, we also obtain the athermal load by extrapolation. As shown in figure 7*b*, a value of is obtained, which then leads to , a narrower lattice trapping range for the crack. We note that the two-dimensional simulations by Pérez & Gumbsch (2000*a*,*b*) on the {110} cleavage cracks also demonstrate orientation-dependent lattice trapping ranges. However, the present three-dimensional study follows a different line of investigation.

The orientation-dependent kink-pair formation energy is closely related to the bond densities along different directions within the (111) plane. The profile of the kinked crack front for the crack at is shown in figure 9, where the distribution of opening displacement across the cleavage plane is plotted for the state with 10 broken bonds on the crack front. It is seen that both the leading and trailing edges of the crack front have a zigzag profile, in contrast to the atomically smooth crack front for the crack as shown in figure 4. The zigzag profile indicates that the second array of bonds immediately adjacent to the crack front also have relatively large opening displacements. This feature can be explained as follows. Assuming the LEFM solution is approximately correct near the crack front, the bond opening displacement (in the *x*_{2} direction) should vary sensitively with the distance (in the *x*_{1} direction) from the crack front. The close-packed arrangement of atoms along the *x*_{1} direction for the crack means a smaller distance of the second array of bonds from the crack front, and, consequently, larger equilibrium opening displacements and a more zigzag profile. Moreover, we believe this is consistent with the crack having a smaller kink formation energy in that one can expect the energy cost of barrier crossing by the crack front kink to become smaller if tensile opening displacements before and after the crossing are closer.

## 5. Discussions

### (a) Effect of interatomic potential

We have selectively repeated the calculations we have just presented using the EDIP potential to study the influence of the interatomic potential on results. The Griffith loads are first determined for two crack orientations. Atomistic calculations give for the crack and for the crack. Similar to the SW potential, a slight variation of *K*_{IG} between the two crack orientations arises due to the difference in the effective compliance within respective *x*_{1}−*x*_{2} planes. However, the SW and EDIP potentials predict different Griffith loads, which is consistent with the fact that they give different values of elastic constants and surface energy as listed in table 1.

Figure 10 shows the energy variations along the MEPs at respective Griffith loads for the two crack orientations calculated using the EDIP potential. Comparing to the SW results shown in figure 8*a*, we find that the crack orientation is also favoured kinetically. Despite the agreement in the trend, there are some quantitative differences in the kink-pair formation and migration energies. Specifically, relative to the SW potential, the EDIP potential gives a lower kink formation energy, but a higher migration energy. Furthermore, our calculations of MEPs at different load levels using the EDIP potential indicate that the lattice trapping range is larger for both crack orientations compared to the corresponding results obtained from the SW potential. This difference in lattice trapping range between the two potentials can be correlated with the different cohesive responses of uniformly cleaving a perfect crystal as shown in figure 2, keeping in mind that the actual crack-tip bond breaking occurs under a non-uniform stress environment with a large stress gradient. Recall that figure 2 shows a significant difference in the slope on the backside of the cohesive response for the two potentials. Previous analyses from a series of simplified interatomic force laws indicated that a steep slope on the backside of the force law will lead to a large lattice trapping range (Thomson *et al*. 1987). Our present calculations show the same trend quantitatively. Finally, we note that a recent two-dimensional study by Bernstein & Hess (2003) on Si shows that the SW potential overestimates lattice trapping barriers compared to the results obtained from a multiscale simulation, i.e. the tight-binding description of bonding near the crack tip embedded in an empirical potential (EDIP) region. It is conceivable that a calculation using a more accurate force field along the lines presented here will give a more realistic estimate on the three-dimensional lattice trapping barriers.

### (b) Other competing transition pathways

In the present work, we focus on studying the lattice trapping barriers along a specific three-dimensional transition pathway, breaking the crack front bonds sequentially. Clearly, for a given *K*-field load, there exist other competing pathways to extend the crack. Among a few possible scenarios that we have considered, the pathway of sequential bond breaking is kinetically the most favourable. We have demonstrated in §3*b* that simultaneous fracture of crack front bonds requires a much higher activation energy. In this subsection, we study another competing process. Given *s* broken bonds on the front, we compare the energetics of breaking a bond on the advancing front versus breaking a bond on the trailing front. The latter corresponds to the scheme of sequential bond breaking we have studied before. Without loss of generality, we choose a starting point as follows: the crack under with *s*=8 broken bonds on the crack front; the EDIP potential is used in the calculation. Figure 11 shows the MEPs of the two competing processes of breaking the ninth bond on the crack front along with the MEPs of breaking the previous eight bonds. It is seen that both the saddle point and final state of the new scheme have higher energies than the corresponding values associated with the scheme of sequential bond breaking. Hence, bond breaking at the advancing front is thermodynamically unfavourable, and there is also a larger kinetic barrier. Moreover, the study of morphological stability of an initially straight crack front in a linear elastic solid by Rice (1985) indicates that for a crack under a displacement-control boundary condition, small perturbations from straightness should decay with time. Assuming this long wavelength analysis is applicable to atomic-scale kinks, an initially straight crack front should maintain straightness during quasi-static extension, which agrees with our atomistic calculations.

## 6. Summary

We have presented previously a detailed account of an atomistic study to quantify the energetics and atomic geometries of bond ruptures along the *three-dimensional* crack fronts in Si. We find the directional anisotropy in cleavage crack propagation can be attributed to a difference in the kink-pair formation energy for different crack orientations. Results from both the SW and EDIP potentials show a consistent trend of orientation-dependent lattice trapping barriers. We note that along the lines presented here, the kinetic process of dislocation emission at the crack front, which involves bond rupture and subsequent bond reformation by thermal activation (e.g. Zhu *et al*. 2004*a*), can also be explored. A quantitative characterization of the two competing processes of cleavage bond rupture and dislocation emission at a three-dimensional crack front is essential to the multiscale modelling of brittle-to-ductile transition in Si; the calculated atomistic energetics can serve as atomistic input for the coarse-grained simulation to study the crack front dynamics at a longer time-scale (e.g. Cai *et al*. 2000).

## Acknowledgments

We thank Ali S. Argon for stimulating discussions. T.Z. and S.Y. acknowledge support by NSF-ITR grant DMR-325553, AFOSR and Lawrence Livermore National Laboratory. J.L. acknowledges support by Honda R & D Co., Ltd, NSF grant DMR-0502711, AFOSR grant FA9550-05-1-0026 and ONR grant N00014-05-1-0504.

## Footnotes

↵† Present address: Woodruff School of Mechanical Engineering, Georgia Institute of Technology, Atlanta, GA 30332, USA.

↵For Si with a diamond-cubic structure, there are two types of {111} plane: one is the shuffle plane, which cuts through single covalent bonds along the direction perpendicular to the {111} plane, and the other is the glide plane, which cuts through triplets of covalent bonds inclined equally to the {111} plane.

- Received July 31, 2005.
- Accepted August 16, 2005.

- © 2006 The Royal Society