## Abstract

Discussed in this paper are several theoretical and computational approaches that have been used to better understand the fracture of both single-crystal and polycrystalline diamond at the atomic level. The studies, which include first principles calculations, analytic models and molecular simulations, have been chosen to illustrate the different ways in which this problem has been approached, the conclusions and their reliability that have been reached by these methods, and how these theory and modelling methods can be effectively used together.

## 1. Introduction and background

As a gem, diamond has been revered for close to five centuries for its beauty, relative rarity and its symbolism for wealth and affection [1]. Diamond gems also have an historical dark side that includes legends of being cursed, cartel-like control of the diamond market, and a history of paying for wars and other violence [2]. Diamonds possess a similar fascination from a materials science viewpoint [3]. Diamond has the highest hardness and thermal conductivity at room temperature of any natural bulk material. These extreme properties come, in part, from the strong covalent bonds in diamond, and from diamond having the highest atomic density of any existing crystal (under ambient conditions). These properties have made diamond attractive as a hard coating, for polishing applications and for thermal management. The colour of a diamond can vary from colourless to include pink, green, yellow, blue, red and purple depending on impurities, defect types and their densities. Diamond also has a relatively high refractive index and optical dispersion, which when combined with an optimized morphology can result in the ‘fire’ that is characteristic of otherwise colourless diamonds. From an electronic structure perspective, diamond is commonly considered to be a wide band gap semiconductor (the band gap is 5.54 eV). This property, combined with a high degree of radiation hardness, the high thermal conductivity, a negative electron affinity and transparency from the ultraviolet to the terahertz range, has motivated considerable research into using diamond for applications involving power electronics, electron field emission, optics and heat management.

Artificial diamond has been created since the 1950s by gas-phase chemical vapour deposition (CVD) [4] and by high-temperature, high-pressure (HTHP) methods [5]. While the CVD method predates HTHP synthesis by a few years, early CVD was a much slower method to produce diamond. This is because CVD methods rely on controlling the kinetics of diamond formation compared with the more thermodynamically stable graphite, while HTHP methods control environmental conditions to thermodynamically favour diamond. Breakthroughs in CVD methods in Japan in the 1970s followed by research starting in the 1980s in the USA and elsewhere led to advances in the deposition of diamond coatings and particles with CVD (figure 1, [6]) [4]. The key to most CVD methods is creating ubiquitous hydrogen atoms and activating the precursor carbon species using methods that include microwave plasma, hot filament, arc-jet and combustion. The hydrogen atoms saturate carbon atoms and create radical sites on the surface for deposition, and preferentially etch non-sp^{3} carbon [4].

Because diamond is a brittle material, its fracture properties are extremely sensitive to initial cracks, flaws and related microstructural elements [7,8,9]. In addition, the cost and inherent small sizes of single-crystal diamond make taking statistically relevant measurements to determine Weibull distributions challenging. Similarly, the wide range of microstructures observed in CVD-grown diamond coupled to a dependence of mechanical properties on film thickness creates a different set of experimental challenges in determining inherent fracture properties of diamond. For example, measurements on single-crystal CVD-grown diamond samples yielded a mean strength of 2860 MPa, with a standard deviation and maximum value of 1200 MPa and 5100 MPa, respectively [4,7]. The fracture strength measured for CVD-grown polycrystalline samples has been found to be substantially smaller than single crystals, with a relatively large dependence on film thickness and growth orientation [10]. It was concluded from these studies that both inter- and trans-granular crack propagation contributes to the fracture of polycrystalline films [4].

The comparable cohesive energy of sp^{2} and sp^{3} bonding in carbon creates a wide range of stable bulk structures with mixed bonding hybridization that have their own unique properties, some of which are comparable to crystalline diamond. These materials are typically characterized by their ratio of sp^{2} to sp^{3} bonding, their hydrogen content, their density, and the thermal, mechanical and electronic properties that arise from these bonding characteristics [11,12,13]. Ultra-nano-crystalline diamond (UNCD) has also received a lot of attention [14,15]. This material is characterized by grain sizes of between about 2 and 10 nm, and typically disordered grain boundary structures compared with those usually found in nano-crystalline metals. This material has been used as hard coatings, for bio-sensors, field-emission cathodes, electro-mechanical resonators and switches, and in biomedical devices [15,16]. Experimental studies have shown that UNCD has mechanical properties that are similar to diamond. UNCD, for example, has a measured Young's modulus of 967 GPa (compared with an average Young's modulus of 1200 GPa for diamond) and fracture strength of 4172 GPa [8]. The latter can be influenced by doping used to enhance the electrical conductivity. Like diamond, the fracture of UNCD also follows Weibull statistics as expected for a brittle material [8,14].

