Influences of carbon concentration on crystal structures and ideal strengths of B2CxO compounds in the B-C-O system

The search for novel superhard materials with special structures and improved thermal stability and hardness remains considerably experimental and theoretical challenges. Recent reports proposed that higher carbon content in ternary B2CxO compounds, which are isoelectronic with diamond, would lead to increased strength and hardness. This notion was derived from the calculated elastic parameters and empirical hardness formulas based on structural and electronic properties of the equilibrium structures. In present work, we introduce three potential ultra-incompressible and thermodynamically stable B2CxO (x ≥ 2) phases via a systematic particle swarm optimization algorithm structure searches. By evaluating the trends of the crystal configuration, electronic structure, and mechanical properties as a function of the C concentration, it is found that the high carbon concentration benefits the formation of the sp3 C-C covalent bonds and leads to the enhanced elastic moduli and ideal strengths in these B2CxO compounds. Studies of strain-stress behavior at large deformation, however, indicate that all these B2CxO compounds possess substantially lower ideal shear strengths than those of diamond and c-BN, suggesting that they may not be intrinsically superhard.

( + + + ) = ( ) · · · n X n Y n Z n 4 1 1 2 3 where n 1 , n 2 , n 3 · · · are the numbers of valence electrons in the atoms of the components, n is an integral positive number, and X, Y, Z · · · are the required concentrations. In contrast to the extensively studied B-C-N system, the ternary B-C-O compounds have become the subject of close attention during the last decade. Experimentally, several B-C-O materials at high pressure and high temperature have been synthesized, such as single crystals of interstitial phases B(C,O) x 19,20 .
Recently, a potential superhard ternary material B 2 CO which is isoelectronic with diamond, stabilized in two tetragonal diamond-like structures (tP4 and tI16) with claimed hardness of 50 GPa was predicted by Li et al. 21 . Moreover, it was also proposed that more ternary compounds which are isoelectronic with diamond exist in B-C-O system, such as B 2 C x O (x ≥ 2) with higher carbon content could possess higher strength and hardness than those of B 2 CO. This is based on the intuitive argument that more carbon content should bring the hybrid compounds closer to diamond in their mechanical properties because of the larger proportion of the strong sp 3 C-C bonds in the structures. It was more recently proposed that an orthorhombic B 2 C 5 O phase within Cmm2 symmetry 22 exhibits higher hardness (66.1 GPa) than that of B 2 CO. These important pioneering works give us broadened views into the essence of atomic binding in solids while inevitably advance our understanding on the structural stability and superior properties of B 2 C x O in the B-C-O system. Therefore, to extend the ultrahard structures for other family members and investigate the effects of the carbon concentration on the mechanical behaviors of B 2 C x O (x ≥ 2) compounds are of great interest and highly desirable. To address these issues, we here perform extensive structure searches to explore the potential energetically stable B 2 C x O (x ≥ 2) phases at ambient pressure using the prevalent developed Crystal structure AnaLYsis by Particle Swarm Optimization package (CALYPSO) 23,24 , unbiased by any known information. This method has been successfully applied to extensive structures which have been confirmed by independent experiments [25][26][27][28][29] . Two typical examples: the experimental crystal structures of cold-compressed graphite phase and synthesized cubic BC 3 phase which have been recently resolved 25,29 . Three novel tetragonal diamond-like structures with group symmetry I4 1 /amd, I-4m2, and P-4m2 for B 2 C 2 O, B 2 C 3 O, and B 2 C 5 O were uncovered, surprisingly, all these new structures can be derived from the previous proposed tP4-B 2 CO supercells. First principles calculations were then executed to characterize the trends of the crystal structures, electronic structures, and mechanical behaviors of these predicted phases.

