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_xN) (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,

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,

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.

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 ultra-soft 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) values which are the averages of the combined Reuss and Voigt values.

For Voigt bulk, B_V and shear, G_V moduli values are calculated as:

respectively, where C_ij are the stiffness constants.

For Voigt bulk, B_R and shear, G_R moduli values are calculated as:

respectively, where S_ij are the compliances.

The bulk and shear values are averaged to give the Hill values as:

The respective Hill values of the Young modulus, E_H and Poisson ratio, n_H are obtained from the expressions:

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.

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₇B, C₃B (or C₆B₂) and C₅B₃ have no negative phonon frequencies at the G 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₇B, C₃B and C₅B₃. 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₅B₃ 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 p 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₇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₅B₃.

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₃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₅B₃ and C₇B. The results are shown in Table 4. Table 4 shows that rhombohedral C₃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 m [SG number index, 227]. This structure has a face-centred cubic Bravais lattice with a Point Group, m 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₇B structure with the Hermann Mauguin Space Group, P 3m  [space  group  number  index, 215], which has a primitive cubic Bravais lattice with the Point Group, 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₇B, in short.

Thirdly, a C₃B (C₆B₂) 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₃B.

Fourthly, a C₅B₃ 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₅B₃.

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₇B, r-C₃B and o-C₅B₃, it is C₇B that seems to be the most promising in terms of hardness. The electronic energy band structure of cubic C₇B was calculated in order to characterize some of its electronic properties.

The k-point path used in this calculation, with reference to Figure 1, was: . Figure 4 shows the energy band diagram of c-C₇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₇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,  and plotted as a function of the applied pressure, as presented in Figure 5. 

The inequalities (Equation 14) are satisfied by the C₃B (or C₆B₂) 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:

The inequality (Equation 15) is satisfied by the C₅B₃ results in Table 5, it follows that o-C₅B₃ 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-yB_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₃B and C₅B₃ 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₅B₃ 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₇B as expected. Interestingly, the same condition is observed for C₃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³ 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³ 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, 2020).

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₇B (c-C₇B), rhombohedral C₃B (r-C₃B) and orthorhombic C₅B₃ (o-C₅B₃) 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 increasing boron content. As seen in Table 6, the materials under study are all brittle with diamond being the most brittle. C₃B and C₅B₃ 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₅B₃ 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₇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₇ is therefore a hole-type 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_xC_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₇B, will qualify for viable niche industrial applications, through experimental evaluation.

The authors thank Institutional facilities of both the Copperbelt and Mukuba Universities used during the course of this research work.