Novel (Super)Hard SiCN from Crystal Chemistry and First Principles

The purpose of this work is to predict the existence of novel equiatomic SiCN based on tetragonal C6 structure the elementary building unit being the 1,4 cyclohexadiene molecule comprising both tetrahedral (sp3) and trigonal (sp2) carbons. From crystal chemistry rationale the structural transformations of C6 to SiC2 and then to SiCN ternary phase were fully relaxed to the ground states using first principles DFT-based calculation. Like early proposed C6 and SiC2, SiCN was found bonding and structurally stable from the elastic properties on one hand, and dynamically stable from the phonons, on the other hand. The Vickers hardness of SiCN is higher than that of cubic silicon carbide, a conventional superabrasive, whereas hardness of tetragonal SiC2 is slightly lower. Besides the abrasive properties, the electronic band structure indicates metal-like behavior of SiCN thus suggesting the potential for heat dissipation in operating conditions, oppositely to insulating SiC.


Introduction
Whereas SiC is a rather well-studied binary compound that is widely used for various applications [1], ternary compounds of the Si-C-N system are ill-studied.
In 1965, equiatomic SiCN was reported as made by SiC crystals grown in nitrogen atmosphere [2]. The stated structure shown in Fig. 1a is a cubic fluorite-like one with face centered Si substructure, and C and N atoms occupying the fluorite F8 cube corners. Later, Chen et al. [3] identified crystalline and amorphous SiCN in the process of thin films growth by microwave CVD; but no crystal structure was reported. The conducted nano-indentation study revealed a hardness value of 30 GPa for the crystalline phase and 22 GPa for the amorphous phase, therefore, the authors' claim that SiCN is a hard material that rivals cubic BN is unjustified. In 1997, Riedel et al. [4] claimed the first ternary crystalline solid with Si 8 C 4 N 16 stoichiometry built of SiN 4 tetrahedra interconnected via carbon atoms, thus leading to N 3 SiN--C--NSiN 3 -like chains of tetrahedra (blue, Fig. 1c) along vertical b-axis, and corner sharing along horizontal planes. The latter feature is also characteristic of the novel SiCN proposed herein, with a mixed N-C coordination polyhedron for Si.
In this work, a novel equiatomic SiCN phase is proposed based on crystal chemistry rationale then quantified by ab initio calculations within the well-established density functional theory (DFT) [5].
Samir F. Matar is a former DR1-CNRS senior researcher at the University of Bordeaux, ICMCB-CNRS, France.
(iii) the tetrahedral coordination (sp 3 -like), common to Si, C, and N which also exhibit linear sp. and trigonal sp 2 types. However, in rare occurrences, Si = Si (double bonds) and Si ≡ Si (triple bonds) are known mainly in organic chemistry as in tetramesityldisilene [6].
As a case study on the interconnectedness between solid state and molecular chemistry regarding carbon, we have shown recently that the rarely occurring tricarbon C 3 molecule can be a building block of 3D covalent ultra-hard carbon networks, with carbon showing two hybridization types: linear-and tetrahedrallike [7]. Earlier, Bucknum and Hoffman in 1994 [7] conceived a dense hexacarbon tetragonal structure with linear C=C (sp 2 ) and tetrahedral C(sp 3 ) based on 1,4-cyclohexadiene. Such molecule depicted in Fig. 2a is used as an effective hydrogen donor for catalytic hydrogenation reactions. Whereas in butadiene there are alternating double and single bonds between the six carbon atoms (Kekulé formulation) in two resonant formulas, 1,4-cyclohexadiene shows a localization of two C=C double bonds. The corresponding 3D (D = dimension) structure called 'glitter' [8] i.e., shining due to the conductive character of the electronic structure as shown in the bad structure section development of this work (cf. Fig. 5). The crystal structure featuring a C6 hexahedron is shown in Fig. 2b. It consists of two types of carbons labeled C1 (tetrahedral) and trigonal C2, letting express the chemical formula as C1 2 C2 4 .
Later, a silicon dicarbide Si 2 C 4 ( Fig. 2c) was theoretically proposed based on tetragonal C 6 , by replacing tetrahedral C1 by Si leading to SiC4 tetrahedra connected through the two C2 characterizing the structure [9]. With the objective of identifying SiCN ternary compound in view of the inefficiency of the cubic setup proposed earlier [2], or of cyanide-like SiCN, we come up with the equiatomic stoichiometry through replacing one of the two trigonal carbons (C2) by N leading to SiC2N2 Fig. 1 Si-C-N ternary structures within preliminary assessments. a Cubic Fluorite-like equiatomic c-Si 4 C 4 N 4 ref. [2], and b Si 4 C 8 derived from cubic Si 4 C 4 N 4 ; c Si 8 C 4 N 16 ref. [4] showing SiN4 tetrahedra in blue color, d optimized SiCN derived from orthorhombic NaCN [15]. Blue, brown, and grey spheres correspond to Si, C and N atoms, respectively tetrahedra. The ad hoc structure submitted to unconstrained geometry relaxation led to a lowering of the tetragonal space group symmetry to P4 2 mc N°105. The ground state equiatomic structure is depicted in Fig. 2d. To the best of our knowledge, no propositions have been made of equiatomic SiCN structure or existence other than the 1965 paper [2].