Rather than further discussing the experimental measurements of the fracture properties of diamond and related materials, this article focuses on the complementary theory and simulation of the intra- and inter-grain diamond fracture from an atomic viewpoint. As discussed below, this remains an area ripe for further studies due to the unique properties of carbon bonding that probably lead to orbital rehybridization and associated changes in chemical bonding during fracture, and the wide range of structures and defects with comparable energies possible in bulk carbon materials compared with other materials.

## 2. Intra-grain cleavage

Diamond exhibits brittle cleavage with a preference for {111} facets [3,17]. This cleavage preference results in natural diamond crystals typically having a roughly octahedral morphology, which is often used by craftsmen to create a rough cut prior to cutting a gem to a shape that optimizes its value. The traditional understanding of this cleavage preference is based on a bond counting argument and the diamond cubic lattice [3]. Cleaving along {111} planes of diamond requires breaking roughly 1.8×10^{19} carbon–carbon single bonds per square metre, while cutting along {100} and {110} planes requires cutting an additional 4×10^{18} and 1.3×10^{19} carbon–carbon bonds per square metre, respectively. Hence the argument is that fewer bonds needing to be broken results in the {111} cleavage preference.

This model for intra-grain cleavage is simple, intuitive and explains experimental observations regarding diamond fracture. Further analysis, however, suggests that there are open questions remaining regarding ideal brittle diamond fracture. Fracture leaves surface carbon atoms with unsatisfied sp^{3} bonds. To satisfy these bonds, and hence lower the surface energy, the orbitals of the affected carbon atoms can rehybridize to form surface pi bonds. For the (111) diamond surface, both theory and experiment suggest that this occurs through a reconstruction which results in the formation of chains of pi-bonded surface carbon atoms (figure 2) [18,19]. The reconstruction requires a topology change involving subsurface bond breaking and reforming. This change in structure also introduces considerable strain energy that partially offsets the gain in energy from the pi bonding. By contrast, on the (110) surface the surface atoms are already bonded along a chain configuration, and therefore forming partial surface pi bonds involves only a simple relaxation with no changes in bonding topology, and a relatively small amount of surface strain (because of the strain the pi bonds are not planar and therefore not full pi bonds) [20]. The (100) surface is intermediate between these two cases. Surface atoms can form surface pi-bonded dimers between what would be second neighbour atoms in the bulk, but at the cost of considerable surface strain and associated energy [18].

Given in the first column of table 1 are energies for the (111), (110) and (100) diamond surfaces derived from the bond-breaking model without any surface relaxation assuming a single bond energy of 3.79 eV (in this case, the cleavage energy is twice the surface energy). Also given in table 1 are surface energies recently calculated from density functional theory (DFT) using the B3LYP density functional [21]. These surface energies include ideal bulk terminated, relaxed and reconstructed surfaces. In contrast to the bond-counting arguments, the DFT calculations suggest that for ideal cleavage formation of (110) surfaces are energetically preferred, followed by (111) and (100) surfaces. As mentioned above, this is because the surface carbon atoms on the (110) surface whose bonds are broken during ideal cleavage have geometries that enable partial pi bonding. Relaxation of the surfaces is predicted by these calculations to produce (111) and (110) surfaces that are energetically similar. Reconstruction of the (111) surface to pi-bonded surface chains substantially lowers the energy of this surface so that it has the lowest energy of the three surfaces, while formation of a (2×1) pi-bond structure in the (100) surface lowers its energy so that this surface has an energy between the (111) and (110) surfaces. Hence the energy arguments are potentially more complicated than simple breaking of inter-planar bonds.

