Exploring configurations and properties of boron carbide by first principle

Based on the first principle, the formation energy, phonon vibration, physical property of three common B4C models were extensively study. Through the calculation of thermodynamic formation energy, it is confirmed po model has the most stable energy configuration. Combined the simulated x-ray and experimental data, it is found that the experimental boron carbide is actually composed of a variety of configurations, the majority of which is po model. Via the analysis of phonon vibration, the highest phonon frequencies of the different configurations were identified as the result of stretching vibrations from the triatomic chain. The research of electrical properties of three B4C models clarify B4C is a semiconductor but will transform to conductor at specific high pressure. The calculation of the mechanical property states that B4C is hard material while the hardness will gradually decrease with pressure increases. Both the relationship of their electrical properties and mechanical properties with pressure illustrate that the po model has the fastest structural change and ch model has the slowest change.


Introduction
Boron carbide (B 4 C) has been a class of materials considered for scientific interest and engineering applications due to the physical and chemical properties [1]. B 4 C has an extensive nuclear application as shielding material, neutron radiation absorbent, neutron detector and control rod in nuclear reactors due to its high neutron absorbing cross-section with self-heal for radiation damage [2][3][4][5][6]. As a p-type semiconductor, B 4 C is identified as a potential candidate material for electronic devices with high operating temperatures [7]. Owing to the Seebeck coefficient as high as 300 μV K −1 , B 4 C is an excellent thermoelectric material [1,2]. Due to its low density (2.52 g cm −3 ) and fantastic crack resistance on ballistic performance [1,8], B 4 C has a military application as bulletproof armor. Also the outstanding properties such as high hardness [1,[9][10][11], extreme abrasion resistance [1], high melting point (2450°C) and thermal stability [12] make B 4 C play an important role in mechanical processing industry as wear-resistant and cutting tool material. B 4 C has a rhombohedral structure (figure 1), which has been determined by x-ray diffraction [8]. Firstly, its structure has been recognized as B12 icosahedra connected by C-C-C chains, namely chain model B 4 C (ch B 4 C), Boron atoms occupy two nonequivalent sets of 6h positions while the carbon atoms in lb and 2c positions form a triatomic chain [8]. The triatomic chain is straight with the angel ψ is 180°. However a linear C-C-C chain with only two covalent bonds on the central C atom, which is difficult to understood through quantum mechanical. And the C-C bond length of 1.435 Å, obviously short than the C-C bond in diamond with length of 1.54 Å. The 60% 1b position was suggested to be occupied for B atoms after the nuclear magnetic resonance (NMR) research of B 4 C [13]. In 1971, NMR study of B 4 C indicates boron atoms occupy the central positions on the triatomic chains and carbon atoms occupy positions on the icosahedra [14]. As two inequivalent sites (polar and equatorial) exist in the icosahedron, B 4 C (B11C icosahedra and C-B-C chains) crystal structural has been evolved as: po B 4 C (carbon atoms take up polar sites forming bonds with the neighboring icosahedra) and eq B 4 C (carbon atoms occupy equatorial sites forming bonds with the C-B-C chains) [15,16]. The icosahedral Any further distribution of this work must maintain attribution to the author(s) and the title of the work, journal citation and DOI. configuration is caused by a tendency for three-centered covalent bonds due to a deficiency of valence electrons [17]. In these rhombohedral structure B 4 C (figure 1(a)), atoms on equatorial sites bonded with atoms in terminal of triatomic chain and atoms on polar sites bonded to neighboring icosahedra. In these structures, all middle atoms in the chain bond two coordination atoms and others in the chain form bonds with four coordination atoms; all other atoms are located on the icosahedra and formed bonds with six coordination atoms.
The ch B 4 C has the largest C-C bond number, eq B 4 C second, and po B 4 C contains no C-C bonds. Among those B 4 C polymorphs, po B 4 C is the most stable structure on thermodynamics at ambient pressure [16,[18][19][20]. The relative enthalpy curves for ch, and po B 4 C do not cross [20], and does not support the idea that po B 4 C can transform into ch B 4 C under pressure [21]. Experimental research reveals B 4 C is consisted of by overwhelming majority po B 4 C, which has also been verified by theoretical study [15,16,18,19,22,23].
In this paper, we executed an extensive research about B 4 C: Firstly, the energy stability of three B 4 C models is revealed through energy analysis, and then the simulated x-ray was compared with experiment data to reveal which configurations are dominant. Next, the phonon vibration modes and their frequency are analyzed to illustrate the fundamental source of phonon vibration difference. Also the structural deformation of B 4 C with pressure was simulated to explain their deformation mechanism under hydrostatic pressure. And the study of electrical and mechanical properties was employed to give a detailed instruction of these B 4 C models.

Computational method
The structural relaxations were implemented in CASTEP module [24] by using local density approximation (LDA) in the frame of density functional theory (DFT). The local exchange-correlation functional in CASTEP is the Perdew and Zunger parameterization of the numerical results of Ceperley and Alde [25]. To give well converged total energies of 0.01eV, plane-wave basis set with an energy cutoff of 750 eV was adopted and k-points separation (2π×0.04 Å −1 ) was assigned to generate a k-points grid using the Monkhorst-Pack grid parameters. To ensure that the structures are dynamically stable, the phonon frequencies were calculated throughout the Brillouin zone using linear response method [26]. Considering LDA generally underestimates the energy gap by 30%∼40% relative to experimental value, when calculating the electrical properties, the Heyd −Scuseria−Ernzerhof hybrid functional (HSE06) was also adopted [27], such method can obtain an accurate prediction of band gaps. The elastic matrix was obtained based on an efficient stress-strain method, with 9 distorted structures generated for each strain model and maximum strain 0.003.

Results and discussion
3.1. Thermodynamic analysis After structure optimization, the densities (g cm −3 ) of po B 4 C (2.633) and eq B 4 C (2.623) are 4.5% and 4.1% higher than the experiment density value 2.52, the slight difference may be caused by LDA underestimating the lattice parameters [28]. While the ch B 4 C has a density of 2.578, quite approach to experimental value.
We also explore the thermodynamic stability of three B 4 C models, which with respect to several separate phases as a function of pressure and can be quantified in terms of the formation enthalpies (DH).
Here α-B was chosen as B source, graphite ( ) DH g and diamond (DH d ) as two possible C sources. Considering that all elements in B 4 C are highly light, it would be important to derive their enthalpies of these phases by considering the zero-point vibrational energy. The zero-point energy can be evaluated as ( ) is the phonon density of states. DH of ch, eq, and po B 4 C (with zero-point vibrational energy or without) was calculated and presented in figure 2.
The calculation of the formation energy revealed that whether or not consider ZPE, all three models are indeed thermodynamic stability. Whatever the reactants chosen, po model has the lowest formation energy, ch model has the highest one (with about 0.40 eV per B 4 C formula higher than that of po) and eq has about 0.18 eV higher that of po, which give the thermodynamic optimality of po model than that of eq and ch models. Whatever the energy sequence or their energy differences, both are in agreement with the previous results [18,29].

X-ray analysis
We simulated x-ray diffraction data of three crystal structures of B 4 C under ambient pressure, and compared them with experimental data [9] as exhibited in figure 3. The x-ray radiation was adopted the copper as target with wavelength of 1.540 Å. During the 2-Theta range 15°to 45°, all diffraction peaks are mainly concentrated in two areas: 17.5°∼ 25°and 30°∼ 40°. Combined with the analysis of the peaks' position and their strength ratio, three B 4 C models have similar diffraction, while po model has the best match with the experimental data than others.

Phonon vibration
The phonon dispersion curves are investigated to analyze the dynamical stability of B 4 C (figure 4). Throughout the whole Brillouin zone, no imaginary frequencies were observed at ambient pressure (figures 4(a)∼(c)), which confirm three B 4 C models' dynamical stability. There exist very significant gaps separate the lattice waves of phonon vibration to two parts: high frequency part (one lattice wave with max frequency) and low frequency part (the others). It is apparent that three B 4 C models have similar frequency range of the low frequency part, while their high frequency parts are quite different. Through vibrational analysis, the max frequency of phonon vibration of ch model is 1722 cm −1 originating from the axial stretching vibration of triatomic chain C-C-C. Both po and eq models have similar the max frequency (1606 cm −1 and 1575 cm −1 , obviously smaller than that of ch model) of phonon vibration originating from C-B-C atomic chains, which due to the strength of the atomic bonds between the chains. According to the energy analysis, po model has the most energy advantages among the three models, and boron carbide is mainly composed of po model [16]. The maximum vibration frequency of po model we calculated is very close to what had been calculated previouslyby Lazzari as 1605 cm −1 [16], which also agreement with the previous results [30].

Structure evolution
There are many theoretical work about structure transformation for B 4 C under non-hydrostatic stresses [31][32][33], hydrostatic stresses [31,32], experimentally as shock-wave experiment [34,35], DAC [36][37][38], et al These works have revealed that uniaxial pressure induced irreversible bending of C-B-C chains [31-33, 36, 38], the central boron atom in C-B-C chain form bonds with boron atoms in the nearby icosahedron [32] and amorphization emerge with structural collapse [31] and carbon atoms rearrangement [36]. Under extreme condition, the central boron atom in the C-B-C chain form bonds with boron atoms in the nearby icosahedron accompany the chain disappear and new 3-atoms triangles forms.
The structural evolution of three B 4 C polymorphisms in the process of pressurization during pressure range 0∼300 GPa are simulated and exhibited in figure 5. Generally speaking, the increase of pressure will lead to a denser atomic accumulation, hence its density increases. In figure 5, the relationship between density and pressure is simulated, and in good agreement with the universal law. For the most common model, the relationship between density of po model and pressure is quite well consistent with the result of Korotaev [39]. Meanwhile, the result of theoretical calculation also matches well with the experiment [37]. It also can be found that the density changes of those three models are different: the densities of po and eq are very similar during the entire pressure range; however ch has the smallest density and the difference in density between ch and po increases with increasing pressure.
The structure of a material will deform under pressure, however the degree of deformation depends on its stiffness. During the process of pressure, all C-C bonds get short (from 1.321 Å at 0 GPa to 1.226 Å to 300 GPa), while the angle ψ in atomic chain of ch model remains 180°unchanged. As for eq and po, their angles ψ in atomic chain both already smaller than that in ch model. Once loading high pressure, the triatomic chain bends, and the angles ψ decreases, in according with the high pressure experimental result as triatomic chain continuous bends to 70 GPa [37,39]. The po model has a bigger angle change than eq model; ψ change from 178.453°at 0 GPa to 121.154°at 300 GPa (decrease about 32.11%) for po model and 179.086°at 0 GPa to 129.406°at 300 GPa (decrease about 27.74%) for eq model. And also because of the bending of triatomic chain, for po model, the center B atom moves toward and then connected the nearby two icosahedrons with B-B bonds length 1.568 Å, and eq model has a similar deformation with B-B bonds length 1.586 Å. The change of angle ψ with pressure reveals that po model has the most deformable structure with structure changes to lower symmetry as amorphous among three B 4 C models, which has been observed experimentally [31,38]. However the ch model is stiffness and non-deformation.

Mechanical properties
Elastic matrix plays an important role in investigation of mechanical properties of materials; the calculation of mechanical properties is always based on the elastic matrix. For instance, based on the elastic matrix, all independent elastic constants (C ij ) can be obtained, bulk modulus (B) and shear modulus (G) can also be acquired based on C ij [40], and then Young's modulus (E), Poisson's ratio (ν) and hardness (Hv) can be further obtained based on modulus as B and G [40,41]. First, the calculated simplest form of elastic matrixes for ch, eq and po models B 4 C are displayed as follow: The numerical unit is unified with GPa. Three elastic matrix are all symmetric about the main diagonal, and these symmetric nonzero data are represented by 'KK', the blank cell represents zero.
Based on the independent elastic constants, both bulk and shear modulus (B and G) of three B 4 C models have been calculated and listed in table 1. The modulus B and G of po and eq models agree well with experimental data. As a measure of the stiffness of a solid material, Young's modulus (E) can be obtained from B and G as E=9BG/(3B+G). Compared to the experimental value, ch model has an obviously smaller E, and po and eq models have larger E, which is consistent with the trend of their respective modulus (B, G) compared with experimental value. And the ductility/brittleness of material can be qualitative by Poisson's ratio (ν, as (3B −2G)/2(3B+G)) with the threshold value of 1/3. Three B 4 C models are all brittle because of their similar values of v are under 1/3, which are in good agreement with experimental data.
As an important parameter of mechanical properties of materials, hardness is often calculated based on the revised empirical scheme [41].
As an important physical parameter, pressure plays an important role in the process of affecting material properties. Here, pressure with range of 0 GPa to 100 GPa is considered, and the mechanical properties including (bulk moduli B, shear moduli G, Young's moduli E, Vickers hardness Hv, Poisson's ratio ν) of three B 4 C models are calculated on their primitive cells and exhibited in figure 6.
From the relationship between mechanical modulus and pressure, it is obvious that all three models have similar B during entire pressure range and they all monotonously increase with pressure, which indicated that all their compression resistance are improved. This is due to the increase of pressure will inevitably leads to the densification of material. As for the shear moduli and Young's moduli, three models have the distinct different trends: ch has a trend as increase at first and then decrease of G and E with pressure, both po and eq have increasing trend of G and E with pressure. The differential change of their E(G) with pressure may cause by their structural transformation with pressure increase. At the same time, as pressure increases, the hardness decrease and Poisson's ratio increase, which means three B 4 C models, will be more ductility at high pressure than at ambient pressure. 3.6. Electrical properties To explore the electrical properties of these B 4 C models, the electronic band structures at ambient pressure are calculated based on DFT-LDA and the band structures through all the high symmetry point in the entire Brillouin region and exhibited in figure 7. It is clear that all three models are semiconductors, which agree well with experimental results [44]. And the indirect band gaps of ch, eq and po models are 1.496 eV, 2.640 eV and 3.083 eV, respectively. However, the exchange-correlation functional LDA tends to give underestimate band gaps of semiconductors [45,46]. Hence, the Heyd-Scuseria-Ernzerhof hybrid functional (HSE06) [27] was employed, which can accurately calculate the electronic band structures of material but with huge computing costs. The electronic band structures calculated based on HSE06 are displayed in figure 8. It is obvious that all B 4 C model are semiconductors with indirect band gap significantly larger than the values based on LDA. For ch, eq and po models, the gaps are 2.947 eV, 3.557 eV and 4.173 eV, respectively. As po B 4 C has insulating feature with broad indirect gaps of 4.173 eV, agreement well with 4.13 eV reported by Ektarawong et al [29], and the result of ch model is completely consistent with the previous results as 2.97 eV [47].  Pressure significantly affects the electrical properties of materials. In general, as the pressure increases, the relative position between atoms decreases, leading to the electron overlap increase. At that moment, electrons will not be limited to one single atom or bond, hence the so-called delocalized electrons form, indicating that the electronic band gap will be reduced or even metallized. According to the classical energy band theory, by loading high pressure, the lattice parameters become smaller, and the Brillouin zone becomes larger, and the energy bands broaden, resulting in the reduction of the band gap. A typical case is the metal hydrogen [48].
Through the calculation of their band gap values (based on LDA) as a function of pressure (displayed in figure 9), we found that with the pressure increase, all band gaps of three models decrease until the band gaps disappear and the system show electrical conductivity. Among these three models, po model will firstly decreases to zero at pressure about 380 GPa and then eq and ch models decrease to zero almost simultaneously (about 500 GPa). In consideration of po model has the largest gap and ch model has the smallest gap at ambient pressure, hence po model has the fastest gap descent rate and ch model has the slowest one, which also revealed that among three models, the structure of po model is most easily deformed, and ch model is the most stiffness one.

Conclusion
By the calculation of the formation energy, the thermodynamic stability of three B 4 C models was revealed that po model is the most stable energy structure. Through the comparison of simulated x-ray and experimental  data, all three models have very similar diffraction peaks to the experiment peaks, especially po model, which declare that experimental B 4 C exists in a variety of configurations, the majority one is po model. Through the analysis of phonon vibration, the difference in vibration were mainly from the stretching vibration of triatomic chain, which give the max vibrational frequency of ch model (1722 cm −1 ) is larger than that of eq (1575 cm −1 ) and po (1606 cm −1 ) models. Under isostatic pressing, both po and eq models will undergo structural deformation with triatomic chain bending and the center B atoms coming near two neighboring icosahedrons. However, there is no similar deformation except for the overall compression density for ch model.
The research of their electrical properties clarifies the semiconductor properties of B 4 C and their band gaps gradually decrease and disappear with electrical conductivity appearing at specific pressure during loading pressure. More accurate hybrid function is used to study the precise electrical properties of B 4 C, and manifest po model has the largest band gap of 4.173 eV, agreement well with other reported values. The calculation of their mechanical property state that all B 4 C models are hard with high hardness, especially, po and eq models are in good agreement with the measured data. The study of mechanical properties with pressure also reveals the hardness of B 4 C will gradually decreases and B 4 C becomes ductility to a certain extent. Both the relationship of their electrical properties and mechanical properties with pressure illustrate that po model has the fastest structural change and ch model has the slowest change.