Results and Discussion
Variable-cell simulations with 1-4 formula units (f.u.) in the unit cell were performed at ambient pressure, and our structure searches uncover the most stable structure for B 2 C 2 O, B 2 C 3 O, and B 2 C 5 O to be tetragonal I4 1 /amd (4 f.u./cell), I-4m2 (2 f.u./cell), and P-4m2 (1 f.u./cell) phase, respectively. All these three structures are schematically showed in Fig. 1, together with the previous proposed tP4-B 2 CO (2 f.u./cell, space group: P-4m2) 21 . As presented in Fig. 1, the structure characteristics of B 2 C x O compounds with increasing carbon content can be revealed as follows. Firstly, there exists a striking structural connections between B 2 C x O (x = 2, 3, 5) and tP4-B 2 CO, that each B 2 C x O compound can be derived from a different tP4-B 2 CO supercell with replacement of the partial O and B atoms by C atoms. For example, the predicted P-4m2-B 2 C 5 O ( Fig. 1(d)) can be constructed by a 1 × 1 × 2 tP4-B 2 CO supercell ( Fig. 1(a)) by replacing two boron atoms at 2g and one oxygen atom at 1c positions to ensure the stoichiometry of B:C:O is 2:5:1. Secondly, similar to the atomic bonding behaviors in tP4-B 2 CO, all nonequivalent atoms in each B 2 C x O compound are all tetrahedrally bonded with sp 3 bonding environment. Moreover, in Figures of 1(b)-1(d) for B 2 C 2 O, B 2 C 3 O, and B 2 C 5 O, the average number of C-C bond in their building blocks (which can be viewed as a pseudo-tetragonal tP4-B 2 CO unit cell, see the dashed cell in Fig. 1) increases from 1.6 to 4 with the increasing of carbon content. Thirdly, compare to those of in B 2 CO and B 2 C 2 O, the increasement of sp 3 C-C tetrahedral covalent bonds can be clearly revealed in B 2 C 3 O and B 2 C 5 O at elevated C concentrations, as shown in Fig. 1(c,d). Therefore, high carbon content might benefits the mechanical properties of these B 2 C x O compounds. The dynamical stability of a crystalline structure requires the eigen frequencies of its lattice vibrations be real for all wavevectors in the whole Brillouin zone. We have thus performed the phonon dispersion calculations for I4 1 /amd-B 2 C 2 O, I-4m2-B 2 C 3 O, and P-4m2-B 2 C 5 O at ambient pressure, as presented in Figure S1. Clearly, no imaginary phonon frequency was detected in the whole Brillouin zone for each predicted tetragonal phase, indicating their dynamical stabilities at ambient conditions. However, the earlier proposed Cmm2-B 2 C 5 O 22 is dynamical unstable with a wide range of imaginary frequency in the Brillouin zone according to our phonon dispersion calculations. The orthorhombic Cmm2 phase thus can be excluded from the candidates of B 2 C 5 O.
The calculated equilibrium lattice constants, average atomic volumes, bulk moduli, and pressure derivatives of bulk moduli for B 2 C x O compounds are listed in Table 1, along with previous theoretical values of tP4-B 2 CO 21 for comparisons. It can be seen that the lattice parameter a of each B 2 C x O (x = 2, 3, 5) is close to that of tP4-B 2 CO, and the lattice parameter c of each B 2 C x O phase is about five, three, and two times that of tP4-B 2 CO, respectively. This is in agreement with crystal configuration of each B 2 C x O pahse as presented in Fig. 1. The bulk moduli and their pressure derivatives are obtained by fitting pressures and volumes with the third-order Birch-Murnaghan equation of state (EOS) 30 . More detailed structural information on the interatomic distances of these predicted phases are presented in Scientific RepoRts | 5:15481 | DOi: 10.1038/srep15481 Supplementary Table S1. From Table 1, one can see that the bulk moduli and densities of B 2 C x O compounds increase continuously with high carbon content, agree well with the decreasement of interatomic distances (see Table S1) in B 2 C x O compounds. Note that the sp 3 C-C tetrahedral covalent bond lengths of B 2 C 3 O (1.585 Å) and B 2 C 5 O (1.575 and 1.568 Å) are slightly longer than that of diamond (1.535 Å), indicating they might possess excellent mechanical properties. Furthermore, it is important to explore the thermodynamic stability for further experimental synthesis. The thermodynamic stability for each B 2 C x O compound, with respect to the separate phases as a function of pressure, is quantified in terms of the formation enthalpies in two possible routes: where the α-B 31 , diamond C, α-O 2