Computational Framework
The search for the ground state structures and their energies such as those devised from crystal chemistry rationale, is a prerequisite to discriminate stable from unstable chemical systems prior to further investigations of the physical properties. Such a task is best carried within quantum mechanical computations based on the DFT [5]. In present work, our investigations, were carried out using the Vienna Ab initio Simulation Package (VASP) code [10,11] and the builtin projector augmented wave (PAW) method [11,12] considering atomic potentials. The DFT exchange-correlation XC effects were considered within the generalized gradient approximation (GGA) [13] found herein to perform like more sophisticated hybrid functional as the GGA-HSE06 [14]. From heavy preliminary calculations, cohesive energies with hybrid functional were found close but slightly larger than with plain GGA results. Such overestimation could lead to validate unstable structures/compositions as the SiCN found unstable also from elastic and dynamic properties in the elaboration of this work.
From the ground state structures the sets of elastic constants C ij were calculated, and the mechanical properties such as hardness and fracture toughness were estimated using contemporary theoretical models. To verify the dynamic stability of the different phases, the phonon dispersion bands were calculated. The approach consists of computing the phonon modes through finite displacements of the atoms off their equilibrium positions to subsequently obtain the forces from the summation over the different configurations. The phonon dispersion curves along the main lines of the tetragonal Brillouin zone (BZ) were obtained using first-principles phonon calculation code "Phonopy" [15]. Finally, the electronic band structures showing small gap insulating to metallic behaviors were assessed.

Energies and Crystal Structures
The structures presented in the introduction were submitted to unconstrained geometry relaxation to obtain the ground state configuration. Energy is the prevailing criterion for a first assessment of the crystal structures.
From the total electronic energies (E tot ) the cohesive energies were obtained by deducting the atomic contributions. In view of the different stoichiometries, only the atom averaged cohesive energies can be compared.