Rather than considering the energy, the maximum in the stress during cleavage can be used to predict cleavage preferences. This stress can be analytically calculated in a way that uses the surface energies via the universal binding energy relation (UBER). The UBER was introduced by Rose and co-workers [22,23,24] to describe a universal energy–distance relation that was observed for a wide range of deformation processes that include stretching of diatomic molecules, symmetric straining of crystalline solids and rigid separation of surfaces. When appropriately scaled by the minimum energy, and the distance and curvature at that energy, it was observed that all of these energy–distance relations fall on a curve described by the function
2.1
In this expression, *s* is a scaled distance given by
2.2
where *a* is the inter-planar distance and *a*^{0}_{ijk} is the minimum energy inter-planar separation. In an original formulation, the distance scaling factor *l*_{ijk} is given by
2.3
where *E*^{0}_{ijk} and (d^{2}*E*/d*a*^{2})_{ijk} are the minimum energy and curvature associated with the inter-planar separation at *a*=*a*_{o}. The maximum force per unit area during ideal cleavage can be taken as the derivative of equation (2.5) with respect to separation
2.4
which has a maximum value of *e*^{−1} at *s*=1. Applying the chain rule to equation (2.4) and converting to strain yields the analytic expression
2.5
where is the maximum stress for cleavage separating along the (*ijk*) plane, *E*^{s}_{ijk} is the energy of surface (*ijk*), and *Y* _{ijk} is Young's modulus in the direction normal to the (*ijk*) plane. This maximum cleavage stress derived from the UBER differs by a multiplicative constant from the long-established Orowan fracture criteria [25]
2.6

Given in table 1 are Young's moduli for the [111], [110] and [100] directions in diamond, the corresponding inter-planar spacings appropriate for fracture, and the estimates for the maximum ideal cleavage stress from equation (2.5) using each of the surface energies in the table. For ideal cleavage assuming no relaxation and using the DFT energies, the maximum stress for cleavage along the (111) plane is about 4% lower than that for cleavage along the (110) plane despite the energy for the (111) plane being higher than that for the (110) plane. For the relaxed surfaces, where the energy of the (111) and (110) surfaces are calculated to be comparable, the difference in the maximum cleavage stress between the (111) and (110) planes increases to over 8%. Hence this analysis suggests that it is the maximum stress and not the energy that determines cleavage preference among these three surfaces for defect-free diamond.

Using detailed semi-empirical and first principles DFT calculations, Paci *et al.* [26] simulated diamond fracture under an applied strain. Owing to the high computational costs of these simulations, a relatively small unit cell of 104 atoms was used that consisted of eight layers of atoms with uniaxial strain applied along the [100] direction. These simulations showed a smooth and continuous energy versus strain relation up to a critical strain, after which the diamond rapidly fractured. The simulations using the semi-empirical total energy method produced a relatively clean fracture with well-defined surfaces. For the simulations using DFT, on the other hand, the strain energy released from fracture resulted in the diamond lattice transforming into an amorphous structure. These simulations produced Young's moduli between 1.56 and 0.955 TPa, which are comparable to the experimental value of 1.08 TPa for diamond. The stress at fracture from the simulations was reported to be in the range 219–277 GPa depending on the details of the total energy method used in the simulations. These values are larger than the UBER results of about 175 GPa for fracture leading to the unreconstructed (100) surface. This discrepancy may be due to several effects, including the high strain rate and the relatively small cell sizes in the simulations, both of which may influence the frequency of the dynamics that initiate the fracture, or from deficiencies in the UBER scaling assumptions. In any case all of these values are substantially higher than those reported from single-crystal diamond, which is likely to be due to existing cracks and flaws in the diamond that are well known to substantially lower fracture strengths in experimental systems.

