Stability and elastic anisotropy of diamond related C8-yBy materials

A number of potentially super-hard materials were examined using ab-initio methods. Low Gibbs free energy polymorphs of diamond-like materials for y = 0 to 7 in the stoichiometric type C8-y By, were identified at absolute zero of temperature. These were proposed as possible super-hard materials with useful applications. The materials with y = 0 to 3, that is, diamond (C), cubic C7B (c-C7B), rhombohedral C3B (r-C3B) and orthorhombic C5B3 (o-C5B3) were found to be dynamically and mechanically stable. A diamond standard was used as a stable comparison. Results of their bulk modulus calculations suggest that these materials were potentially super-hard in character. Systematic trends were established, the hardness was observed to reduce with increasing boron content. The materials under study were all determined as being brittle with diamond being the most brittle, C3B and C5B3 are the least brittle with B/G values of 1.32. Of the materials studied, diamond was determined to have the lowest degree of elastic anisotropy with a Universal Elastic Anisotropy Index of only 0.041 while C5B3 had the highest anisotropy of 1.160, making it the most susceptible to micro-cracks. Our electronic band structure studies of c-C7B, which was predicted to be the hardest in the C8-y By system after diamond, showed that the top of the valence band was about 1.7 eV above the Fermi level with a band gap between the valence and conduction bands, making c-C7B a hole-type conductor having a likely increase in conductivity with increased applied hydrostatic pressure. 
 
 Key words: Phase stability, elastic anisotropy, ultra-hard material.


INTRODUCTION
Super-hard materials have important applications in high speed machining tools for cutting and drilling as well as abrasives and wear-resistant coatings because of their strength. Diamond is extremely hard but at high temperatures and in the presence of oxygen it becomes unstable due to oxidation reactions (John et al., 2002). It is not suitable for machining steel and other alloys of iron because of its redox reactions with iron and some other metallic elements (Nassau and Nassau, 1979) at temperatures exceeding 80 K. Cubic boron nitride (c-BN) is a super abrasive which is more thermally stable than diamond and is better suited for machining steel. However, despite the high oxidation resistance temperature and high chemical inertness of c-BN, it is only about half as hard as diamond (Singh, 1986). The possible replacement of both diamond and c-BN by better performing super abrasives has motivated a lot of researcher interest.
An important frontier in high-pressure science and technology research has been the synthesis of novel ultra-or super-hard materials. Diamond anvil cells (DAC) are reported (Stavrou et al., 2016) to have been successfully used in synthesizing such materials under high pressure and temperature conditions. Covalent bonds are strong and highly resistant to plastic and elastic deformations (Zhao et al., 2016). The predominant atomic binding mechanism in super-hard materials is therefore covalent bonding. The constituent atomic elements in these materials are often light, like C, B, O, and N (Habanyama et al., 2018) which are chemically suited for covalent bonding (Hu et al., 2016). Since covalent bonds are generally directional and short, the light elements are able to form highly shear resistant 3-dimentional networks.
A method of reducing diamond"s susceptibility to oxidation and reaction with ferrous alloys is to dope it with boron. However, the highest possible doping concentration of boron in boron-doped diamond (BDD) is very low. The synthesis of new diamond-like compounds of carbon and boron is being considered as an effective method of achieving a much higher boron content. There have been recent developments in computational discovery of super-hard materials (Kvashnin et al., 2019) and ab-initio calculations of elastic constants (Mazhnik and Oganov, 2019;Niu et al., 2019). Recent theoretical work has also been carried out on diamond-like boron carbon nitrides (d-BC x N) (He et al., 2019;Gao et al., 2018).
The Vickers hardness (H V ) of ultra-hard materials is known to be more than 40 GPa (Solozhenko and Gregoryanz, 2005). Materials with super-hard characteristics are known to have a value of the bulk modulus that exceeds 250 GPa (Lowther, 2000). Both the bulk and shear moduli have been used as scales to estimate the level of hardness of materials (Clerc, 1999).
A single crystal seldom exhibits a totally isotropic elastic response, essentially all known crystals exhibit some form of elastic anisotropy, meaning that the elastic moduli are generally dependent on the different crystal orientations. The degree of anisotropy in the properties of materials is very important in their application. The generation of micro-cracks and lattice distortions in materials is often related to the elastic anisotropy. Ledbetter and Migliori (2006) point out the effects of elastic anisotropy in dislocation dynamics, phase transformations and other crystal phenomena.
Energetic or thermodynamical stability is determined by the relative value of the system"s Gibbs free energy, (1) The internal energy is E, volume is V, pressure is P, temperature is T and the system"s entropy is S. The most stable compound phase maximizes the chemical bond strength and has the lowest value of G. Our calculations involve ground state structures of compounds determined under a thermal limit conditions (that is, neglecting zero point motion) at zero-temperature (T = 0 K). In this case, the minimum Gibbs free energy is equivalent to the minimum enthalphy, H where, (2) In this work, we study the energetic, dynamical, mechanical and anisotropic properties of diamond-like materials for y = 0 to 7 in the stoichiometric type C 8-y B y .