Cubic SiCN Structure Templates
The structure of cubic SiCN (space group F43 m, N°216) is shown in Fig. 1a resembling fluorite-type is built of face centered Si substructure with C and N occupying cube corners [2]. Replacing N by C to model cubic SiC 2 for the sake of comparison in present work leads to higher symmetry fluorite-type (Fm3m N°225) shown in Fig. 1b. Both cubic SiCN and SiC 2 were submitted to unconstrained geometry optimization calculations of the energy. The numerical results of the cohesive energies obtained by subtracting the atomic contributions from the total electronic energy are shown in the first lines of Table 1. They reveal small positive magnitudes leading to doubt the validity of cubic SiCN structure and to search for alternatives.
The ordered structure of orthorhombic equiatomic sodium cyanide NaCN [16] was considered by replacing Na by Si and operating full geometry optimizations. The resulting geometry consisted of a layered structure with square planar SiC2N2 motifs (Fig. 1d). Such coordination is unusual for Si. The cohesive energy, while being negative, i.e., favored with respect to cubic SiCN, is of a lower magnitude versus the other configurations based of C 6 Glitter as shown in next paragraph to feature tetrahedral coordinated silicon. Such layered structure hypothesis with square planar SiC2N2 motifs was subsequently abandoned.
Lastly, the other alkaline cyanides ACN (A = Li, K, Cs) characterized solely by A-N connections were also considered in the preliminary calculations but all led to lower cohesive energies versus square-planar NaCN-type SiCN (Fig. 1c). Such results can be chemically expected since ACN cyanides are A+ CN-, i.e.,"ionic", whereas SiCN cannot be a cyanide but rather a carbonitride, whence the interest in proposing a stable phase with tetrahedral SiC2N2 motifs.

Glitter-Derived SiCN
The second part of Table 1 shows cohesive chemical systems all larger than the preliminary results at upper part of the table. C 6 which contains both tetrahedral and trigonal carbon types is characterized with E coh. /at. = −2.037 eV, i.e. less cohesive than diamond made of purely C(sp 3 ) -cf. [17]. Si 2 C 2 N 2 is cohesive with E coh. /at. magnitude smaller but close to that of Si 2 C 4 . Crystallographically, the replacement of C1 by Si in C1 2 C2 4 ("glitter" C 6 , Fig. 2b) leads to hypothetic silicon dicarbide Si 2 C 4 (Fig. 2c) [9] which shows the same tetrahedral stacking as in C 6 . Both structures are defined in tetragonal space group P4 2 /mmc, N° 131 (cf . Tables 1 and  2). Since the two C2 belong to the same Wyckoff site (4i), the replacement of one carbon by nitrogen leads to Si 2 C 2 N 2 where Si is tetrahedrally coordinated with trigonal N and C. The initial tetragonal space group P4 2 /mmc N° 131 symmetry is lowered and subsequent crystal structure determinations after geometry relaxation led to define the ternary phase in a new space group P4 2 mc N°105 (Fig. 2c), i.e., with the loss of a mirror (m) symmetry due to the occupation of the fourfold (4i) ½, 0, z position in Si 2 C 4 ( Table 2) by two different atoms, N and C, which are now are in twofold (2c) ½, 0, z positions ( Table 2). The trigonal carbon interatomic distances increase from C 6 with d(C2-C2) = 1.34 Å to ~1.37 Å in Si 2 C 4 due to large Si atom resulting in the increase of a and c lattice constants. Si 2 C 2 N 2 structure has larger d(C-N) = 1.38 Å while d(Si-C) = 1.82 Å and d(Si-N) = 1.81 Å become close in magnitudes. Note that the tetrahedral angle ∠C-Si-N = 109.76° presents an interesting feature of sp 3 -like environment. This is a relevant result in view of the symmetry lowering to tetragonal space group P4 2 mc N°105. Consequently, SiCN is likely to exist in a structure with mixed N, C coordination expressed in regular SiC2N2 tetrahedra connected via C-N along the tetragonal c-axis direction. Figure 2e shows a 2 × 2 × 1 supercell of SiCN highlighting the SiC2N2 tetrahedra shown in the same blue color as Si.

