Traction curves for the decohesion of covalent crystals

Interest in the cleavage of covalent crystals has recently emerged in the semiconductor industry after demonstrations that controlled spalling techniques can be used to separate as-grown device materials from single crystal substrates.^1–4 Such spalling techniques allow for the reuse of expensive substrate materials and thereby promise significant cost savings in the manufacture of integrated circuits and photovoltaic devices. The approach exploits the fracture mechanics of tensile strained films on brittle substrates,^5–8 where cracks can be induced that run parallel to and just below the film/substrate interface.² The spallation crack follows a path where the mode II shear stress is zero and thus propagates under a purely mode I tensile stress.⁵ While brittle fracture is a complex multi-scale phenomenon, it ultimately depends on the physics and chemistry of bond breaking. In the case of controlled spallation of thin films, this bond breaking at the atomic scale occurs through the clean tensile cleavage of two half crystals.² 

In landmark papers, Rose et al.^9,10 studied the binding between two flat semi-infinite solids for a large number of metallic systems. They found that the curves describing the energy, u, as a function of separation distance, δ, or traction curves, exhibit a universal behavior. All the curves could be made to collapse by normalizing the energy, and expressing the separation in terms of an effective length λ. Furthermore, they constructed a simple expression for such Universal Binding Energy Relationship (UBER) with the form

where u₀ is the binding energy at equilibrium (δ = 0) and γ is the surface energy at infinite separation.

This relationship implies that the elastic constant in the direction of separation, C, and the surface energy, γ, are related. If d₀ is the equilibrium separation between the atomic planes, then

By evaluating the ultimate theoretical stress, the previous equation can be rearranged to eliminate the microscopic material parameter λ. Therefore, if UBER holds the elastic constant in the stretch direction C, the surface energy γ and the ultimate theoretical stress $σUTS$ are not independent but related by

Equation (1) works remarkably well in describing decohesion in metallic and bimetallic systems.^11,12 In fact, it can even be applied to a broader range of problems including chemisorption, the binding of diatomic molecules and even binding at the atomic nucleus.¹³ Banerjea and Smith¹³ elaborated an explanation of the origin of the universality of this relationship. In its most basic form, the main argument they present is that in all these problems, binding is due to the overlapping of the electron densities, which at some distance starts following an exponential decay. This basic argument explains how UBER works well in systems with non localized electron wave functions.

Here, we are interested in whether these ideas can be applied or extended to covalent crystals. Like in the original work by Rose et al.,^9,10 we will consider the rigid separation of two halves of the crystal. Local atomic relaxations can be included, as in Refs. 14 and 15, but the main objective of this work is to determine if any type of universality like the one found in metals exists. There are two points that we address. First, if there is any form of universality in the traction curves for different covalent crystals or cleavage directions, and second, whether there is a simple mathematical expression that can describe such traction curves.

To study these questions, in this work we perform periodic boundary condition, ab initio density functional theory (DFT) calculations within the generalized gradient approximation with the Perdew, Burke, and Ernzerhof¹⁶ (PBE) exchange correlation functional as implemented in the VASP code.^17–20 Calculations are carried out for different cell sizes ranging from 12 to 48 atoms, and the number of lattice planes perpendicular to the fracture direction are 24 for the (111) direction, 16 for the (100) direction, and 8 for the (110) direction. Energy cutoffs are set at 400 eV, while the k-space grid is set with a single point in the decohesion direction and 8 in the perpendicular directions, using a Gaussian smearing with a width of $0.05 eV$⁠. This k-space grid yields total bulk energies with an error of less than $1meV/atom$⁠. While maintaining the atoms in their relative bulk configuration a displacement δ is introduced in the middle plane and the energy-displacement curves are evaluated.

Figure 1(a) shows the traction curves for the decohesion of SiC. Different curves correspond to different cleavage planes: (100), (110), and (111). In the case of the (111) plane, in the zincblende structure of SiC there is a sequence of two types of inter-atomic planes that differ in the number of bonds per atom linking the atomic planes. There is a weak (w) interplanar space where atoms are linked by a single perpendicular bond and there is a strong (s) interplanar space where three oblique bonds link the atoms. For all the cleavage planes, the calculations yield binding curves that exhibit the same type of qualitative behavior as the crystal is pulled apart. However, unlike metallic crystals, a simple scaling of the energy and length scales does not result in a collapse of the calculated points onto a single curve.