There remain several unanswered questions with respect to intra-grain cleavage of ideal diamond that could, in principle, be addressed by theory and modelling. First, it is not known where along the cleavage path bond rehybridization and relaxation of surface structures occur, nor is the influence these processes may have on the stress at cleavage [27]. Similarly, whether release of strain energy during cleavage can provide sufficient thermal energy to facilitate reconstruction, and if so what role this may have on cleavage preference, is also unknown. Finally, it is well established that chemisorption of hydrogen can inhibit (and even reverse) the surface reconstructions discussed above by terminating unsatisfied bonds without the need for additional surface stresses associated with surface relaxation and rebonding [18]. Whether chemical dynamics of this type plays a role in cleavage dynamics of diamond is unknown. Additional scaling relations have been proposed that reportedly provide a better description of ideal cleavage, for example by accounting for relaxation during deformation [28,29]. Whether such expressions can also account for chemical effects has not been established.

## 3. Inter-grain fracture

Like the intra-grain studies discussed in the previous section, there have been different approaches taken for atomic-scale modelling and theory of inter-grain fracture. Using the same semi-empirical and first principles methods used to simulate fracture of an ideal diamond lattice, Paci *et al.* [26], for example, simulated the fracture of an (100) interface containing a *Σ* 13 twist grain boundary. These simulations used a periodic unit cell containing 208 carbon atoms (16 planes of 13 atoms) with two twist boundaries. Upon straining in the [100] direction, fracture occurred abruptly at either one or both of the grain boundaries depending on the computational method used. These simulations yielded substantially lower stresses of 100–163 GPa at fracture compared with the companion single-crystal diamond studies mentioned above. It was also found that the stress and strain at fracture were somewhat dependent on the straining conditions used in the simulations [30]. While it is clear that existing flaws in polycrystalline diamond dominate the measured fracture stresses, these simulations demonstrate that weakening at grain boundaries can play a large role in the overall fracture behaviour of diamond polycrystals and UNCD.

Shenderova, Brenner and co-workers [31,32] carried out a similar but more extensive set of simulations using a reactive empirical bond-order (REBO) potential energy function. This expression gives the potential energy and inter-atomic forces as a function of relative atomic positions and allows for bond breaking and formation with appropriate changes in atomic hybridization, including pi bonding. The calculations require a fraction of the computational cost of electronic structure total energy methods and hence much larger systems can be explored, but at the potential cost of accuracy and transferability [33]. The simulations reported below were carried out using non-commercial molecular simulation code that has been thoroughly tested. Simulation step sizes and other details can be found in the original references.

Using the REBO potential, molecular dynamics simulations were carried out to explore the role of grain boundary structure, energy and the presence of existing cracks in the cleavage dynamics of diamond [27,34,35,36]. To study grain boundary structure without pre-existing cracks, these simulations used bi-crystals with grain boundaries of different types and orientations. Periodic boundaries were applied in two directions, while a uniaxial tensile strain was applied in the third direction by slowly straining layers of outer atoms whose relative positions were held constant. Similar to the simulations using the semi-empirical and first principles forces discussed above, fracture was accompanied by an abrupt drop in stress and energy at the critical applied strain. However, the large-scale disordering that was observed in the DFT calculations after fracture was not observed in these simulations. This may be due to the larger systems studied with the REBO potential (approx. 4000 atoms versus approx. 200 atoms) that are better able to dissipate kinetic energy generated from the release of the strain at fracture.

Plotted in figure 3 are stresses at fracture for uniaxial strain of bi-crystals containing symmetric tilt grain boundaries around the 〈001〉 tilt axis as a function of the tilt angle. In each case, the strain is applied in the direction normal to the interface. The solid symbols are from the molecular simulations using the REBO potential containing approximately 4000 atoms [27]. The stresses at 0° and 90° tilt angles correspond to bulk fracture of ideal diamond along the 〈110〉 and 〈100〉 directions, respectively. The dotted line is the maximum stress for fracture along the 〈111〉 direction. The maximum stresses from these simulations from inter-grain fracture (i.e. tilt angles between 0° and 90°) are lower than those for intra-grain fracture. Hence these simulations predict that, for ideal cleavage with only grain boundaries and no other defects such as cracks, inter-grain fracture should dominate over intra-grain fracture. The accuracy of this conclusion is questioned below.

The open triangles in figure 3 correspond to predictions of the critical fracture from the Orowan expression (2.6), while the open squares indicate predictions from the UBER expression (2.5), which, as shown above, is a scaled version of equation (2.6). While the Orowan expression significantly overestimates the maximum stress, the new scaled form from the UBER provides a reasonable approximation to the simulation results, especially for inter-grain fracture.