32
, and I4 1 /amd-B 2 O 21 are chosen as the reference phases. As shown in Figure S2, the calculated formation enthalpies obtained by this two reaction routes are all negative, which indicate that the formations of these three phases are exothermic at ambient as well as at high pressures. Thus these B 2 C x O compounds are stable against the decomposition into the mixture of B + C + O or B 2 O + C and the syntheses of these structures are highly desirable at very readily attainable pressures. Figure 2 and 3 show the band structures, total density of states (DOS), and partial DOS for these B 2 C x O compounds. As shown in Fig. 2(a-d), the indirect semiconductor nature of B 2 CO and B 2 C 5 O are characterized by energy gaps of 1.96 and 2.42 eV between M and Γ points, while B 2 C 2 O and B 2 C 3 O are direct semiconductors with band gaps of 2.14 and 2.40 eV at Γ points in Fig. 2 (b,c). The increasement of band gaps of these B 2 C x O compounds can be clearly disclosed with increasing carbon contents. From inspection of their partial DOS curves in Fig. 3, it shows that the peaks from − 24 to − 21 eV have predominantly O-2s character. From Fig. 3(a-d), the weight of the C-2s and C-2p states increases gradually and the increased overlap of C-2s and C-2p in energy of − 19 to − 7 eV (B-2p and C-2p in energy of − 7 eV to 0 eV ) leads to an enhanced hybridization interactions between C-C atoms (B-C atoms). To analyze the chemical bonding character of different atoms in these B 2 C x O compounds, electron localization function (ELF) 33 and Bader charge analyses 34 are calculated as presented in Figure S3 and listed in Table 2. For the selected 3D ELF distributions (ELF = 0.9) of B 2 C 3 O and B 2 C 5 O, the high electron localization can be seen in the regions of C-C bond and B-C bond, which is indicative of strong covalent bonding. However, at the same ELF value for B-O bond, the ELF attains local maximum at the O sites compared to B sites, reflecting the ionicity of B-O bond. From Bader charge analysis results listed in Table 2, the strong covalent nature of the C-C and B-C bonds in B 2 C x O were quantitatively revealed by the evidences of large charge densities at their bond critical points with negative Laplacian values. Especially for B 2 C 5 O, the charge density at C1-C2 bond critical point is 1.549 electrons/Å 3 with a Laplacian value of − 11.068, which is close to that of C-C bond in diamond 35   For the tetragonal B 2 C x O phase, six independent elastic constants C ij were determined at GGA level 36 from the stress of the strained structure with a small finite strain. The elastic stability, incompressibility, and rigidity of each tetragonal B 2 C x O compound are thus studied based on the calculated elastic constants and derived Hill elastic moduli 37 , as listed in Table 3. Firstly, one can see that the mechanical stabilities of three B 2 C x O compounds satisfy the Born-Huang criterion for a tetragonal crystal 38 indicating that they are mechanically stable at ambient pressure. For B 2 C 5 O, we find its unusually high  incompressibility along a-direction, as demonstrated by the extremely large C 11 value (889 GPa) which is slightly larger than experimental data of c-BN (C 11 : 820 GPa) 39 and is comparable to that of known superhard diamond (C 11 : 1076 GPa) 40 . Secondly, our calculated elastic parameters for tP4-B 2 CO are in a good accordant with the previous theoretical results 21 using the same exchange-correlation functionals: local density approximation (LDA) 41 . Moreover, the calculated Hill bulk moduli agree well with those directly obtained from the fitting of the Birch-Murnaghan EOS (see Table 1), which further demonstrates the accuracy of our elastic constants calculations for three B 2 C x O compounds. Thirdly, a parallel increasement of elastic moduli and a reduction of Poisson's ratios can be clearly revealed with increasing carbon contents, as shown in Table 3. Since the hardness is deduced from the size of the indentation after deformation, a hard material typically requires a high bulk modulus to support the volume decrease created by the applied pressure, and a low Poisson's ratio (ν) or high shear modulus (so that the material will not deform in a direction different from the applied load) 42 . Among these ternary compounds, the calculated bulk modulus for B 2 C 3 O and B 2 C 5 O is 322 GPa and 345 GPa, respectively, indicating their highly incompressible nature. According to Teter 43 , the shear modulus is a significantly better qualitative predictor of hardness than the bulk modulus, governing the indentation hardness. The shear moduli of B 2 C 3 O and B 2 C 5 O are 302 and 351 GPa, and they are expected to withstand shear strain to a large extent. In view of the large bulk and shear moduli of these B 2 C x O compounds, the hardness calculations are of great interest. By using three different hardness models proposed by Gao et al. 44 , Chen et al. 45 , and Šimůnek et al. 46 , the theoretical hardness of these B 2 C x O phases were estimated and listed in Table 4. One can see that our calculated theoretical hardness for tP4-B 2 CO is in excellent with previous result performed by Li et al. 21 using the Gao's model. For other new predicted B 2 C 2 O, B 2 C 3 O, and B 2 C 5 O phase, the calculated theoretical hardness is in the range of 43~57GPa, 47~62GPa, and 51~68 GPa, respectively, suggesting their potential superhard nature.
The elastic anisotropy of crystal can exert great effects on the properties of physical mechanism, such as anisotropic plastic deformation, crack behavior, and elastic instability. Therefore, the elastic anisotropies of tetragonal B 2 C x O compounds were systematically studied for their further engineering applications. We here calculated the orientation dependences of the Young's modulus E and shear modulus G. For each tetragonal B 2 C x O phase, the Young's modulus E is described by the following equation 47  where α, β, and γ are the direction cosines which determine the angles between the a-, b-, and c-axis of a crystal and a given direction [uvw]. And s 11 , s 33 , s 44 , s 66 , s 12 , and s 13 are the elastic compliance constants    Fig. 4, the distance from the origin of system of coordinate to this surface equals to the Young's modulus in a given direction. For a perfectly isotropic medium this three-dimensional surface should be a sphere, however, all these B 2 C x O compounds exhibit a well-pronounced anisotropy in Fig. 4. In more detail, Fig. 5  . As plotted in Fig. 6, the orientation dependence ofthe shear modulus G of each B 2 C x O compound was also calculated  49,50 to understand the structural deformation, strength, and hardness have been extensively applied to strong solids under specified loading strains [51][52][53][54] . It is the ideal strength of a material which is defined as the stress at which a perfect crystal becomes mechanically unstable, that sets an upper bound for material strength. Studies of the strain-stress relations and the underlying atomistic bond-breaking processes can provide critical insights into the atomistic mechanism for the structural deformation and failure modes. Figure 7 presents the calculated strain-stress relations for these four B 2 C x O compounds under tensile strains in four principal symmetry crystallographic directions. It can be seen that these B 2 C x O compounds have strong stress responses in the < 100 > , < 001 > , and < 110 > directions with peak tensile stresses above 50 GPa. Especially in < 110 > directions, the peak tensile stresses of these compounds (tP4-B 2 CO: 137.  consistent with their body-diagonal alignment of the weak B-O bonds. One can see that an additional 37.3 wt% of carbon content in B 2 C 5 O over that in B 2 CO leads to 75.7% of enhancement to its ideal tensile strength. However, even with this large increase the ideal tensile strength of B 2 C 5 O is still below those of diamond and c-BN. The weakest peak tensile stress occurs in the < 111 > directions, which indicates that under tensile loadings, these B 2 C x O compounds would first cleave in the (111) plane. The critical shear stresses in these compounds are then calculated by applying < > 110 , < > 112 and < > 112 shear deformations in the (111) easy cleavage plane perpendicular to the weakest tensile direction. As shown in Fig. 8, for each B 2 C x O compound along (111) < > 110 shear direction, the lowest ideal shear strength is determined to be 2.6 GPa, 11.1 GPa, 11.4 GPa, and 16.3 GPa, which is lower than that of the lowest tensile strength, respectively. This means that the failure mode in each B 2 C x O phase is dominated by the shear type. Since indentation hardness measurements produce volume and shape changes beneath the indenter, it is expected that shear strength is more closely related to the intrinsic hardness of a material. In fact, it has been shown that the limit of the structural stability is correlated with the maximum shear strength 42 . The present ideal shear strength calculations indicate that these B 2 C x O members are likely not to be superhard materials despite their increased C content. Take B 2 C 5 O for example, we next explore the atomistic structural deformation mode along (111) < > 110 shear direction to understand this result. Figure 9 shows the lengths of the B-O bonds in B 2 C 5 O that connect the shear planes as a function of strain, and the structural snapshots of the unit cell at the critical steps near the bond-breaking points are also plotted. In Fig. 9, the equivalent B-O lengths (1.637 Å at equilibrium state) are stretched and split into four nonequivalence B-O lengths denoted as d1, d2, d3, and d4 under shear loadings. The bond lengths of d2 and d3 decrease nearly synchronously at each strain. On the contrary, d1 and d4 increase conformably with increasing strains and d1 break at the critical shear strain of 0.056, and d4 then decreases abruptly along with the breaking of d1 bonds. Such a bond-breaking can also be seen from the selected crystal structures (at shear strains of γ 1, γ 2, γ 3, and γ 4) before and after shear instability. Therefore, the origin of the lattice instability of B 2 C 5 O under large shear strain at the atomic level during shear deformation can be attributed to the breaking of weak ionic B-O bonds, and this is also the case for other B 2 C x O members. According to the results discussed above, the increased number of stronger sp 3 C-C bonds indeed contribute to larger elastic moduli, hardness, and idea strengths to these ternary B 2 C x O compounds isoelectronic with diamond. However, because of their significantly lower ideal shear strengths originated from weak B-O bonds than those of superhard diamond and c-BN, all these B 2 C x O compounds may not be intrinsically superhard.