In Fig. 1(b), we apply the scaling procedure of Rose et al.^9,10 First, the energies are scaled by the surface energy (⁠$2γ$⁠) for each crystallographic direction. Then, the separation distances (δ) are scaled so the curvature for small displacements coincide; for each curve, we define a scaling length λ as

Fig. 1(b) shows that after rescaling none of the curves for the different crystallographic directions follow the shape dictated by UBER (Eq. (1)) nor do they collapse to a single curve. As argued by Banerjea and Smith,¹³ the marked localization of the electron density in covalent crystals is likely the main reason why UBER fails to describe covalent crystals with the same accuracy as metals. Lazar and Podloucky²¹ reported results for diamond that show that UBER is not suitable for this system and argued that the directionality of the electronic bond is the reason behind the discrepancy.

From a practical point of view, it remains desirable to have a mathematical expression that can describe these traction curves. In principle, any type of interpolation will work, but it would be more useful to follow an argument along the lines of Rose et al.^9,10 and Banerjea and Smith¹³ to construct an extended form of UBER. To this end, we start by considering the expression of the force per unit area as derived from the UBER

We can interpret this equation as the multiplication of two terms: a linear force describing an elastic response for small displacements times an exponential tail. Following Banerjea and Smith, we can argue that the exponential tail is an essential term that has its origins in the shape of the electron wave functions of density functional theory. But we can then phenomenologically improve on the linear force term by adding higher order nonlinear terms. Therefore, we propose that the force is written as

This expression can be integrated to obtain the extended form of the binding energy relationship, which then becomes

Notice that the extended binding relationship, which we will refer to as xUBER, is written in a manner that is proportional to an energy scale and that it has a single length scale once the α_i coefficients are set, taking a form that can be useful to address universality, as we will see below. This traction curve therefore does not attempt to introduce several length scales for a single decohesion curve,²² and has recently been applied to model fracture in non-metal systems.^15,23 We can also point out that in this modified version of UBER the relationship between the elastic constant and the internal parameters changes from the form given by Eq. (2) to the expression

Figure 2 shows how this expression can be used to fit the data from Fig. 1. Here, the expansion has been taken up to the fifth power (that is, with parameters α₂ to α₅).

Even though there is no universality in the traction curves for different crystallographic directions, we can still ask ourselves whether there is universality when considering different covalent crystals with similar structure undergoing cleavage on a particular plane. Figure 3 shows the traction curves for (100) cleavage for a family of covalent crystals with the diamond and zincblende structures. For these materials, ranging from C to Ge, including the compounds SiC and GaAs, we observe a wide range of binding energies, but when the curves are normalized following the procedure of Rose et al.^9,10 as described above, we do observe a collapse of the data. Hence, we find that chemistry determines γ and λ as in the UBER, while crystallography and the geometry of decohesion determines the additional coefficients α_n with n > 1 of the xUBER. This is a weaker form of universality with the traction curves for covalent crystals depending on the geometry of the atomic structure but not on the exact chemical nature of its components. It should be noted that all these compounds exhibit qualitatively similar bonding and we have not established whether this universality will extend to chemistries that are not stable in the diamond (or zincblende) crystal structure. Furthermore, this universality is unlikely to hold if the nature of bonding across the cohesive zone changes as a function of separation distance, as it would happen during rehybridization. For example, in the (111) decohesion of carbon in the diamond crystal structure, our DFT calculations predict a nonmonotonic traction curve indicating a change in the nature of the bonding as the separation between the planes increases.

The calculations presented in the current work are based on the rigid separation of two halves of the crystal. Atomic relaxations can be taken into account following the procedures delineated in Ref. 15. A large fraction of atomic relaxations during uniform decohesion is in the form of elastic stretching within the crystalline material adjacent to the decohering region that has a negligible effect on bond breaking.^14,15 In the case of SiC, it was observed that upon subtraction of the elastic response of the adjacent crystals, the resulting traction curve can also be represented with the xUBER relation and shows no significative qualitative change in shape compared to the traction curve for rigid separation.¹⁵ This suggests that the universality relationship should also hold with the inclusion of internal atomic relaxations.

This work was funded by DOE under NEUP of the Idaho National Laboratory, Award No. 103195, and Project No. 10-924.