Plotted in figure 4 as the open triangles and open squares are the cleavage energies (corresponding to the left vertical axis) and direction-dependent Young's moduli (corresponding to the right vertical axis), respectively, that were used to calculate the maximum stresses plotted in figure 3 as the open symbols. Young's moduli were calculated directly from the REBO using small strain, elastic deformations along the respective directions. These cleavage energies are the sum of the grain boundary energy and the energies of the two surfaces that are created by the fracture. The solid circles are cleavage energies for each tilt orientation without a grain boundary (i.e. twice the surface energy for each orientation). For the REBO potential, the (110) surface energy is 3.4 J m^{−2}, which is substantially lower than the 6.5 J m^{−2} given from the DFT calculations discussed in the previous section. For the relaxed (111) surface, the REBO potential yields a surface energy of 5.4 J m^{−2}, which is relatively close to the DFT values of 6.5 J m^{−2}. Hence the REBO probably underestimates the maximum stress for ideal cleavage along the (110) plane. This is discussed more below.

Defining the inter-grain cleavage energy requires not only surface energies, but also the grain boundary energies prior to fracture. From a computational viewpoint, energies for short-period structures such as the twist grain boundary discussed above [26,30] can be calculated straight-forwardly using semi-empirical and DFT methods. However, accurately calculating energies for grain boundaries with long periods can be challenging [34,35]. This is because not only are the unit cells extended in the directions in the plane of the grain boundary, but these structures tend to have extended stress fields in the direction normal to the interface compared with short-period structures. This implies a need for relatively large unit cells that can capture stress effects on the interface energy. Grain boundaries with long periods also tend to be the higher energy interfaces, and therefore may be more prone to fracture than lower energy structures. Adding to this is the relatively large number of degrees of freedom associated with tilt and twist grain boundaries (5 d.f.) and the possibility of different closely iso-energetic bonding structures along an interface [35]. While large systems can be simulated using interatomic potentials such as the REBO, their accuracy and transferability may not be sufficient to produce qualitative fracture data. As discussed above, an example is the energy for the diamond (110) surface, which the REBO predicts to be smaller than the results of DFT calculations.

With these challenges in mind, Shenderova and co-workers applied a first principles disclination structural units model (FP-DSUM) to the calculation of symmetric tilt grain boundaries in diamond across an entire range of tilt angle for rotation about the 〈001〉 [37] tilt axis [37,38]. Based on work by Vitek and co-workers [39,40,41], the general concept is to model the structure of a grain boundary at any arbitrary angle as a linear combination of structures of special short-period grain boundaries that delimit the grain boundary of interest. The energy of the intermediate grain boundary can then be taken as a linear combination of the energies of the delimiting structures (properly weighted by the length of the respective grain boundaries), plus other terms that may be neglected or incorrectly accounted for by the combination of structural unit energies. In the disclination approach used previously [37,38], these terms are composed of an elastic energy that arises from the long-range interactions of different structural units, and a core energy term that accounts for intersections of different structural units.

The DSUM approach has several advantages compared with direct calculation of the energies of grain boundaries at many tilt angles. First, the delimiting grain boundaries can be taken as short-period structures that require relatively small unit cells compared with the structures at intermediate angles. The energies for these structures can therefore be calculated relatively easily with first principles methods. Second, in the disclination approach the core–core and elastic energy terms require only bulk elastic properties and one parameter that can be fitted to the energy of any carefully chosen intermediate angle grain boundary. Finally, the formalism can be used to analyse the contributions of each of the energy terms to the total interface energy as a function of tilt angle [42].