Conclusion
In summary, we have systematically explored the crystal structures and properties of ternary B 2 C x O (x ≥ 2) compounds at ambient conditions by using unbiased structure searching techniques in combination with first-principles calculations. Three novel tetragonal diamond-like structure with group symmetry I4 1 /amd, I-4m2, and P-4m2 for B 2 C 2 O, B 2 C 3 O, and B 2 C 5 O is uncovered, respectively. They are all dynamically stable and can be synthesized at ambient conditions according to the phonon dispersions and formation enthalpies calculations. The elastic anisotropy of each B 2 C x O phase has been demonstrated by the distributions of Young's and shear moduli along different crystal orientations. A good relation between the mechanical behavior and carbon concentration of B 2 C x O was established by detailed evaluating the variations of the crystal configurations, electronic structures, and mechanical properties. The present results suggest that the high carbon content benefits that the formation of the covalent sp 3 C-C polyhedral stacking structure and contributes to larger elastic moduli and hardness. The ideal strength calculations, however, indicate that all these B 2 C x O compounds may not be intrinsically

Methods
The crystal structure searches were performed based on a global minimization of energy surfaces merging ab initio total-energy calculations as implemented in CALYPSO code 23,24 . CALYPSO code was designed to predict stable or metastable crystal structures requiring only chemical compositions of a given compound at given external conditions (e.g., pressure). Here, using the CALYPSO code in combined with Vienna ab initio simulation package (VASP) 55 , variable cell structure searches for each B 2 C x O (x = 2, 3, 5) compound containing 1-4 f.u. in the simulation cell were systematically performed at ambient pressure. During the structure searchs, the 60% structures of each generation with lower enthalpies were selected to generate the structures for the next generation by Particle Swarm Optimization (PSO) operation, and the other structures in new generation were randomly generated to increase the structural diversity. The underlying structure relaxations and electronic calculations were carried out using density functional theory with all-electron projector-augmented wave (PAW) method 56 to describe the electron-ion interactions, as implemented in the VASP code. The exchange-correlation functional was treated by the generalized gradient approximation (GGA) 36 using functional of Perdew Burke Ernzerhof (PBE) 57 . The electronic wave functions were expanded in a plane-wave basis set with cutoff energy of 600 eV, and the Brillouin zone integration was employed using Monkhorst-Pack k point meshes 58 with a grid of 0.03 Å −1 for all cases to ensure the total energies converged to be better than ~1 meV/atom. The phonon frequencies were calculated by the direct supercell approach, which uses the forces obtained by the Hellmann-Feynaman theorem calculated from the optimized supercell 59 . The strain-stress method 60 was used in calculating the single elastic constants. A set of given strains with a finite variation between − 0.01 and + 0.01 were applied on the optimized structure and the atomic position was fully optimized. Then, the elastic constants were obtained from the stress of the strained structure. The polycrystalline bulk modulus and shear modulus were thus derived from the Voigt-Reuss-Hill averaging scheme 61 . The quasistatic ideal strength is calculated by incrementally deforming the modeled cell in the direction of the applied strain and controlling the specific strain components, and simultaneously relaxing both the other strain components, as well as the atoms inside the unit cell, at each step.