Charge Density Projections
Further illustration of electronic ↔ crystal-structure relationship is provided by the charge density projections. Figure 3 shows the charge density yellow volumes in the different panels. In "glitter" C 6 ( Fig. 3a) the charge density is expectedly different between tetrahedral C1 with sp 3 -like 4-lobes yellow shapes versus trigonal C2 (sp 2 -like with linear-like yellow volumes) thus confirming the difference of chemical behavior between carbons belonging to the two crystal sites. Table 2 Calculated crystal structures parameters a) C 6 , space group P4 2 /mmc (N° 131). (Fig. 2b). ) Si 2 C 2 N 2 (SiCN) space group P4 2 mc (N° 105). (Fig. 2d). Upon replacing C1 by Si the charge density changes significantly with a transfer from Si to C2. C2-C2 groups show a larger charge envelop versus C2-C2 in C 6 (Fig. 3b).
Nevertheless, Si keeps a tetrahedral coordination with C, forming SiC4 tetrahedra.
The replacement of one C2 by N, leading to Si 2 C 2 N 2 equiatomic compound results in a differentiated charge envelop along the C-N pairs with the largest charge volume around N (Fig. 3c). These results are chemically sound is so far that they go along with the electronegativity trends: χ(N ) = 3.04 > χ(C) = 2.55 > χ(Si) = 1.80 which coincide with the trends of charge transfers from Si to both C and N, i.e., δ(N) = −2.12, δ(C) = −1.72, and δ(Si) = +3.84.

Elastic Constants
The elastic constants C ij are needed to derive the mechanical properties. They were determined by performing finite distortions of the lattice and deriving C ij from the strainstress relationship. Then the bulk B and shear G moduli are obtained by averaging the single-crystal elastic constants using, here, Voigt's method [18] based on a uniform strain. The calculated sets of elastic constants are given in Table 3. All values are positive, which is the first indication of mechanical stability. Further proofs were obtained from C ij combinations obeying rules related to mechanical stability: C ii (i = 1, 3, 4, 6) > 0; C 11 > C 12 , C 11 + C 33 − 2C 13 > 0; and 2C 11 + C 33 + 2C 12 + 4C 13 > 0.
The equations providing the bulk B V and shear G V moduli as function of the elastic constants C ij are as follows for the tetragonal system: C 6 characterized by shortest interatomic distances has the highest elastic constants. In Si 2 C 4 obtained by the replacement of C1 by larger Si, the resulting bulk and shear moduli lead to drastically lower C ij , as well as B tetr Voigt and G tetr Voigt . B tetr Voigt value agrees with literature [9]. The replacement of carbon by silicon possessing a larger radius and the overall resulting larger volume can explain the larger compressibility of Si 2 C 4 . Turning to the ternary B tetr Voigt = 1∕9 2C 11 + C 33 + 2C 12 + 4C 13 .
G tetr Voigt = 1∕15 2C 11 + C 12 + 2C 33 − 2C 13 + 6C 44 + 3C 66 . Fig. 3 Charge density projections. a tetragonal C 6 ; b tetragonal Si 2 C 4 (SiC 2 ); c tetragonal equiatomic Si 2 C 2 N 2 (SiCN); with highlighting of blue SiC2N2 tetrahedra. Blue, brown and grey spheres correspond to Si, C and N respectively phase Si 2 C 2 N 2 , the calculations of the two moduli B tetr Voigt and G tetr Voigt lead to significant increase, explained by the smaller d(Si-C) and d(Si-N) versus d(Si-C) in the dicarbide (Table 2). Such results should be further rationalized with the calculation of the hardness in following section.