Plotted in figure 5*c* are the energies previously calculated from the FP-DSUM [37] for symmetric tilt grain boundaries in diamond rotated around the 〈100〉 axis as a function of tilt angle. Energies for the delimiting structures, which are from DFT calculations, are indicated by the solid squares. Plotted in figure 5*b* are the cleavage energies as a function of tilt angle for bulk cleavage (dashed line) and cleavage along grain boundaries (solid line) based on the DFT energies in the table 1. For comparison, the dotted line corresponds to intra-grain cleavage along the (111) plane. Plotted in figure 5*a* are the maximum stresses for inter-grain and intra-grain cleavage predicted from equation (2.5) using the first principles cleavage energies in figure 5*b* and Young's moduli from the table. As given in the previous section, intra-grain cleavage is preferred along the (111) plane. However, based on both cleavage energy and maximum stress, the first principles-based calculations predict that inter-grain cleavage is preferred for low tilt angles near a (110) orientation. Based on energy alone, this region is between 0° and approximately 48°, while the maximum stress criteria shorten this range of tilt angle to between approximately 5° and 38°. By contrast, the full simulations with the REBO potential suggest inter-grain fracture is preferred at all tilt angles.

Shenderova *et al.* [27] also performed similar molecular simulations of strained tilt grain boundaries using the REBO except that a 30 Å surface crack was initially introduced into the interface before the system was strained. These simulations contained approximately 50 000 atoms. As expected for a brittle material, the stresses at which the crack began to propagate in the simulations were two to three times smaller than those for ideal fracture. Fracture toughness was calculated from the simulations using the expression
3.1
where *σ*_{MD} is the stress in the simulations at which the crack begins to propagate, and *l* is the initial crack length. Toughness *K*_{c} was also calculated analytically via the expression
3.2
where *Y* _{ijk} and *E*_{ijk} are Young's modulus and surface energy, respectively, along the straining direction [*ijk*]. The fraction toughness values estimated from the dynamics simulations as a function of tilt angle for the symmetric tilt grain boundaries are plotted in figure 6 as the open squares, while the open triangles are estimates from equation (3.2) using moduli and surface energies from the REBO. Overall, the agreement between the two sets of values is reasonably good (less than about a 15% difference) and both sets of data show similar trends. Plotted as the dashed line is the fracture toughness from equation (3.2) as a function of tilt angle using the FP-DSUM data in figure 5. In this case, the REBO and DFT derived data are reasonably close to one another, and in good agreement with the value of 3.4 MN m^{−3/2} reported by Field & Freeman [43].

## 4. Conclusion

Discussed in this paper are several theoretical and computational approaches that have been used to better understand the fracture of both single-crystal and polycrystalline diamond at the atomic level. The studies discussed are not exhaustive but have instead been chosen to illustrate the different ways in which this problem has been approached, the conclusions (and their reliability) that have been reached by these methods, and how these methods can be effectively used together. For example, a modification of the traditional Orowan fracture criterion was derived from the UBER, and it was shown that this new expression reproduces the results of molecular simulations using an analytic REBO potential energy function. It was also shown that the limitations of the REBO in terms of its accuracy for surface energies compared with first principles results can lead to incorrect results, in this case regarding the stress at fracture. With the UBER-derived fracture criteria validated within the classical simulations, it could be parametrized with first principles data to lead to presumably more accurate conclusions regarding inter- versus intra-grain fracture. In particular, it was found that for ideal fracture of single-crystal diamond the (111) surface has the lowest maximum stress, but that inter-grain fracture has an even lower maximum stress. These conclusions are consistent with experimental results regarding diamond, and lend new insights into the origin and details of these observations.

This combination of methods has been used to examine ideal fracture for a very limited number of possible grain boundary orientations, and therefore there is much more in terms of computation to be done to fully characterize fracture of diamond from an atomic perspective. Similarly, the role of bond breaking and forming with associated changes in bonding type in the cleavage process has not been well characterized, nor has the role of processes such as hydrogen and oxygen termination of unsatisfied bonds created during fracture in the fracture process. Although not discussed above, there have also been modelling studies that have examined the role of nitrogen near grain boundaries in the fracture process and strength of UNCD [30,44–47]. Much more can be done in terms of theory and computation to understand the role of these types of lattice defects in the fracture of different forms of diamond-based materials.

## Acknowledgements

We wish to acknowledge the US Office of Naval Research, which originally supported our research efforts in this area.

## Footnotes

One contribution of 19 to a theme issue ‘Fracturing across the multi-scales of diverse materials’.

- © 2015 The Author(s) Published by the Royal Society. All rights reserved.