COMPUTATIONAL METHODS
The method of calculation applied in this work is based on the density functional theory (DFT) (Kohn and Sham, 1965). The Quantum Espresso (Giannozzi et al., 2009) software package implementation of DFT was adopted. The exchange-correlation interaction between the electrons was modeled using the Generalized Gradient Approximation as parameterized by Perdew, Burke and Ernzerhof (GGA-PBE) (Perdew et al., 1996). The ultrasoft pseudo-potential method was used to calculate the interaction between the electrons and the ion cores. The cut-off energy used for the plane-wave functions was 50 Ry, sampling of the k-point mesh in the Brillouin zone was 666 Monkhorst Pack (Monkhorst and Pack, 1976). The convergence threshold of the self-consistent field (SCF) was within 10 -3 eV/atom. The initial structure consisted of eight atoms in a unit with boron or carbon atoms placed at diamond lattice positions in a ratio that represented a particular stoichiometric value of y. The lattice positions and atomic types at each position were stored in matrix arrays. A C++ template library for linear algebra called Eigen (http://eigen.tuxfamily.org/) was used to randomise the atomic placement. Variable-cell dynamics calculations were performed on the starting cell to optimize the atomic position geometry and cell parameters using the Quantum Espresso (Giannozzi et al., 2009) package. This was achieved by performing some relaxation operations which allowed the positions of the atoms to adjust in accordance with the relative forces between them while allowing the unit cell to vary; equilibrium atomic structures were thereby achieved. The space-group and symmetry operations of these equilibrium configurations were identified using the utility program, SGROUP (Yanchitsky and Timoshevskii, 2001).
The procedure of randomising the atomic placement (using Eigen) in the starting cell, variable-cell relaxation (using Quantum Espresso) and identification of the final converged configurations (using SGROUP), was carried out several times for each value of y in C 8-y B y where y = 0 to 7. The crystal structures that resulted from the same value of y had similar stoichiometry and were polymorphs. Calculations were carried out on these structure to determine the configurations with the lowest value of the Gibbs free energy, G for each value of y in C 8-y B y where y = 0 to 7, at absolute zero of temperature. These eight configurations were then subjected to dynamical stability tests.
A material"s vibrational spectrum needs to be analyzed in order to establish its dynamical or vibrational stability. The criterion for this aspect of stability is the nonappearance of imaginary phonons, that is, frequencies of vibration that are not positive in the dispersion calculations on phonons or lattice vibrations at reciprocal lattice vectors in the Brillouin zone. Phonon dispersion frequencies were obtained for specific modes of vibration using the Quantum Espresso implementation of the density functional perturbation or "linear response" theory (Baroni et al., 1987;Giannozzi et al., 1991).
The structures which were found to be dynamically stable were constructed and visualized by the Xcrysden software package (Kokalj, 2003) after which they were tested for mechanical or elastic stability which considers the second order elastic constants of a material. Elastic stiffness and compliance constants were calculated using the elastic package (Golesorkhtabar et al., 2013).
Three types of averaging calculations were carried out to determine the shear and bulk moduli of the various compounds: the Voigt (1928) calculation where the strain is taken to be uniform, Reuss and Angew (1929) calculation where the stress is taken to be uniform and Hill (1963)  The elastic constant results presented for all materials in this work are Hill-averaged, unless indicated otherwise.
The bulk modulus for crystals that exhibit cubic symmetry is the same in all directions. Cubic elastic anisotropy is therefore determined by the shear anisotropy alone. There are both bulk and shear anisotropic contributions in crystals that are not cubic. An appropriate way to quantify the combined bulk and shear degree of anisotropy is by using the universal elastic anisotropy index, A U (Ranganathan and Ostoja-Starzewski, 2008) which is given by, . ( If a single crystal is isotropic in both bulk and shear then, B V = B R and G V = G R , making A U equal to zero. In cubic crystals B V = B R but G V  G R . The degree of anisotropy is determined by the extent of the departure of A U from zero.

RESULTS
The crystal structures that resulted from the same value of y had similar stoichiometry and were polymorphs. Calculations were carried out on these structure to determine the configurations with the lowest value of the Gibbs free energy, G for each value of y in C 8-y B y where y = 0 to 7, at absolute zero of temperature. The resulting eight configurations are listed in Table 1. The structures were then tested for dynamical stability. Our phonon calculation results at the Brillouin zone center (the Г point) are also presented in Table 1. The value, y = 0 represents a diamond standard, used as a stable comparison. The calculations of the results shown in Table 1 covered 24 frequencies for each material but the table only shows the lowest and highest 6 results, all other frequencies were positive. This table shows that the materials with y = 0 to 3, that is, diamond (C), C 7 B, C 3 B (or C 6 B 2 ) and C 5 B 3 have no negative phonon frequencies at the  point of the Brillouin zones, indicating possible dynamical stability. All the other structures showed some negative frequencies indicating dynamical instability. Further phonon calculations, to complete the test for dynamical stability, were therefore only carried out on C 7 B, C 3 B and C 5 B 3 . Further phonon calculations for diamond are not included here since the stability of diamond is well established. Figure 1 shows some high symmetry lines and points in the first Brillouin zone of the orthorhombic lattice structure.
Phonon calculations were carried out at all the eight high symmetry points as shown in Figure 1, for the orthorhombic C 5 B 3 structure, the results are presented in Table 2. The calculations in Table 2 covered 24 frequencies for each symmetry point, but as in Table 1, we only picked the lowest and highest 6 results, all other frequencies were positive. Each frequency is presented in units of both THz (the upper value) and cm -1 (the lower value). The k-space coordinates of the Brillouin zone indicated in Table 2 are in integral units of a reciprocal lattice multiplying factor of  divided by the corresponding lattice parameters. It is seen in Table 2 that there are no negative phonon frequencies at any of the eight high symmetry points of the Brillouin zone; this structure is therefore dynamically stable. In Table 2, acoustic mode vibrations correspond to the low frequencies while the higher frequencies indicate vibrations in the optical mode.  If all sides of the first Brillouin zone of the orthorhombic lattice were drawn with equal lengths then Figure 1 would represent the Brillouin zone of a cubic lattice. Phonon calculations were carried out at all the eight high symmetry points as shown in Figure 1, for the cubic C 7 B structure, the results are presented in Table 3. The results in Table 3 show the lowest and highest 6 frequencies, with all other frequencies being positive, as explained in the case of C 5 B 3 . It is seen in Table 3 that there are no negative phonon frequencies at any of the eight high symmetry points of the Brillouin zone; this structure is therefore also dynamically stable.
Phonon calculations for C 3 B were carried out at eight k-   point coordinates in the Brillouin zone of the rhombohedral structure of the trigonal system, in a similar way as was done for C 5 B 3 and C 7 B. The results are shown in Table 4. Table 4 shows that rhombohedral C 3 B is dynamically stable. The details of the crystal systems for the materials that were found to be dynamically stable are shown in the following.
Firstly, the carbon-cubic diamond structure which has the Hermann Mauguin Space Group, Fd 3 m [SG number index, 227]. This structure has a face-centred cubic Bravais lattice with a Point Group, m 3 m. The "F" in the space group signifies that it is a face-centred lattice type. Diamond is already well studied and is only being used here as a standard to which the other materials will be compared.
Secondly, a C 7 B structure with the Hermann Mauguin Space Group, P 4 3m [space group number index, 215], which has a primitive cubic Bravais lattice with the Point Group, 4 3m. The "P" in the space group signifies that it is a primitive lattice type. We will refer to this cubic material as c-C 7 B, in short.
Thirdly, a C 3 B (C 6 B 2 ) structure with the Hermann Mauguin Space Group, R3m [space group number index, 160], which has a rhombohedral structure of the trigonal Bravais lattice crystal system with the Point Group, 3m. The "R" in the space group signifies that it is a rhombohedral lattice type. We will refer to this rhombohedral material as being r-C 3 B.
Fourthly, a C 5 B 3 structure with the Hermann Mauguin Space Group, Imm2 [space group number index, 44], which has a body-centred orthorhombic Bravais lattice with the Point Group, mm2. The "I" in the space group signifies that it is a body-centred lattice type. We will refer to this orthorhombic material as being o-C 5 B 3 . Table 3. Results of phonon calculations for the cubic C 7 B structure at , Z, T, Y, S, X, U and R high symmetry points of the Brillouin zone.  1.0,1.0,1.0 The k-space coordinates indicated are in units of  divided by the corresponding lattice parameters.
Crystal diagrams of these four structures were constructed and visualized by the Xcrysden software package (Kokalj, 2003); they are as shown in Figure 2a and b followed by Figure 3a and b, respectively.
Previous theoretical studies (Sun et al., 2001;Nozaki and Itoh, 1996), on covalent bonding dominated compound structures of the ternary B-C-N system (which includes the binary B-C system), indicated that no B-B bonds were expected to exist in the crystal structures. This is because B-B bonds would effectively increase the total energy of the system hence making the structures less stable. It can be seen in Figures 2 and 3 that our diagrams are consistent with the results of these previous studies.
Elastic stiffness and compliance constants were calculated using the Elastic package (Golesorkhtabar et al., 2013) for the four materials that were identified as being dynamically stable; the results are shown in Table  5.
The materials in Table 5 were determined to be mechanically stable, as discussed subsequently. This qualified these materials for the determination of their bulk and other moduli, Hill values of these moduli are shown in Table 6. The table also shows the B/G ratio and the Poison ratio. Table 7 presents the Voigt and Reuss bulk and shear moduli together with the universal elastic anisotropy index, calculated using Equation 11. Table 7 shows that of the three predicted materials, that is, c-C 7 B, r-C 3 B and o-C 5 B 3 , it is C 7 B that seems to be the most promising in terms of hardness. The electronic energy band structure of cubic C 7 B was calculated in order to characterize some of its electronic properties.   1.0,1.0,1.0 The k-space coordinates indicated are in units of  divided by the corresponding rhombohedral lattice parameters.
The k-point path used in this calculation, with reference to Figure 1, was: Figure 4 shows the energy band diagram of c-C 7 B under zero applied pressure where the origin of the energy axis is placed at the Fermi level, E F.
The electronic band structure in Figure 4 shows that the highest valence energy band of the material crosses over the Fermi level up to about 1.7 eV above it. On the other hand, Figure 4 also shows a band gap between the top of the valence band and the bottom of the conduction band. This type of band structure configuration is indicative of hole-type conductivity. The material, c-C 7 B is therefore a hole-type conductor.
The evolution of the Fermi energy, E F at varying applied hydrostatic pressures from 0 to 2500 kbars was calculated relative to the Fermi energy at zero pressure, ) 0 ( F E and plotted as a function of the applied pressure, as presented in Figure 5. It is seen in Figure 5 that the relative Fermi energy increases monotonically with increased applied pressure. However, the curvature in the graph suggests that saturation would eventually be arrived at with further increase in applied pressure. An increase in Fermi energy with applied pressure has been associated with increased conductivity (Ghosh et al., 2016).

DISCUSSION
Materials with super-hard characteristics are known to have a value of the bulk modulus that exceeds 250 GPa (Lowther, 2000). Table 6 indicates that our materials with y = 0 to 3 could potentially have super-hard characteristics.
The original conditions for mechanical stability were proposed by Born-Huang (1954). In the current work, we adopt the modified conditions of Mouhat and Coudert (2014). A necessary although insufficient Born condition is that there should be no negative diagonal elements, that is, . It is seen in Table 5 that all the diagonal elements are positive for the compounds studied.
A second condition which is necessary but insufficient for Born stability is, .
(12) Table 5. Independent stiffness constants, C ij and compliances, S ij for C 8-y B y materials where y = 0 to 3.

Material Crystal Structure
Stiffness matrix elements, C ij (GPa) and compliances, S ij (10 -5 GPa -1 )   We find from Tables 5 that all the elastic constants obtained for the different materials satisfy the condition given by inequality (Equation 12).
The cubic system only has 3 independent elastic constants and as a sufficient condition for mechanical stability, these constants have to satisfy the inequalities: ; .
The inequalities (Equation 13) are satisfied by the diamond and C 7 B results in Tables 5, it follows that these two cubic structures are mechanically stable.
The rhombohedral class has 6 independent elastic constants. Sufficient Born conditions for this class as modified by Mouhat and Coudert (2014), are given by the inequalities: (14) , plotted against the applied pressure, for c- The inequalities (Equation 14) are satisfied by the C 3 B (or C 6 B 2 ) results in Tables 5, it follows that this rhombohedral structure is mechanically stable.
The orthorhombic system has 9 independent elastic constants. A sufficient Born condition for this system (Born and Huang, 1954) is given by the inequality: . 15) The inequality (Equation 15) is satisfied by the C 5 B 3 results in Table 5, it follows that o-C 5 B 3 is mechanically stable. Figure 6 is a graph of the values of the Young, shear  bulk moduli plotted against the values of y in C 8-y B y materials. The three moduli are seen to reduce in magnitude as the boron content increases. The brittleness of a material is reflected by the ratio B/G. As reported by Pugh (1954), the critical value of this ration is B/G=1.75, for ductile materials B/G>1.75. The ratio is less than 1.75 for brittle materials. As seen in Table 6, the four materials under study are all brittle with diamond being the most brittle. C 3 B and C 5 B 3 are the least brittle with a B/G value of 1.32.
The Young modulus relates a stress to the resultant strain in the direction of the stress. The lateral and axial strains are related by the Poison ratio, which gives an indication of the stability of a material in relation to resistant shear and also the magnitude of the change in volume of a crystal during uniaxial deformation. The plasticity of a material increases with an increase in the Poisson ratio. The results in Table 6 indicate that the plasticity of the materials generally increases with increased boron content. Table 7 shows that diamond has the lowest degree of anisotropy with a Universal Elastic Anisotropy Index of only 0.041. C 5 B 3 has the highest anisotropy (A U = 1.160), this material is therefore likely to be the most susceptible to micro-cracks. The table shows that B V = B R for the two cubic structures, C and C 7 B as expected. Interestingly, the same condition is observed for C 3 B which is not cubic, indicating that this material is also isotropic in terms of the bulk modulus alone.
There are a combination of mechanisms responsible for the observed reduction in magnitude of the three moduli of elasticity as the boron content increases. These same mechanisms are responsible for the general increase in plasticity and decrease in brittleness with increased boron content. The underlying mechanism is the sp 3 hybridization in the orbitals of the carbon atoms.
The creation of hybrid orbitals in carbon atoms during material synthesis is illustrated in Figure 7. Figure 7a shows the pairing of 1s and 2s spins, while 2p electron spins are not paired and are available for covalent bonding. A 2s electron can be excited, during material synthesis, to the p level resulting in four similar hybridized (mixed) sp 3 orbitals as shown in Figure 7b. It is these orbitals that form the strongest covalent bonds in nature; they are highly resistant to plastic and elastic deformations. In general, the more the hybrid C-C bonds you have in a boron/carbon material, as opposed to C-B bonds, the higher the magnitude of the three moduli of elasticity (Sun et al., 2001;Nozaki and Itoh, 1996). However, there are other factors that may come into play. The shorter the bonds and the denser they are the higher the magnitude of the three moduli, the crystal structure of the materials affects the length and density of the bonds. A combination of these factors affects the plasticity and brittleness of the materials with increased boron content. Factors like the iconicity of the bonds also play a part. Some theoretical models have been developed to relate these fundamental factors to the hardness of materials (Bao et al., 2018(Bao et al., , 2020.

Conclusion
Low Gibbs free energy polymorphs of diamond-like materials of the type C 8-y B y where y = 0 to 7, were identified at absolute zero of temperature. It has been determined that the four materials: diamond (C), cubic C 7 B (c-C 7 B), rhombohedral C 3 B (r-C 3 B) and orthorhombic C 5 B 3 (o-C 5 B 3 ) are all both dynamically and mechanically stable. The bulk modulus results in Table 6 indicate that these compounds are potentially super-hard; the values of their bulk moduli being greater than 250 GPa (Lowther, 2000). Figure 6 suggests that the hardness reduces with Samukonga et al. 27 increasing boron content. As seen in Table 6, the materials under study are all brittle with diamond being the most brittle. C 3 B and C 5 B 3 are the least brittle with B/G values of 1.32. Table 7 shows that diamond has the lowest degree of anisotropy with a Universal Elastic Anisotropy Index of only 0.041 while C 5 B 3 has the highest anisotropy (A U = 1.160), which makes the latter material most susceptible to micro-cracks. Our electronic band structure calculations for cubic C 7 B, which was predicted to be the hardest in the C 8-y B y system after diamond, indicate that the highest valence energy band of the material crosses over the Fermi level up to about 1.7 eV above it. On the other hand, there is a band gap between the top of the valence band and the bottom of the conduction band. The material, c-BC 7 is therefore a holetype conductor, similar to heavily boron doped diamond (Lee and Pickett, 2004). This is understandable because boron is an acceptor with one less electron in the B x C y valence matrix, hence functioning as a hole-dopant. When carrying out a cost-benefit ratio analysis for the industrialization of the materials studied in this article, it should be noted that our results and the available literature indicates that the synthesis of materials with hardness exceeding that of diamond is highly unlikely, if not impossible (Solozhenko and Godec, 2019). The economical approach should therefore be to synthesize or design materials for niche applications, materials that are harder than cubic boron nitride and more useful than diamond in terms of thermal and chemical stability. Such novel materials have great prospects for the creation of new technologies required for emerging applications and their world market is inexhaustible (Solozhenko and Godec, 2019). It is hopeful that the materials studied during this project, particularly cubic C 7 B, will qualify for viable niche industrial applications, through experimental evaluation.