I. INTRODUCTION
Superhard material is a new functional material in science research and industrial application and is a much-discussed topic. Diamond is considered as the hardest natural superhard material. However, the application of diamond has many limits because it oxidizes in high temperature and undergoes chemical reaction with metals such as iron and nickel. Therefore, it is imperative to find a new superhard material to replace diamond. For example, cubic boron nitride (c-BN) is a typical artificial superhard material with antioxidation property.Reference Kaner, Gilman and Tolbert1–Reference Solozhenko3 Because the crystal structure of diamond and c-BN is similar, it is anticipated that B–C–N compound possibly has high hardness and antioxidation under the condition of high temperature or pressure. Increasing attention is being paid on the synthesis of cubic B–C–NReference Nakano, Akaishi, Sasaki and Yamaoka4–Reference Zhao, He, Daemen, Shen, Schwarz, Zhu, Bish, Huang, Zhang and Shen8 and theoretical prediction.Reference Lambrecht and Segall9–Reference Zhu and Gao15 Solozhenko and Zhao et al. have reported the synthesis of cubic BC2N and they found that the hardness of BC2N is next to diamond.Reference Solozhenko, Dub and Novikov7,Reference Zhao, He, Daemen, Shen, Schwarz, Zhu, Bish, Huang, Zhang and Shen8 Moreover, by theoretic methods, Sun et al. have found that BC2N has a higher volume modulus and elastic shear modulus.Reference Sun, Jhi, Roundy, Cohen and Louie13
From the results of theoretical research, BC2N has seven unequivalence structures, and Sun has found that the crystal, with the P2221 space group (group number is 17), has the lowest total energy, which is possibly the stable and high density phase reported by Zhao et al. in the experiment.Reference Zhao, He, Daemen, Shen, Schwarz, Zhu, Bish, Huang, Zhang and Shen8 The unit cell of P2221-BC2N is composed of one nitrogen (N), one boron (B), and two carbon atoms (denoted by C1 and C2), which is shown in Fig. 1. The atoms of N, B, C1, and C2 form four types of covalent bonds, including C1–C2, C1–B, B–N, and C2–N. The structure is similar to the carbon tetrahedron.Reference Azevedo16 Moreover, the distribution of these four bonds in different directions is even, which would lead to high hardness and elastic modulus. Thus, BC2N has potential application in industry and science research.Reference Zhang, Sun and Chen17
Crystal structure of BC2N and its electron density difference under the ground state.
b–a, c–a data, and crystal volume as a function of temperature. (a) Linear increasing of b and c with a; (b) V–P curve.
B-M EOS fit of E–V and P–V data.
Nowadays, the research on BC2N is focused on its crystal and electron structures, and its mechanical. Nevertheless, in the field of engineering application of superhard materials, thermodynamic properties under extreme conditions are quite significant, which would greatly affect the life of the material. However, study on thermodynamic properties of BC2N at high temperature and pressure is barely reported.
In this study, based on the first principle and quasi-harmonic Debye model, thermodynamic property in extreme conditions for BC2N has been studied in detail. Our results will guide the theoretical and experimental research on B–C–N superhard system.
II. THEORETICAL METHOD
Based on first-principle calculations, total energy (E) and electron structures have been obtained using CASTEP code.Reference Payne, Teter, Allan, Arias and Joannopoulos18,Reference Milman, Winkler, White, Pickard, Payne, Akhmatskaya and Nobes19 In addition, thermodynamic properties such as thermal capacity (C v) and Debye temperature (θ) under extreme conditions have been calculated from quasi-harmonic Debye model (GIBBS package).Reference Blanco, Francisco and Luana20
For the purpose of obtaining the best exchange correlation function and testing the feasibility of first-principle calculations in the present study, exchange correlation functions including LDA-CA-PZ,Reference Perdew and Zunger21 GGA-PBE,Reference Perdew, Burke and Ernzerhof22 GGA-RPBE,Reference Hammer, Hansen and Nørskov23 GGA-PW91,Reference Perdew, Chevary, Vosko, Jackson, Pederson, Singh and Fiolhais24 GGA-WC,Reference Wu and Cohen25 GGA-PBESOLReference Perdew, Ruzsinszky, Csonka, Vydrov, Scuseria, Constantin, Zhou and Burke26 have been tested. The results are shown in Tables I and II. The test results show that LDA-CA-PZ is the optimal exchange correlation function. In the present study, ultrasoft pseudopotentialReference Vanderbilt27 has been used. The cut-off energy is 310 eV. Pseudo atomic calculations are performed for B: 2s 22p 1, C: 2s 22p 2, N: 2s 22p 3; self-consistent convergence of the total energy is 5 × 10−6 eV/atom. For the Brillouin-zone sampling, we use the 4 × 4 × 4 Monkhorst-Pack meshes.
Parameters of c-BC2N under the ground state calculated by different exchange correlation functions together with many experimental and theoretic data.
Atom charge population and bond population of c-BC2N under the ground state calculated by LDA method.
a Ref. Reference Hill32.
As for quasi-harmonic Debye model calculation, the nonequilibrium Gibbs function G*(V;P,T) of solidReference Blanco, Francisco and Luana20,Reference Ren, Wang, Zheng, Xu and Sun28 is shown in Eq. (1).
$${G^ * }(V;P,T) = E(V) + PV + {A_{{\rm{vib}}}}[\theta (V),T]\quad {\rm{,}}$$where E(V) is the total energy of the unit cell, P denotes static pressure, V and T are cell volume and temperature, respectively; θ(V) is Debye temperature; A vib is vibration Helmholtz free energy. In the quasi-harmonic Debye model and considering density of phonon's states, A vib is expressed in Eq. (2).Reference Francisco, Recio, Blanco, Pendás and Costales29–Reference Flórez, Recio, Francisco, Blanco and Pendás31
$${A_{{\rm{vib}}}}\left( {\theta ;T} \right) = nkT\left\{ {{{9\theta } \over {8T}}\; + 3{\rm{ln}}\left[ {1 - {\rm{exp}}\left( { - {\theta \over T}} \right)} \right] - D\left( {{\theta \over T}} \right)} \right\}\quad {\rm{,}}$$where n is the atom number in formula unit, k is Boltzmann constant, and D(θ/T) is the Debye integral.Reference Blanco, Francisco and Luana20 Debye temperature is shown in Eq. (3).Reference Flórez, Recio, Francisco, Blanco and Pendás31
$$\theta = {\hbar \over k}{\left( {6{\pi ^2}{V^{{1 \over 2}}}n} \right)^{{1 \over 3}}}f\left( \sigma \right)\sqrt {{{{B_{\rm{S}}}} \over M}} \quad .$$In Eq. (3), M is relative molecular mass in unit cell (in our study, it is 48.839). B S is adiabatic bulk modulus, σ is Poisson ratio, ħ is Plank constant. B S and f(σ) are expressed in Eqs. (4) and (5).Reference Francisco, Recio, Blanco, Pendás and Costales29–Reference Flórez, Recio, Francisco, Blanco and Pendás31
$${B_{\rm{S}}} \approx B\left( V \right) = V{{{\mathop{\rm d}\nolimits} E(V)} \over {{\mathop{\rm d}\nolimits} {V^2}}}\quad {\rm{,}}$$
$$f(\sigma ) = {\left\{ {3{{\left\{ {2{{\left[ {{{2\left( {1 + \sigma } \right)} \over {3\left( {1 - 2\sigma } \right)}}} \right]}^{{3 \mathord{\left/ {\vphantom {3 2}} \right. \kern-\nulldelimiterspace} 2}}} + {{\left[ {{{1\left( {1 + \sigma } \right)} \over {3\left( {1 - \sigma } \right)}}} \right]}^{{3 \mathord{\left/ {\vphantom {3 2}} \right. \kern-\nulldelimiterspace} 2}}}} \right\}}^{ - 1}}} \right\}^{{1 \mathord{\left/ {\vphantom {1 3}} \right. \kern-\nulldelimiterspace} 3}}}\quad {\rm{.}}$$According to Ref. Reference Poirier30, Poisson ratio is chosen to be 0.25. By minimization method, the function of G* with V can be obtained as shown in Eq. (6).
$${\left. {{{\partial {G^ * }(V;p,T)} \over {\partial V}}} \right|_{p,T}} = 0\quad {\rm{.}}$$By solving Eq. (6), equation of states of solid V(p,T) will be obtained. Moreover, isochoric heat capacity C v in different conditions can be calculated by Eq. (7), which is given in Ref. Reference Hill32.
$${C_{\rm{v}}} = 3nk\left[ {4D\left( {{\theta \over T}} \right) - {{3{\theta \mathord{\left/ {\vphantom {\theta T}} \right. \kern-\nulldelimiterspace} T}} \over {{\rm{exp}}({\theta \mathord{\left/ {\vphantom {\theta T}} \right. \kern-\nulldelimiterspace} T}) - 1}}} \right]\quad {\rm{.}}$$Moreover, V(T), the cell volume as a function of temperature, could be obtained from the results of GIBBS. More detailed calculation method on thermodynamic properties can be found in Ref. Reference Blanco, Francisco and Luana20.
According to Eq. (8), we could obtain the expansion coefficients. Where a, b, and c are the crystal constants along different axes; and a 0, b 0, c 0 are the crystal constants under the ground state, T 0 is 0 K. For crystals with different lattice constants (a = b ≠ c or a ≠ b ≠ c), thermal expansion coefficients along different axes will possibly not be the same. In our study, we find that BC2N is an orthorhombic crystal, and it is necessary to clarify the expansion coefficient along different axes. According to Eq. (8), the key point of expansion coefficient calculation is obtaining crystal constants at different temperatures. By geometry optimization method of first-principle calculations, a 0, b 0, and c 0 could be attained easily. Although one can execute geometry optimization calculation at different pressures and obtain the crystal constants conveniently, the crystal constants at high temperature could not be obtained by geometry optimization directly. Then, it is difficult to obtain expansion coefficients from first-principle calculations.
$${\alpha _a} = {a \over {{a_0}(T - {T_0})}},\;{\alpha _b} = {b \over {{b_0}(T - {T_0})}},\;{\alpha _c} = {c \over {{c_0}(T - {T_0})}}\quad {\rm{.}}$$In the present study, we execute a series of geometrical phase analysis for geometry optimization ranging from −70 to 250 GPa (negative pressure corresponding to expansion and positive pressure corresponding to compression), and obtain crystal constants under different conditions. Choosing a as an independent variable, b and c dependence on a can be expressed as b = f 1(a) and c = f 2(a), as shown in Fig. 2(a). In our study, α = β = γ = 90°, cell volume V = a·f 1(a)·f 2(a) is attained. f 1(a), f 2(a), and V in this study are expressed in Eqs. (9)–(11).
$$b = 1.012a - 0.004\quad {\rm{,}}$$
$$c = 1.019a - 0.0065\quad {\rm{,}}$$
$$V = 1.031228{a^3} - 0.069856{a^2} + 0.000626a\quad {\rm{.}}$$From the quasi-harmonic Debye model calculation, we attain the V = V(T) = a·f 1(a)·f 2(a), as shown in Fig. 2(b). By solving Eq. (11), a(T) with temperature can be obtained, and so b(T) and c(T). By this method, expansion coefficients of BC2N in different temperatures have been studied in detail.
III. RESULTS AND DISCUSSION
At present, it is difficult to fabricate pure BC2N by experimental methods; normally, the samples are a compound with many nanostructures.Reference Zhao, He, Daemen, Shen, Schwarz, Zhu, Bish, Huang, Zhang and Shen8,Reference Solozhenko, Andrault, Fiquet, Mezouar and Rubie33,Reference Guo, Liu, Luo, Yu, He, Tian, Sun and Wang34 Therefore, no accurate measurements for lattice constant, cell volume, bulk modulus, and so on, can be achieved. In consequence, theoretical results are significant to guide experimental work.
Here, our calculated results using the LDA-CA-PZ exchange correlation function are shown in Tables I and II. The average lattice constant state (
$\bar a$) and bulk modulus (B 0) under the ground are 3.57 Å and 396.8 GPa, respectively. Compared with the results of Sun (ab-initio, LDA),Reference Sun, Jhi, Roundy, Cohen and Louie13 GuoReference Guo, Liu, Luo, Yu, He, Tian, Sun and Wang34 (CASTEP, GGA-PBE), and FanReference Fan, Wu, Shen and Kuo35 (VASP, LDA), the results of the present study have little deviation. Moreover, the bond length and bond population values of our study are also close to the results of Guo.Reference Guo, Liu, Luo, Yu, He, Tian, Sun and Wang34 Therefore, we choose the LDA-CA-PZ exchange-correlation function for the latter calculation.
Third-order Birch–Murnaghan equation of states (third B-M EOS) plays an important role in describing crystal macroscopic properties, which are shown in Eqs. (12) and (13). By using CASTEP, geometry optimizations have been executed ranging from −70 to 250 GPa. In addition, we have obtained cell volume and total energy (E) at different pressures, followed by attaining the pressure and total energy as a function of cell volume (P–V and E–V data). Fitting P–V and E–V by third B-M EOS,Reference Francis36 we obtained volume V 0, bulk modulus B 0 in ground state, and first-order derivative of bulk modulus (B′), as shown in Fig. 3. The calculated results of B 0, B′, and V 0 are 393 GPa, 3.72, and 45.9 Å3, respectively. They are in agreement with the theoretical results reported by Sun,Reference Sun, Jhi, Roundy, Cohen and Louie13 Guo,Reference Guo, Liu, Luo, Yu, He, Tian, Sun and Wang34 Fan,Reference Fan, Wu, Shen and Kuo35 and experimental data by KnittleReference Knittle, Kaner, Jeanloz and Cohen5 (as shown in Table I).
$$E(V) = {E_0} + {{9{V_0}{B_0}} \over {16}}\left[ {{{\left( {{{\left( {{{{V_0}} \over V}} \right)}^{{2 \over 3}}} - 1} \right)}^3}{{B'}_0} + {{\left( {{{\left( {{{{V_0}} \over V}} \right)}^{{2 \over 3}}} - 1} \right)}^2}\left( {6 - 4{{\left( {{{{V_0}} \over V}} \right)}^{{2 \over 3}}}} \right)} \right]\quad {\rm{.}}$$
$$P(V) = {3 \over 2}{B_0}\left( {{{\left( {{{{V_0}} \over V}} \right)}^{{7 \over 3}}} - {{\left( {{{{V_0}} \over V}} \right)}^{{5 \over 3}}}} \right) \cdot \left( {1 + {3 \over 4}\left( {{{B'}_0} - 4} \right)\left( {{{\left( {{{{V_0}} \over V}} \right)}^{{2 \over 3}}} - 1} \right)} \right)\quad {\rm{.}}$$Figures 4 and 5 show the tendency of isochoric molar heat capacity (C v) and Debye temperature (θ) at temperatures and pressures, respectively. As shown in Fig. 4, although the temperature is lower than 300 K, C v decreases with pressure similar to a minus exponent tendency. Moreover, with the decrease in temperature, the tendency of C v with pressure becomes obvious. It implies that temperature and pressure have opposite effect on C v. From Fig. 4(a), it can be seen that C v is quite small at very low temperature, which is because of weak phonon oscillation. With the increase in temperature, C v increases gradually. In the whole temperature range, two conclusions are anticipated by the standard continuum elastic theory, shown in Ref. Reference Francis36. One is that although the temperature is lower than 300 K, C v comes from thermostimulation of phonon's acoustic mode and it is proportional to the third power of temperature (Debye T 3 law).Reference Debye37 The other is, at high temperature, C v has a tendency of Dulong–Petit limit.Reference Petit and Dulong38 However, in the medium temperature region, C v is decided by thermovibration of atoms. As shown in Fig. 4(a), although the temperature is lower than 100 K, C v obeys T 3 law at any pressure. With increasing temperature, C v increases gradually but slower obviously, but it does not obey T 3 law any more. From Fig. 4(a), it can be seen that in the medium temperature region, the increase of C v is very obvious, where the atoms' vibration plays an important role in C v of BC2N. Under the condition of 300 K and 0.1 GPa, C v = 23.15 kJ mol−1 K−1. With further increase in temperature, C v reaches the Dulong–Petit limit finally,Reference Francis36,Reference Petit and Dulong38 with the value about 74.3 kJ mol−1 K−1. In high temperature, C v reaching the Dulong–Petit limitReference Petit and Dulong38 implies that the vibration of atoms reaches the limit. If the vibration continues to increase, the crystal melts. As shown in Fig. 4(b), C v decreases with pressure at any temperature. In the low temperature region, the decrease of C v is not obvious because pressure barely has any effect on phonons. In the medium temperature region, the decrease of C v with pressure is the most obvious, which means that high pressure can control the atoms' vibration effectively. In the high temperature region, the decrease of C v becomes weak, implying that pressure has no important influence on the atoms’ vibration at high temperature.
Heat capacity of BC2N at different pressures and temperatures. (a) C v versus temperature; (b) C v versus pressure.
Debye temperature (θ) of BC2N as a function of temperature and pressure. (a) θ versus temperature; (b) θ versus pressure.
Figure 5 shows Debye temperature (θ) calculated by quasi-harmonic Debye model. Figure 5(a) shows that θ decreases with the temperature at any pressure. However, although the temperature is lower than 300 K, the variation of θ is quite weak but at high temperature, it becomes obvious. Moreover, with the increase in pressure, the variation of temperature with θ becomes weaker gradually. Figure 5(b) implies that at any temperature, θ increases with pressure rapidly. We also found that the higher the temperature is, the larger the increase of θ with pressure is.
Thermal expansion coefficient of solids α is an important parameter in the application. Using Eqs. (8)–(11), the linear expansion coefficients along a, b, and c axes (αa, αb, αc) in BC2N with P2221 structure have been calculated in the temperature range 0–3000 K, under the condition of 0–200 GPa, which have been plotted in Fig. 6. As shown in Fig. 6(a), α increases with temperature. Although the temperature is greater than 1500 K, the rise of α becomes weak especially in high pressure. Figure 6(b) shows that α decreases with pressure. When pressure is higher than 70 GPa, the decrease of α slows, which implies that pressure plays an important role in controlling temperature.
Thermal expansion of BC2N along a, b, and c axes as a function of temperature and pressure. (a) Thermal expansion versus temperature; (b) thermal expansion versus pressure.
From Fig. 6(a), it can also be seen that at a temperature lower than 500 K, αa, αb, αc along a, b, and c axes are almost the same. With further increase in temperature, expansion coefficients along a and b axes are still similar. However, αc is larger than αa and αb obviously, especially in high temperature. Thus, P2221-BC2N has an anisotropy in thermal expansion coefficients.
To reveal the intrinsic mechanism of anisotropy in thermal expansion coefficients, the increase of the four covalent bonds at high temperature has been described in detail (C1–C2, C1–B, C2–N, B–N). Figure 7 is the variation of four types of covalent bond distance with temperature. It can be seen that at a temperature lower than 500 K, the tendency of increase of these four types of covalent bonds with temperature is the same. With further increase in temperature, the increase of C2–N and B–N becomes larger than that of C2–N and B–N, implying C2–N and B–N are more sensitive to temperature in the high temperature region. This leads to an easier expansion in C2–N and B–N. The population and electron density difference analysis (as shown in Fig. 1 and Table II) show that C2–N and B–N bonds are not absolutely covalent but include the component of ionic bond partly. Especially for B–N bond, the covalent property is weaker than C1–C2 and C1–B. Compared with diamond, BC2N is not a perfect covalent bond crystal. As shown in Fig. 1, C2–N and B–N bonds are not symmetric in BC2N crystal with P2221 structure yet. Along c axis, the contribution of C2–N and B–N is larger, which leads to the obvious asymmetry of expansion coefficients. Moreover, the expansion coefficient along c axis is larger than a and b axes.
Increment ratio of length of the four types of covalent bond at different temperatures.
IV. CONCLUSION
In conclusion, thermodynamic property of BC2N with P2221 structure has been obtained from first-principle theory and quasi-harmonic Debye model. From first-principle calculations, the lattice constants a 0, b 0, c 0 and bulk modulus B 0 and B′ of BC2N are 3.553 Å, 3.601 Å, 3.557 Å, 393 GPa, and 3.72, respectively. From quasi-harmonic Debye model, isochoric heat capacity, Debye temperature, and expansion coefficient have been studied at different temperatures and pressures for BC2N. In normal temperature and pressure, C v = 23.15 kJ mol−1 K−1. Heat capacity of BC2N increases with the temperature but decreases with the pressure. It obeys the Debye T 3 law in the low temperature region, but in the high temperature range, it reaches to Dulong–Petit limit. In the medium temperature region, it is determined by the vibration of atoms. In the low temperature region, pressure has a weak effect on phonon, which leads to no obvious decrease of heat capacity. In the medium range of temperature, high pressure can control the vibration of atoms effectively. In the high temperature region, the decrease of heat capacity becomes increasingly smooth. Debye temperature decreases with temperature. In the low temperature region, it barely has any variation. With the increase in temperature, the decrease in Debye temperature becomes obvious. Debye temperature increases with the increase in pressure. In higher temperature, the increase in Debye temperature becomes steep. In BC2N, the strength of C2–N and B–N bonds is weak and not symmetric, which leads to a larger expansion along c axis. It is the origin of asymmetry of heat expansion coefficients of BC2N.
ACKNOWLEDGMENTS
The project was supported by the Open Foundation of Laboratory for Extreme Conditions Matter Properties, China (Grant No. 11zxjk05); First Funding Scheme of Young Teachers in Institution of Higher Education, Chongqing, China; Key Science and Technology Project of Ministry of Education 2012, China.
