Prediction of Novel Ultrahard Phases in the B–C–N System from First Principles: Progress and Problems

The modern synthesis of superhard and, especially, ultrahard phases is a fascinating area of research that could lead to the design of new, industrially important materials. Computational methods built within the well-established quantum mechanics framework of density functional theory (DFT) play an important role in the search for these advanced materials and the prediction of their properties. The close relationship between the physical properties of carbon and boron nitride has led to particular interest in the B–C–N ternary system, characterized by the small radii of the elements, resulting in short interatomic distances and reduced volumes—the parameters being ‘recipes’ for very high hardness in three-dimensional structures. The purpose of this review is to provide a brief outline of recent developments and problems in predicting novel ultrahard carbon allotropes as well as binary and ternary compounds of the B–C–N system with particular emphasis on the analysis of the models used to evaluate the hardness of the theoretically predicted structures.


Introduction
Historically, research in the field of ultrahard materials (usually defined as having Vickers hardness H V ≥ 80 GPa) was initiated a century ago, triggered by the increasing industrial application of diamond, the hardest (H V up to 120 GPa) known material. However, the use of diamond in mining and tooling raised problems due to the high cost of natural diamond, on the one hand, and its relatively low stability at even moderate operating temperatures, on the other hand. The first problem was solved by the development of the high-pressure synthesis of diamond [1], but its high reactivity with oxygen and ferrous metals remained a problem that required further efforts to find reliable substitutes for diamond.
An interim solution was the synthesis of cubic boron nitride (cBN) [2], which is half as hard as diamond, but has much higher thermal and chemical stability. It should be noted that in terms of electronic structure, BN is equivalent to 2C (vide infra). The close relationship between the physical properties of carbon allotropes and BN polymorphs has facilitated the search for ultrahard phases in the B-C-N ternary system constituted by light elements with small radii, resulting in short interatomic distances and reduced volumes-all the parameters being 'recipes' for high strength and hardness.
Until very recently, diamond was the only known material that was ultrahard. However, in 2001, cubic BC 2 N, a ternary compound that is halfway between diamond and BN in composition, was synthesized [3]. Vickers hardness of 76(4) GPa [4] makes it the second member of the ultrahard phases family. In 2009, the third ultrahard phase, diamond-like BC 5 , was discovered [5], which possesses Vickers hardness of 71 (8) GPa, an unusually high for superhard materials fracture toughness (~10 MPa·m

Elements
Carbon and boron (in particular, their dense allotropes, i.e., diamond and γ-B12 [19]) are the hardest elements. We shall not discuss boron here (the problem of searching for its superhard forms was already addressed in a recent review [20]), but will rather concentrate on carbon. The main research efforts with respect to carbon are focused on the search for new, dense allotropes with mechanical properties close to those of diamond. Diamond is known to exist in two forms: the cubic one that belongs to space group Fm−3m (No. 227), and a rare hexagonal form called 'lonsdaleite' (space group P63/mmc, No. 194). Recently, however, the existence of lonsdaleite as a discrete phase has been questioned, and it has

Elements
Carbon and boron (in particular, their dense allotropes, i.e., diamond and γ-B 12 [19]) are the hardest elements. We shall not discuss boron here (the problem of searching for its superhard forms was already addressed in a recent review [20]), but will rather concentrate on carbon. The main research efforts with respect to carbon are focused on the search for new, dense allotropes with mechanical properties close to those of diamond. Diamond is known to exist in two forms: the cubic one that belongs to space group Fm−3m (No. 227), and a rare hexagonal form called 'lonsdaleite' (space group P6 3 /mmc, No. 194). Recently, however, the existence of lonsdaleite as a discrete phase has been questioned, and it has been interpreted as a cubic diamond dominated by extensive twins and stacking faults [21].
It should be noted that information on 3-periodic carbon allotropes extracted from the scientific literature is gathered and indexed in the "SACADA" database [22]. Currently, it counts 524 unique carbon allotropes. From the structural point of view, and using chemistry terms, diamond is characterized by tetrahedral (sp 3 ) carbon, as in the methane gas molecule CH 4 , and its crystal structure demonstrates three-dimensional stacking of corner-sharing C4 tetrahedra. In contrast, graphite is characterized by a layered structure where carbon is sp 2 hybridized, as in the ethene gas molecule C 2 H 4 . While diamond is a large-band-gap insulator with E gap ≈ 5 eV, graphite is a semiconductor with a very small gap. Mixed sp 2 -sp 3 hybridizations carbon was recently reported for a new, stable metallic allotrope, hex-C 18 (called 18H carbon), which belongs to space group P6/mmm (No. 191) [23].
Carbon in a linear triatomic C-C-C arrangement is found in carbon suboxide C 3 O 2 , where carbon is sp 2 hybridized. In the solid state, the linear configuration is kept for carbon suboxide, which is a molecular solid where weak (van-der-Waals-like) interactions prevail between separate molecules. The rarely occurring tricarbon C 3 molecule was observed in interstellar space, mainly in the tails of comets (e.g., Hale-Bopp, C/1995O1), and experimentally identified by spectroscopic measurements [24]. Recently [25], we considered linear C-C-C in two novel structures: (i) rh-C 3 (or hex-C 9 ) (space group R−3m, No. 166) based on rhombohedral sodium azide characterized by the presence of a linear N 3 fragment (Figure 2a), and (ii) hex-C 6 (space group P6 3 /mmc, No. 194) derived from lonsdaleite through the insertion of one extra carbon atom along the c-axis at 2d (2/3,1/3, 1 4 ) Wyckoff position (Figure 2b). From geometry optimization to the energy ground state within DFT (cf. Appendix A), energies and the energy-derived quantities were found within the range of diamond and lonsdaleite. been interpreted as a cubic diamond dominated by extensive twins and stacking faults [21].
It should be noted that information on 3-periodic carbon allotropes extracted from the scientific literature is gathered and indexed in the "SACADA" database [22]. Currently, it counts 524 unique carbon allotropes. From the structural point of view, and using chemistry terms, diamond is characterized by tetrahedral (sp 3 ) carbon, as in the methane gas molecule CH4, and its crystal structure demonstrates three-dimensional stacking of corner-sharing C4 tetrahedra. In contrast, graphite is characterized by a layered structure where carbon is sp 2 hybridized, as in the ethene gas molecule C2H4. While diamond is a large-band-gap insulator with Egap ≈ 5 eV, graphite is a semiconductor with a very small gap. Mixed sp 2 -sp 3 hybridizations carbon was recently reported for a new, stable metallic allotrope, hex-C18 (called 18H carbon), which belongs to space group P6/mmm (No. 191) [23].
Carbon in a linear triatomic C-C-C arrangement is found in carbon suboxide C3O2, where carbon is sp 2 hybridized. In the solid state, the linear configuration is kept for carbon suboxide, which is a molecular solid where weak (van-der-Waals-like) interactions prevail between separate molecules. The rarely occurring tricarbon C3 molecule was observed in interstellar space, mainly in the tails of comets (e.g., Hale-Bopp, C/1995O1), and experimentally identified by spectroscopic measurements [24]. Recently [25], we considered linear C-C-C in two novel structures: The elastic properties of 'tricarbon' allotropes rh-C3 and hex-C6 were determined by performing finite distortions of their lattices. The elastic constants Cij were derived from the strain-stress relationship. Indexes i and j represent directions: when i = j, the elastic constants correspond to the application of unidirectional stress (as in the case of C33 elastic constant discussed below), and when i ≠ j, the elastic constants are relevant to applying shear stress.  The elastic properties of 'tricarbon' allotropes rh-C 3 and hex-C 6 were determined by performing finite distortions of their lattices. The elastic constants C ij were derived from the strain-stress relationship. Indexes i and j represent directions: when i = j, the elastic constants correspond to the application of unidirectional stress (as in the case of C 33 elastic constant discussed below), and when i = j, the elastic constants are relevant to applying shear stress. For both allotropes, all calculated C ij values are positive, and their combinations obey the rules pertaining to the mechanical stability of the phases. The bulk (B V ) and shear (G V ) moduli were calculated from the elastic constants following the Voigt method [26], based on a uniform strain.

GPa MPa·m
As has been reported earlier [36], in the case of ultrahard compounds of light elements, the thermodynamic model shows surprising agreement with available experimental data. Moreover, its use is preferable in the case of hybrid dense carbon allotropes, for which the Lyakhov-Oganov model gives underestimated hardness values, whereas the empirical models are not reliable. For this reason, both new 'tricarbon' allotropes should be considered as ultrahard phases.
Although the hardness and elastic moduli of rh-C 3 and hex-C 6 are somewhat lower than those of diamond, strong anisotropy was found for both 'tricarbon' allotropes, with exceptionally large C 33 values along the hexagonal c-axis, i.e., C 33 = 1636 GPa for rh-C 3 ( Figure 2a) and C 33 = 1610 GPa for hex-C 6 (Figure 2b), exceeding the corresponding value for lonsdaleite, C 33 = 1380 GPa. The Vickers hardness predicted using four theoretical models (Tables 1 and 2) points to slightly lower H V values for these new carbon allotropes compared to diamond (both cubic and hexagonal), but much higher than the hardness of the vast majority of recently predicted carbon allotropes, such as C 14 , C 16 , C 24 , C 36 , etc. [42][43][44][45]. Thus, both rh-C 3 and hex-C 6 have exceptional mechanical properties and can be considered as prospective ultrahard phases [46].
The dynamic stability of the 'tricarbon' allotropes was confirmed by phonon calculations. All frequencies are positive, with the particular feature of a gap in the highest-frequency optical phonon domain, not observed in lonsdaleite, and caused by the rigidly aligned C 3 unit. A remarkable consequence of the presence of two carbon hybridizations (sp 2 and sp 3 ) is the occurrence of a metallic character in the electronic band structure of both 'tricarbon' allotropes, similar to that previously observed for hexagonal C 18 [23] and monoclinic C 12 [47] characterized by mixed sp 2 -sp 3 and sp 1 -sp 2 hybridizations, respectively.
Another novel ultrahard carbon allotrope, body-centered tetragonal C 6 (space group I−4m2, No. 119) presenting mixed sp 2 /sp 3 hybridizations, has been proposed very recently via a crystal chemistry approach and studied for the ground state structure and stability using DFT calculations [29]. Since C4 tetrahedra are in-plane stacked with corner-sharing and connected out-of-plane with C-C trigonal carbon (Figure 2c), a close relationship with so-called 'glitter', a hypothetical dense carbon network invented in 1994 [48], is apparent, and thus the new allotrope was named 'neoglitter'. Besides the mechanical stability (positive values of elastic constants and their combinations), 'neoglitter' is also dynamically stable, as it follows from its phonon band structure. The novel allotrope reveals exceptional mechanical properties, i.e., very high hardness and elastic moduli (see Tables 1 and 2), being conductive due to the metallic-like electronic structure, which is mainly caused by the itinerant role of trigonal carbon π-electrons.
Since the Vickers hardness calculated in the framework of the thermodynamic model exceeds 80 Gpa for the three novel carbon allotropes described above, they all should be attributed to the family of ultrahard phases.

Binary Compounds
The hardness of dense compounds of the binary B-N system-cubic BN [2], rhombohedral B 13 N 2 [49,50] and tetragonal B 50 N 2 [51]-does not exceed H V = 62 Gpa for single-crystal cBN [4], i.e., they all belong to the group of superhard phases. Special mention should be given to nanocrystalline cBN [52], with Vickers hardness up to 85 Gpa [53], mainly due to the Hall-Petch effect, i.e., nanosize effect, which restricts dislocation propagation through the material.
We will focus in more detail on two other binary systems, i.e., C-N and B-C, in which compounds with very high hardness have been predicted.

Carbon Nitrides
The main interest in studying the C-N system is a result of numerous (but unsuccessful) attempts to synthesize hypothetical ultrahard C 3 N 4 . In the 1990s, Liu et al. [30,54] and Teter and Hemley [31] predicted a number of dense low-compressibility carbon nitrides of C 3 N 4 stoichiometry that were claimed to exhibit bulk moduli and hardness higher than those of diamond because of the short length and high covalence of the C-N bonds. However, our analysis in the framework of the thermodynamic model of hardness reveals that the Vickers hardness of the densest hypothetical cubic (P−43m [30] and I−43d [31]) and pseudocubic (P−42m [31]) polymorphs of C 3 N 4 does not exceed 73 GPa.
Besides the carbon nitrides of C 3 N 4 stoichiometry that are isoelectronic with diamond [55], carbon subnitrides of C 11 N 4 stoichiometry were also studied [56] as they allow the modeling of CN x films with less than 30 at% N, which usually form when vapor phase deposition techniques are used [57]. In this context, structural models of C 11 N 4 phases accounting for the low nitrogen content were derived from diamond being isoelectronic with it, through creating defects. Indeed, diamond expressed as 2C in primitive cells has eight valence electrons, and C 11 N 4 has 11 × 4 + 4 × 5 = 64 = 8 × 8 electrons, i.e., an integer multiple of 8.
In the present paper, a novel (ultra)hard tetragonal C 11 N 4 (space group P−4m2, No. 115) was derived from C 16 , a 2 × 2 × 1 cell of body-centered tetragonal C 4 diamondlike structure [58] (Figure 3a). In this template, a defect was created at the center of the tetragonal cell by removal of the yellow carbon atom (see Figure 3a) and replacing the four surrounding carbon atoms with nitrogen (blue spheres in Figure 3b). The resulting fully relaxed carbon subnitride C 11 N 4 was analyzed for the cohesive energy E coh obtained from subtracting the atomic energies of 11 C and 4 N atoms from the calculated total energy. C 11 N 4 was found to be cohesive, with E coh = −1.93 eV/atom, which is lower than the corresponding value for pristine C 16 (E coh = −2.49 eV/atom, the value identifying diamond). This is quite expected, since C 11 N 4 results from the defect diamond-like structure of C 16 . Such observations are also valid for other binary compounds resulting from the perturbation of the diamond lattice. nitrides of C3N4 stoichiometry that were claimed to exhibit bulk moduli and hardness higher than those of diamond because of the short length and high covalence of the C-N bonds. However, our analysis in the framework of the thermodynamic model of hardness reveals that the Vickers hardness of the densest hypothetical cubic (P−43m [30] and I−43d [31]) and pseudocubic (P−42m [31]) polymorphs of C3N4 does not exceed 73 GPa.
Besides the carbon nitrides of C3N4 stoichiometry that are isoelectronic with diamond [55], carbon subnitrides of C11N4 stoichiometry were also studied [56] as they allow the modeling of CNx films with less than 30 at% N, which usually form when vapor phase deposition techniques are used [57]. In this context, structural models of C11N4 phases accounting for the low nitrogen content were derived from diamond being isoelectronic with it, through creating defects. Indeed, diamond expressed as 2C in primitive cells has eight valence electrons, and C11N4 has 11 × 4 + 4 × 5 = 64 = 8 × 8 electrons, i.e., an integer multiple of 8.
In the present paper, a novel (ultra)hard tetragonal C11N4 (space group P−4m2, No. 115) was derived from C16, a 2 × 2 × 1 cell of body-centered tetragonal C4 diamond-like structure [58] (Figure 3a). In this template, a defect was created at the center of the tetragonal cell by removal of the yellow carbon atom (see Figure 3a) and replacing the four surrounding carbon atoms with nitrogen (blue spheres in Figure 3b). The resulting fully relaxed carbon subnitride C11N4 was analyzed for the cohesive energy Ecoh obtained from subtracting the atomic energies of 11 C and 4 N atoms from the calculated total energy. C11N4 was found to be cohesive, with Ecoh = −1.93 eV/atom, which is lower than the corresponding value for pristine C16 (Ecoh = −2.49 eV/atom, the value identifying diamond). This is quite expected, since C11N4 results from the defect diamond-like structure of C16. Such observations are also valid for other binary compounds resulting from the perturbation of the diamond lattice.   found to be mechanically stable, with the whole set of elastic constants being positive, as well as being dynamically stable, with positive acoustic and optic phonon frequencies.
Calculated values of the hardness and elastic moduli of tet-C 11 N 4 are listed in Tables 1 and 2. Although the bulk and shear moduli of the novel carbon subnitride are lower than the corresponding values for high-density C 3 N 4 , its Vickers hardness H V = 76 GPa is higher than that of all hypothetical high-density polymorphs of C 3 N 4 [30,31,54].
The electronic band structure of the novel tetragonal C 11 N 4 is shown in Figure 3c, characterizing an insulator with a gap value slightly below 5 eV, from Γ (valence band) to Z (conduction band), similar to that of diamond, as a result of both structures being isoelectronic, as mentioned above.

Boron Carbides
Rhombohedral boron carbide, B 4 C (B 12 C 3 ), is the most important, well-studied and widely used compound of the B-C binary system; however, its Vickers hardness does not exceed 37 Gpa [59]. The synthesis of diamond-like BC 5 with hardness above 70 GPa [5] has stimulated interest in the study of carbon-rich compounds of this system. Thus, four polymorphs of BC 5 were predicted from first-principles structural optimizations [32], for two of which (tetragonal I−4m2 and triclinic P−1) Vickers hardness of approximately 80 Gpa was claimed. However, our assessments (the results for the densest tetragonal BC 5 are presented in Tables 1 and 2) showed that the hardness of all predicted phases was overestimated by~15% as a result of the use of unreliable empirical correlations between the shear modulus and hardness.
A similar situation is observed in the case of five predicted BC 7 polymorphs [33]: the claimed hardness values (e.g., 81 GPa for orthorhombic Pmm2 and 78 GPa for cubic P−43m phases) are also overestimated by 10-15% as a result of using the empirical microscopic hardness model (see our results for cubic BC 7 in Tables 1 and 2), and thus these phases, as well as BC 5 polymorphs, cannot be considered ultrahard. However, one may expect that a further decrease in boron content will be accompanied by a hardness increase in the formed B-C binary compound(s).
Very recently, the latter has been confirmed by the prediction of trigonal BC 11 (space group P3m1, No. 156) [34], produced by the substitution of carbon with boron in the diamond-like hex-C 12 template [60], which led to the lowering of the crystal symmetry down to trigonal.
In the context of the energy criterion, it was interesting to position the novel carbonrich BC 11 among other B-C binary compounds. Comparison of trig-BC 11 's cohesive energy with those reported for trig-BC 5 [32] and trig-BC 7 [33] shows a clear trend of decreasing E coh /atom with increasing boron content: −2.49 eV (diamond) < −2.33 eV (BC 11 : 8.3 at% B) < −2.24 eV (BC 7 : 12.5 at% B) < −2.16 eV (BC 5 : 16.7 at% B). Besides being more cohesive than the two other binary compounds, trig-BC 11 was also found to be mechanically (elastic constants) and dynamically (phonon band structures) stable.
The crystal structures of template hex-C 12 and trig-BC 11 are shown in Figure 4. In hex-C 12 , the carbon network is perfectly covalent in all dimensions. Changes are observed for trig-BC 11 , featuring a large covalent part where the carbon networks remain as in hex-C 12 , but not in the surrounding area of boron atoms, whose charge density is transferred to carbon due to the larger electronegativity of carbon. Regarding its electronic band structure, bands belonging to boron states were found crossing the Fermi level E F , signaling a metallic character arising from one electron-less B (2s 2 ,2p 1 ) versus C (2s 2 ,2p 2 ) in the wide gap insulating diamond.  The hardness of trig-BC11 calculated using four models, as well as other mechanical properties, is given in Tables 1 and 2. Although the introduction of boron atoms into the diamond crystal lattice predictably lowers the hardness, it remains high enough (> 80 GPa) to consider trig-BC11 as an ultrahard phase, in contrast to other reported binary compounds of the B-C system [32,33,61].

Ternary Compounds
Interest in the search for possible 'hybrid' structures of carbon and boron nitride and prediction of their properties (mechanical, in particular) has especially grown after the synthesis of ultrahard cubic BC2N [3]. Over the past 20 years, several dozen papers have been published on the subject, but here we will focus only on those that claim the 'ultrahardness' (HV ≥ 75 GPa) of the predicted phases. In addition to different BC2N structures The hardness of trig-BC 11 calculated using four models, as well as other mechanical properties, is given in Tables 1 and 2. Although the introduction of boron atoms into the diamond crystal lattice predictably lowers the hardness, it remains high enough (>80 GPa) to consider trig-BC 11 as an ultrahard phase, in contrast to other reported binary compounds of the B-C system [32,33,61].

Ternary Compounds
Interest in the search for possible 'hybrid' structures of carbon and boron nitride and prediction of their properties (mechanical, in particular) has especially grown after the synthesis of ultrahard cubic BC 2 N [3]. Over the past 20 years, several dozen papers have been published on the subject, but here we will focus only on those that claim the 'ultrahardness' (H V ≥ 75 GPa) of the predicted phases. In addition to different BC 2 N structures [35,36,62,63], ultrahard phases of BCN [64], BC 4 N [ [65][66][67], BC 6 N [68,69] and BC 10 N [70] compositions have been reported.
Our assessments for orthorhombic (P222 1 ) and trigonal (P3m1) BC 2 N [35], as well as for novel rhombohedral (R3m) BC 2 N [36] (Figure 5a), show that they all have Vickers hardness of the order of 75 GPa (see Tables 1 and 2; the H V values calculated using the thermodynamic model are the most reliable), i.e., almost the same as experimental value 76(4) GPa for cubic BC 2 N [3,4]. As for the 79.7 GPa hardness claimed for the calculated low-energy zinc-blended BC 2 N [63], it seems to be overestimated due to the use of the empirical hardness model suggested by Gao et al. [71] (see Table 3). A similar situation is observed for trigonal (P3m1) [65] and orthorhombic (Imm2) [66] BC4N, and tetragonal (P−42m) [68], rhombohedral (R3m) [68] and monoclinic (Pm and Cm) [69] BC6N. In all these cases, the use of Gao's model results in a 4-11% hardness overestimation compared to the values that we obtained in the framework of the thermodynamic model of hardness (see Table 3).
A similar situation is observed for trigonal (P3m1) [65] and orthorhombic (Imm2) [66] BC 4 N, and tetragonal (P−42m) [68], rhombohedral (R3m) [68] and monoclinic (Pm and Cm) [69] BC 6 N. In all these cases, the use of Gao's model results in a 4-11% hardness overestimation compared to the values that we obtained in the framework of the thermodynamic model of hardness (see Table 3).
The use of another empirical so-called 'microscopic' model of hardness suggested by Tian et al. [9] for trigonal BC 4 N [67] and BC 10 N [70] leads to even higher (15-20%) overestimations of hardness compared to the thermodynamic model (see Table 3).
Below, we discuss in more detail the tetragonal (P4 2 mc) equiatomic boron carbonitride ( Figure 5b) recently proposed by us using crystal chemistry rationale and DFT calculations [37].
Compounds containing CN − anions (such as ionic sodium cyanide Na 1+ CN 1− ) are called 'cyanides'. Since boron is a metalloid (i.e., halfway between a metal and non-metal), its combination with nitrogen leads to the equiatomic boron nitride BN that could be expressed as 'B 3+ N 3− ' considering its ionic-like nature. However, BN is rather a polar covalent compound, as can be inferred from the Pauling electronegativity difference (∆χ). In the framework of the crystal chemistry approach, we considered three template structures for BCN: octahedral (CoCN), square-planar (NaCN) and linear (CuCN). The square-planar BC2N2 coordination was found to be the most stable among all three templates in terms of cohesive energy, but despite the relative stability, it remains a non-compact 2D-like structure. Therefore, as a 3D template allowing tetrahedral coordination for B with C and N, and connected BC2N2 tetrahedra, we used tetragonal hexacarbon C 6 , so-called 'glitter' [48], which possesses two types of carbon, tetrahedral C1 and trigonal C2, the latter forming C2-C2 pairs that separate the C1C2 4 tetrahedra. The structure shown in Figure 6a featuring the charge density projections reveals the two types of carbons in the 'glitter' structure, and the corresponding C1C2 4 tetrahedra (C1: sp 3 -like carbon). With appropriate substitutions of carbon for boron and nitrogen leading to BCN, the ground state energy configuration of the derived 3D structure was found to be more cohesive than the 2D-like candidate mentioned above. pact 2D-like structure. Therefore, as a 3D template allowing tetrahedral coordination for B with C and N, and connected BC2N2 tetrahedra, we used tetragonal hexacarbon C6, socalled 'glitter' [48], which possesses two types of carbon, tetrahedral C1 and trigonal C2, the latter forming C2-C2 pairs that separate the C1C24 tetrahedra. The structure shown in Figure 6a featuring the charge density projections reveals the two types of carbons in the 'glitter' structure, and the corresponding C1C24 tetrahedra (C1: sp 3 -like carbon). With appropriate substitutions of carbon for boron and nitrogen leading to BCN, the ground state energy configuration of the derived 3D structure was found to be more cohesive than the 2D-like candidate mentioned above.
(a) (b) Figure 6. Crystal structures of template 'glitter' C6 (a) and tetragonal BCN (b) with charge density distributions. Boron, carbon and nitrogen atoms are shown by green, brown and grey spheres, respectively. The resulting BCN (B 2 C 2 N 2 ) structure sketched in Figure 6b shows BC 2 N 2 tetrahedra replacing the C1C2 4 tetrahedra in 'glitter' C 6 , and large differences in charge density distribution as compared to C 6 , with a charge density concentration skewed toward C-N bonds, and larger intensity on N versus C. Such B → C → N charge transfers are expected from the Pauling electronegativities: χ(B) = 2.04 < χ(C) = 2.55 < χ(N) = 3.44. In the case of BN, three electrons depart from B to N, leading to B 3+ N 3− . In BCN, we equally observe B 3+ , but the negative '3-' charge is now distributed between C and N, with a larger value on N due to its larger electronegativity versus C, and one obtains B 3+ C 0.316− N 2.684− . As follows from our results, cohesive tetragonal BCN is mechanically (elastic constant) and dynamically (phonon) stable. An interesting feature of this phase is that, at a relatively low density (2.783 g/cm 3 ), it is characterized by very high hardness, H V = 65 GPa (i.e., harder than single-crystal cubic boron nitride, with density of 3.486 g/cm 3 [72]), highly likely due to the presence of both tetrahedral (sp 3 ) and trigonal (sp 2 ) carbons in its crystal structure. Such mixed hybridizations in tetragonal BCN lead to its weakly metallic behavior, as illustrated by the electronic band structure in Figure 5d, exhibiting a few bands crossing the Fermi level E F .

Conclusions
The modern high-pressure synthesis of superhard and, especially, ultrahard phases is a fascinating area of research that could lead to the production of industrially important new materials. However, this field is still in its infancy, and a large number of new super-and ultrahard phases still remain to be discovered. Theoretical predictions play an important role in the present search for advanced materials with desired properties (mechanical, in particular). In this review, we have illustrated, with selected examples, the wealth of (ultra)hard allotropes and phases in the B-C-N ternary system, the theoretical (crystal chemistry considerations combined with quantum mechanics calculations) study of which is a very active area of research.
At the same time, precise calculations of the mechanical properties of superhard materials (hardness, in particular) often lie beyond the capabilities of the most advanced and modern techniques. Thus, we should state that neither widely used empirical models [7][8][9][10]71] nor machine learning [14,70,73] allow us to reliably estimate the hardness of newly predicted superhard and, especially, ultrahard phases. The only model that seems to work in these cases is the thermodynamic model [27]. Moreover, not all theoretically predicted structures exist or can be synthesized. A vivid illustration is the case of hypothetical cubic C 3 N 4 with a bulk modulus claimed to be higher than that of diamond [31]. Despite enormous efforts (new attempts are still being undertaken), this phase has not been synthesized so far, and its expected ultrahardness has never been demonstrated.
Finally, it should be noted that the search for new ultrahard phases is indeed at the frontier of fundamental science and promises great prospects for the creation of new materials that are needed for existing and prospective applications. However, the recent advances in this field clearly show that phases with hardness exceeding that of diamond are highly unlikely, or even impossible [40].
Funding: This research received no external funding.

Data Availability Statement:
The data presented in this study are available on request.

Conflicts of Interest:
The authors declare no conflict of interest.

Appendix A. Computational Framework
Accurate energy-based studies are primarily achieved within the framework of quantum mechanics. The most successful framework is density functional theory (DFT), whose theoretical basis was presented by Hohenberg and Kohn in 1964 [74]. One year later, Kohn and Sham devised the so-called KS equations to efficiently solve the wave equation of the system [75] using computational codes based on DFT. Among the numerous programs, we used the plane-wave Vienna Ab Initio Simulation Package (VASP) code [76,77] with the projector augmented wave (PAW) method [77,78] for the atomic potentials with all valence states, especially in regard to light elements such as boron, carbon and nitrogen. The exchange-correlation effects inherent to DFT were considered with the generalized gradient approximation (GGA) following Perdew, Burke and Ernzerhof (PBE) [79]. The PBE scheme was affirmed by Kilmes et al. as the best suited for quantum calculations of carbon structures [80]. A conjugate-gradient algorithm [81] was used in this computational scheme to relax the atoms onto the ground state. The tetrahedron method with Blöchl et al. corrections [82] and the Methfessel-Paxton scheme [83] were applied for both geometry relaxation and total energy calculations, respectively. Brillouin-zone integrals were approximated using a special k-point sampling scheme by Monkhorst and Pack [84]. The optimization of the structural parameters was performed until the forces on the atoms were less than 0.02 eV/Å, and all stress components below 0.003 eV/Å 3 . The calculations were converged at an energy cut-off of 500 eV for the plane-wave basis set concerning the k-point integration in the Brillouin zone, with a starting mesh of 6 × 6 × 6 up to 12 × 12 × 12 for the best convergence and relaxation to zero strains. In the post-treatment process of the ground state electronic structures, the electron localization function and the electronic and phonon band structures were computed, visualized and assessed. Calculations of phonon dispersion curves were also carried out to verify the dynamic stability of the proposed structures. The phonon modes were computed via finite displacements of the atoms off their equilibrium positions to obtain the forces from the summation over different configurations. The phonon dispersion curves along the direction of the Brillouin zone were subsequently obtained using the "phonopy" code based on the Python language [85].
Further investigations of the electronic band structure and electron density of states were carried out with the full-potential DFT-built augmented spherical-wave ASW method [86]. Figure 4c,d show the results of band structure calculations for template diamond-like hex-C 12 and the resulting trigonal BC 11 .