Hardness
Vickers hardness (H V ) was predicted using two contemporary theoretical models. The thermodynamic (T) model [19] is based on thermodynamic properties and crystal structure, and the Lyakhov-Oganov (LO) approach [20] considers topology of the crystal structure, strength of covalent bonding, degree of ionicity and directionality.
The fracture toughness (K Ic ) was evaluated using Mazhnik-Oganov model [21] which is the only model that allows such predictions. The results relevant to the mechanical properties of the considered phases of the Si-C-N system are presented in Table 4. It is worth noting that hardness and bulk moduli (B 0 ) calculated in the framework of the thermodynamic model are in perfect agreement with available experimental data [3,[22][23][24][25]. Such results confirm furthermore the validity of the used model. Regarding novel Si 2 C 2 N 2 , our findings shown that the Vickers hardness higher than that of cubic silicon carbide c-SiC (Table 4), a conventional superabrasive, whereas the hardness of tetragonal Si 2 C 4 is slightly lower.

Dynamical Stabilities from the Phonons
Further stability criteria can be obtained from the phonon dispersions. Following the method presented in the 'Computational framework' section, the obtained phonon dispersion curves for C 6 , Si 2 C 4 and Si 2 C 2 N 2 are shown in Fig. 4. In each panel the bands run along the main lines of the tetragonal Brillouin zone.
Along the vertical direction the frequency ω is given in units of Terahertz (THz). Since no negative energy magnitudes are observed, the three structures can be considered as dynamically stable. There are 3 N-3 optical modes at high energy and 3 acoustic modes. The acoustic modes start from zero energy (ω = 0) at the Γ point, center of the Brillouin Zone, up to a few Terahertz. They correspond to the lattice rigid translation modes in the crystal (two transverse and one longitudinal). The remaining bands correspond to optic modes. Remarkably, in C 6 (Fig. 4a) the highest bands are at 50 THz, a magnitude higher than in diamond where ω max. ~40 THz [26] observed with Raman vibration spectroscopy. Such high frequency is assigned to the C2-C2 stretching vibrations. The phonons of Si 2 C 4 (Fig. 4b) have similar shape but the highest magnitude is now at ω max ~ 45 THz. This is concomitant with the overall smaller elastic constants especially C 33 and the resulting B and G values as well as the hardness. The phonons dispersions of the equiatomic Si 2 C 2 N 2 show different features versus the two binaries with many flat bands resulting from the unique structure of SiC 2 N 2 tetrahedra and the strong C-N interaction.

Electronic Band Structures
The electronic band structures were calculated with the Augmented Spherical Wave ASW method [27] based on the DFT. The crystal input data of Table 1 were used. The band structures are depicted in Fig. 5. The energy along the vertical direction is with respect to the Fermi level E F since all systems are metallic-like with bands crossing E F ; note that Si 2 C 4 shows vanishingly small band crossing at E F . Beside used GGA XC functional, additional calculations using hybrid XC functional were carried out. They didn't change the metallic character identified for C 6 and SiCN, Table 4 Mechanical properties of some phases of the Si-C-N system: Vickers hardness (H V ), bulk modulus (B), shear modulus (G), Young's modulus (E), Poisson's ratio (ν) and fracture toughness (K Ic ). All values are in GPa except the Poisson ratio; and K1c in MPa·m ½ a Thermodynamic model [19] b Lyakhov-Oganov model [20] c E and ν values calculated using isotropic approximation d Mazhnik-Oganov model [21] Phase (space group)  while bringing little changes to SiC 2 that could be considered as semi-conducting with a very small band gap. However, all three phases are far from the behavior of electronically insulating diamond characterized with a large band gap of 5 eV. Also, it needs to be mentioned that SiC and other non-oxide compounds with silicon are characterized with small band gaps.

Conclusions
The present work focused on the identification of a viable structure of SiCN equiatomic following a protocol involving crystal chemistry rationalization of structures submitted in a second step to unconstrained geometry optimization with DFT-based quantum mechanical calculations. Among different hypotheses, the most stable structure consists of SiC 2 N 2 tetrahedra with a stacking resembling tetragonal "glitter" C 6 and the "glitter"-derived silicon dicarbide. Confirmation of the stabilities for "glitter" C 6 , silicon dicarbide and SiCN were inferred from the cohesive energies. The mechanical (elastic properties) and dynamic (phonon dispersions) stabilities have been confirmed. From the quantitative hardness assessments SiCN was found to be (super)hard, harder than SiC. This result makes it a prospective superabrasive. Furthermore, metallic properties of SiCN revealed from the electronic band structure provide the potentiality of evacuating heat during operation (machining, etc.), oppositely to SiC which is an insulator.