A New Superhard Phase and Physical Properties of ZrB₃ from First-Principles Calculations

Abstract

Using the first-principles particle swarm optimization algorithm for crystal structural prediction, we have predicted a novel monoclinic C2/m structure for ZrB₃, which is more energetically favorable than the previously proposed FeB₃-, TcP₃-, MoB₃-, WB₃-, and OsB₃-type structures in the considered pressure range. The new phase is mechanically and dynamically stable, as confirmed by the calculations of its elastic constants and phonon dispersion curve. The calculated large shear modulus (227 GPa) and high hardness (42.2 GPa) show that ZrB₃ within the monoclinic phase is a potentially superhard material. The analyses of the electronic density of states and chemical bonding reveal that the strong B–B and B–Zr covalent bonds are attributed to its high hardness. By the quasi-harmonic Debye model, the heat capacity, thermal expansion coefficient and Grüneisen parameter of ZrB₃ are also systemically investigated.

1. Introduction

Superhard materials are considerably used in various industrial applications due to their superior properties of low compressibility, thermal conductivity, refractive index, chemical inertness, and high hardness. Previously, it was generally accepted that the superhard materials are strong covalently bonded compounds formed by light atoms B, C, N, and O, such as diamond [1], c-BN [2], B₆O [3], etc. These materials are expensive and rare, for they can be synthesized only under high-temperature and high-pressure conditions. Recently, it has been reported that the incorporation of light atoms (B, C, N and O) into heavy transition metals (TMs) with high valence electron densities provides a good candidate for designing new hard materials. The compounds formed by TM and light atoms usually possess high valence electron density and directional covalent bonds, and these covalent bonds are strong enough to inhibit creation and movement of dislocations, which greatly improve their mechanical properties and create high hardness. Based on this design criterion, the recent design of new intrinsically potential superhard materials has focused on TM borides, e.g., ReB₂ [4], OsB₂ [5], WB₄ [6,7], CrB₄ [8] and so on. The obtained results exhibit that these materials all have large bulk and shear moduli. Moreover, TM borides can be synthesized under ambient pressure, which results in a low-cost synthesis condition and is beneficial to their application. Thereby, these pioneering studies open up a new route for pursuing new superhard materials.

Up to now, in the Zr-B system, three binary phases [9,10], i.e., ZrB, ZrB₂, and ZrB₁₂, have been synthesized experimentally. Among them, ZrB with the NaCl-type cubic structure has a high melting point and thermal conductivity, and its thermal stability is in the temperature range from 1073 to 1523 K [9]. ZrB₂ is an AlB₂-type hexagonal structure and has properties of both ceramics and metals, and its melting temperature has been reported to be either 3313 [9] or 3523 K [11,12]. Meanwhile, ZrB₂ also possesses high hardness, electrical conductivity, thermal conductivity, and chemical corrosion resistance [13,14]. ZrB₂ ceramic products have been considerably used in making high-temperature structural and functional materials [15,16]. ZrB₁₂ with an UB₁₂-type cubic structure was formed by the peritectic reaction at T = 2303 K [11], and it presents a wide range of interesting characteristics such as high temperature resistance, superconductivity, and large hardness. The crystal structure and related physical properties of ZrB, ZrB₂, and ZrB₁₂ have been investigated theoretically [17,18,19,20,21,22,23,24,25,26]. By using the first-principles calypso algorithm for crystal structure prediction, Zhang et al. [27] predicted two orthorhombic Cmcm and Amm2 structures of ZrB₄, and they also showed that these two phases are potentially superhard materials because of large shear modulus (229 and 236 GPa) and high hardness (42.8 and 42.6 GPa). Very recently, Chen et al. [28] systemically investigated the elasticity, hardness, and thermal properties of ZrB_n (n = 1, 2, 12) by using the first-principles calculations of plane-wave ultra-soft pseudo-potential technology based on the density functional theory. However, there is no experimental report on the synthesis of zirconium triborides; besides, the theoretical research works of zirconium triborides are also few so far.

In this work, the crystal structures of ZrB₃ are extensively investigated by the first-principles particle swarm optimization algorithm (PSO) on crystal structural prediction [29,30,31]. A stable monoclinic C2/m structure has been uncovered for ZrB₃, which is much more energetically preferable than the earlier proposed FeB₃-, TcP₃-, MoB₃-, WB₃-, and OsB₃-type structures in the pressure range of 0–100 GPa. The lattice parameters, total energy, formation enthalpy, phonon frequency, phase stability, elastic properties, density of states, and electronic localization function are then performed to study this novel monoclinic phase by using the first-principles calculations. To further study ZrB₃, the thermodynamic properties are also investigated by the quasi-harmonic Debye model.

2. Computational Methods

To search for potential structure, the PSO technique in crystal structure analysis using the CALYPSO code [32] has been implemented at 0 GPa with one to four formula units (f.u.) in each simulation cell. The underlying ab initio calculations are performed using density functional theory within the Perdew–Burke–Ernzerhof (PBE) generalized gradient approximation (GGA) of the exchange-correlation energy, as implemented in the Vienna ab initio simulation package (VASP) [33]. The all-electron projector augmented wave (PAW) method [34] is employed with 2s²2p¹ and 4s²4p⁶4d²5s² treated as the valence electrons for B and Zr, respectively. Geometry optimization is performed using the conjugate gradient algorithm method with a plane-wave cutoff energy of 520 eV. The calculations are conducted with 11 × 7 × 4, 9 × 11 × 6, 3 × 11 × 10, 8 × 8 × 4, 8 × 8 × 5, and 13 × 13 × 7 for the predicted C2/m-ZrB₃ structure and considered FeB₃-type (No, 11, P2₁/m), TcP₃-type (No, 62, Pnma), MoB₃-type (No, 166, R-3m), WB₃-type (No, 194, P6₃/mmc), and OsB₃-type (No, 187, P-6m2) structures. For hexagonal structures, Γ-centered k mesh is used, and for other structures, appropriate Monkhorst-Pack k meshes [35] are used to ensure that all the total energy calculations are well converged to better than 1 meV/atom. The phonon calculation is carried out by using a supercell approach as implemented in the PHONOPY code [36,37]. Single-crystal elastic constants are calculated by a strain-energy approach, i.e., applying a small strain to the equilibrium lattice and fitting the dependence of the resulting change in energy on the strain. The bulk modulus, shear modulus, Young’s modulus, and Poisson’s ratio are derived from the Voigt–Reuss–Hill approximation [38]. The theoretical Vickers hardness is estimated by the empirical formula proposed by Chen et al. [39]. Moreover, the quasi-harmonic Debye model [40], which is constructed from the Helmholtz free energy at the temperature below the melting point in the quasi-harmonic approximation, is used to obtain thermodynamic properties of ZrB₃.

3. Results and Discussion

Through the PSO technique, we perform a variable-cell structure prediction simulation for ZrB₃ including 1–4 f.u. in the simulation cell at 0 GPa, and we successfully predict a novel monoclinic structure with space group C2/m, as depicted in Figure 1. For the monoclinic C2/m structure, it contains four ZrB₃ in a unit cell with the equilibrium lattice parameters of a = 3.163 Å, b = 5.440 Å, c = 8.773 Å, α = γ = 90°, and β = 93.633°, in which three inequivalent Zr, B1, and B2 atoms occupy the Wyckoff 4i (0.4677, 0, 0.2994), 8j (0.8915, 0.1694, 0.0932), and 4h (0, 0.1668, 0.5) sites, respectively. Figure 1a presents the polyhedral view for the predicted C2/m structure, which shows an intriguing B–Zr–B sandwich stacking order along the crystallographic c-axis. In this structure, B atoms form parallel hexagonal planes, and each Zr atom is coordinated by 12 neighboring B atoms, forming edge-shared ZrB₁₂ hexagonal columns which are connected by boron hexagonal planes. For each ZrB₁₂, the calculated B–B distances are 1.809 (×4), 1.844 (×2), 1.822 (×4), and 1.816 (×2) Å, while the calculated Zr–B distances are 2.729 (×2), 2.563 (×2), 2.434 (×2), 2.472 (×2), 2.528 (×2), and 2.599 (×2) Å.

At zero temperature, a stable crystalline structure requires all phonon frequencies to be positive. Therefore, we calculate the phonon dispersion curves for the predicted C2/m-ZrB₃ at 0 and 100 GPa, respectively. As shown in Figure 2, no imaginary phonon frequency appears in the whole Brillouin zone, indicating its dynamic stability at ambient pressure and high pressure. It is well known that the thermodynamic stability of a compound can be indicated by the energy of its most stable phase. Thus, we calculate the total energy per f.u. as a function of the volume for the predicted C2/m structure, as shown in Figure 3a. For comparison, the previously known five structures of FeB₃, TcP₃, MoB₃, WB₃, and OsB₃ are also considered for ZrB₃. As seen from Figure 3a, our predicted C2/m structure for ZrB₃ has a lower energy minimum than the FeB₃-, TcP₃-, MoB₃-, WB₃-, and OsB₃-type structures. This implies that the predicted C2/m structure is the ground-state phase at 0 GPa. Figure 3b shows the enthalpy curves of ZrB₃ with six different structures relative to the FeB₃-ZrB₃ phase. One can see clearly that the predicted C2/m phase is much more energetically favorable than the previously proposed FeB₃-, TcP₃-, MoB₃-, WB₃-, and OsB₃-type structures in the whole pressure range.

We also calculate the formation enthalpies of the considered structural candidates of ZrB₃ in the pressure range of 0–100 GPa. The formation enthalpy of ZrB₃ with respect to the separate phases can be examined by following the reaction route $\Delta H=H_{{\text{ZrB}}_3}-H_{\text{Zr}}-3H_{\mathrm{B}}$, where ∆H represents the formation enthalpy and the hexagonal Zr and α-B are chosen as the reference phases. Figure 4a presents the calculated formation enthalpies of ZrB₃ with different structures under pressure. It can be clearly seen from this figure that the predicted C2/m structure is the most stable phase of all the considered structures up to 100 GPa. All the considered structures are thermodynamically stable with a negative value of the formation enthalpy below 100 GPa. Among them, the predicted C2/m phase has the lowest formation enthalpy of −2.715 eV at 0 GPa. Therefore, the C2/m phase is more easily synthesized at ambient conditions, and meanwhile this also indicates that the predicted C2/m phase is more stable than the reference phase mentioned above. Figure 4b presents the relative enthalpy-pressure curves of C2/m-ZrB₃ and its respective competing phases with respect to elemental Zr and B. Here we choose ZrB, ZrB₂, and ZrB₁₂ as the competing phases because they have been synthesized in experiments and are thermodynamically stable. From Figure 4b, the predicted C2/m-ZrB₃ is much more energetically stable for decomposing into elements (Zr + B), compounds ZrB + B and ZrB₁₂ + Zr in the pressure range from 0 to 100 GPa, whereas the competing phase ZrB₂ + B is the most stable phase within the given pressure range. It is noteworthy that the relative enthalpy difference between the C2/m-ZrB₃ and ZrB₂ + B becomes smaller and smaller with the increasing pressure. This indicates that the ZrB₃ may be the most stable phase against which to decompose into compounds ZrB₂ + B when the pressure exceeds 100 GPa.

Mechanical stability is a necessary condition for the existence of a crystal, and the mechanical properties (elastic constants and elastic moduli) are important for potential technological and industrial applications. Accurate elastic constants help to understand the mechanical properties and also provide useful information for estimating the hardness of a material. By using a strain-energy method, we obtain the zero-pressure elastic constants (C_ij) of the C2/m-ZrB₃. The calculated elastic constants are listed in

Table 1. The calculated elastic constants C_ij (GPa), bulk modulus B (GPa), EOS fitted bulk modulus B₀ (GPa), shear modulus G (GPa), Young’s modulus E (GPa), G/B ratio, Poisson’s ratio v, Debye temperature Θ_D (K), and hardness H_v (GPa) for the C2/m-ZrB₃.

Structure	Work	C11	C22	C33	C44	C55	C66	C12	C13	C23	B	B0	G	E	G/B	ν	ΘD	Hv
C2/m-ZrB3	this work	576	564	466	207	244	250	45	115	109	238	237	227	517	0.95	0.14	998	42.2
Fm-3m-ZrB	this work	363	-	-	60	-	-	60	-	-	161	160	88	223	0.55	0.27	529	10.5
	theory [17]	282	-	-	111	-	-	100	-	-	161	-	103	254	-	0.24	-	-
	theory [18]	350	-	-	43	-	-	70	-	-	163	-	71	186	-	0.31	-	5.5
	theory [19]	353	-	-	51	-	-	64	-	-	160	-	78	202	-	0.29	-	11.8
	theory [28]	274	-	-	42	-	-	51	-	-	125	-	63	162	0.50	0.28	450	7.9
	experiment [41]	-	-	-	-	-	-	-	-	-	148	-	-	244	-	0.23	-	-
P6/mmm-ZrB2	this work	565	-	432	252	-	255	54	121	-	239	238	233	527	0.98	0.13	930	44.1
	theory [17]	504	-	427	241	-	-	91	112	-	229	-	211	484	-	0.15	-	-
	theory [21]	551	-	436	252	-	-	65	121	-	239	-	229	520	-	0.14	921	-
	theory [28]	502	-	374	225	-	229	43	85	-	199	-	211	468	1.06	0.11	884	43.5
	experiment [20]	581	-	445	240	-	261	55	121	-	245	-	243	554	-	-	-	-
	experiment [42]	-	-	-	-	-	-	-	-	-	220	-	225	502	-	0.11	910
Amm2-ZrB4	theory [27]	554	576	454	223	243	254	50	122	113	241	-	229	522	0.95	0.14 a	-	42.6
Fm-3m-ZrB12	this work	412	-	-	251	-	-	144	-	-	233	233	195	458	0.84	0.17	933	32.5
	theory [26]	413	-	-	244	-	-	141	-	-	232	-	193	453	0.83	0.174	1206	-
	theory [28]	420	-	-	249	-	-	123	-	-	222	-	202	466	0.91	0.15	1231	35.4
	experiment [43]	443	-	-	265	-	-	129	-	-	234	-	-	-	-	0.23	1260	-
	experiment [44]	-	-	-		-	-	-	-	-		-	-	-	-	-	1260	-
Pnma-RuB3	theory [45]	497	462	471	224	219	187	111	229	-	267		185	451	0.69	0.219	-	26.3
R-3m-MoB3	theory [46]	602	-	420	247	-	-	106	166	-	276	277	222	526	0.80	0.18	-	31.8
P63/mmc-WB3	theory [47]	656	-	479	277	-	-	-	-	-	291	-	252	588	-	0.168	-	42
P-6m2-OsB3	theory [48]	525	-	751	186	-	-	125	221	-	317	-	195	485	0.62	0.245		36.9
P21/m-OsB3	theory [45]	674	525	584	135	291	128	97	247	-	303	-	186	463	0.61	0.245		25.1

^a Our calculated result from Reference [27] on the basis of v = (3B − 2G)/(6B + 2G).

along with the theoretical values and available experimental data of other zirconium borides, MoB₃, WB₃, RuB₃, and OsB₃. For a stable monoclinic crystal, the independent elastic stiffness tensor consists of 13 components, C₁₁, C₂₂, C₃₃, C₄₄, C₅₅, C₆₆, C₁₂, C₁₃, C₂₃, C₁₅, C₂₅, C₃₅, and C₄₆, and its mechanical stability criteria can be given by:

$$\begin{array}{c}{C}_{11}>0,C_{22}>0,C_{33}>0,C_{44}>0,C_{55}>0,C_{66}>0,\left[C_{11}+C_{22}+C_{33}+2\left(C_{12}+C_{13}+C_{23}\right)\right]>0,\\ (C_{33}{C}_{55}-C_{35}^2)>0,(C_{44}{C}_{66}-C_{46}^2)>0,(C_{22}+C_{33}-2C_{23})>0,\\ [C_{22}(C_{33}{C}_{55}-C_{35}^2)+2C_{23}{C}_{25}{C}_{35}-C_{23}^2{C}_{55}-C_{25}^2{C}_{33}]>0,\\ g=C_{11}{C}_{22}{C}_{33}-C_{11}{C}_{23}^2-C_{22}{C}_{13}^2-C_{33}{C}_{12}^2+2C_{12}{C}_{13}{C}_{23},\\ \{2[C_{15}{C}_{25}\left(C_{33}{C}_{12}-C_{13}{C}_{23}\right)+C_{15}{C}_{35}\left(C_{22}{C}_{13}-C_{12}{C}_{23}\right)+C_{25}{C}_{35}\left(C_{11}{C}_{23}-C_{12}{C}_{13}\right)]-\\ [C_{15}^2(C_{22}{C}_{33}-C_{23}^2)+C_{25}^2(C_{11}{C}_{33}-C_{13}^2)+C_{35}^2(C_{11}{C}_{22}-C_{12}^2)]+C_{55}g\}>0.\end{array}$$

(1)

As shown in Table 1, the elastic constants of the predicted C2/m-ZrB₃ completely fulfill the elastic stability criteria for a monoclinic crystal, indicating that it is mechanically stable at ambient conditions. The values of C₁₁, C₂₂, and C₃₃ for the new phase are all larger than 400 GPa, suggesting its strong incompressibility along the a, b, and c axis, respectively. Additionally, the values of C₁₁ and C₂₂ of the predicted C2/m phase are much higher than that of C₃₃, which reflects that the bond strengths along the [100] and [010] directions are much stronger than that of the [001] direction. This is because the covalence of the B–B bond is higher than that of the Zr–B bond, as indicated by Figure 8. Thus, the relatively small C₃₃ with regard to C₁₁ and C₁₂ is understandable. C₄₄ is an important indicator for the hardness of a material. The large C₄₄ value (207 GPa) of the predicted C2/m-ZrB₃ indicates its relatively strong strength against the shear deformation. Using the obtained elastic constants, the polycrystalline bulk modulus B and shear modulus G are thus determined by the Voigt–Reuss–Hill approximation. The Young’s modulus E and Poisson’s ratio v are derived by the following formulas: E = 9BG/(3B + G) and v = (3B − 2G)/(6B + 2G). The calculated bulk modulus, shear modulus, Young’s modulus, and Poisson’s ratio of the C2/m phase together with the reference materials mentioned above are listed in Table 1. As seen from Table 1, the predicted C2/m phase has a larger bulk modulus of 238 GPa, which is larger than those of ZrB [17,18,19,28,41] and ZrB₁₂ [26,28,43] but comparable to those of ZrB₂ [17,20,21,28,42], ZrB₄ [27], and other TMB₃ (TM = Mo, W, Ru, Os) [45,46,47,48]. This indicates that the C2/m phase is difficult to compress. Moreover, the bulk modulus (B = 238 GPa) for C2/m-ZrB₃ agrees well with that directly obtained from the fitting results (B₀ = 237 GPa) of the third-order Birch-Murnaghan equation of states, which further verifies the reliability of our elastic calculations. To compare the incompressibility of the C2/m-ZrB₃, other zirconium borides, c-BN, and diamond under pressure, the volume compressions as a function of pressure are shown in Figure 5. Clearly, all the considered zirconium borides except for ZrB have almost the same incompressibility due to their very close bulk moduli, but their incompressibility is lower than those of c-BN and diamond. The shear modulus of a material is a measure of the ability to resist shape change at a constant volume and it plays an important role in hardness compared with the bulk modulus. Interestingly, the shear modulus of the predicted C2/m-ZrB₃ is almost the same as those of ZrB₂ [17,20,21,28,42] and ZrB₄ [27], suggesting that it is a potential superhard material. Except the bulk modulus and shear modulus, the Young’s modulus can also provide a good measure of the stiffness of materials. Generally speaking, the larger the Young’s modulus a material has, the harder it is to deform. Since the Young’s modulus of the C2/m-ZrB₃ is similar to those of the ZrB₂ [17,20,21,28,42] and ZrB₄ [27], it is conceivable that the C2/m-ZrB₃ is a superhard material. Further, Poisson’s ratio v is a crucial parameter to describe the degree of directionality for the covalent bonding. From Table 1, the v value of the C2/m-ZrB₃ (0.14) is a little larger than that of ZrB₂ [42], the same as ZrB₄ [27], and much smaller than those of the other TM borides mentioned above [45,46,47,48], which means there is a strong degree of covalent bonding due to the presence of the planar six-membered ring boron network. To describe the brittleness or ductility of materials, Pugh [49] proposed an experimental criterion by the ratio G/B. The critical value which separates ductile and brittle materials is about 0.57. If G/B < 0.57, the material behaves in a ductile manner; otherwise, the material behaves in a brittle way. It can be seen from Table 1 that the G/B value for the predicted phase is 0.95, implying it is rather brittle nature. All the results strongly support that the C2/m-ZrB₃ is a potential candidate for an ultra-incompressible and superhard material.

The Debye temperature is a very important parameter of materials, and it is closely related to many physical properties such as specific heat, elastic constants, melting temperature, hardness and so on [50]. It is used to differentiate between high- and low-temperature regions for a solid. If the temperature T > Θ_D, all modes are expected to have the energy of k_BT; if T < Θ_D, the high-frequency modes are expected to be frozen, namely the vibrational excitations originate only from the acoustic vibrations. The Debye temperature Θ_D for our studied C2/m-ZrB₃ is estimated from the average sound velocity (v_m) by the following expression:

$${\mathrm{\Theta }}_{\mathrm{D}}=\frac{h}{k}{\left[\frac{3n}{4\mathsf{\pi }}\left(\frac{\mathsf{\rho }{N}_A}{M}\right)\right]}^{\frac{1}{3}}{v}_m$$

(2)

where h is Planck’s constant, k is Boltzmann’s constant, N_A is Avogadro’s number, n is the number of atoms per formula unit, M is the molecular mass per f.u., and ρ is the density. The average sound velocity v_m is given by:

$$v_m={\left[\frac{1}{3}\left(\frac{2}{v_t^3}+\frac{1}{v_l^3}\right)\right]}^{-\frac{1}{3}}$$

(3)

where v_t and v_l are the transverse and longitudinal elastic wave velocities of the polycrystalline material, which is determined by Anderson’s method [51]. Figure 6 presents the dependence of the Debye temperatures of the C2/m-ZrB₃ on the pressure. At zero pressure and zero temperature, our calculated Debye temperature is 998 K for the C2/m-ZrB₃, which is much larger than the values of the known ultra-incompressible RuB₂ (780 K) [52], ReN₂ (735 K) [53], ReO₂ (850 K) [54], and ReB₂ (858.3 K) [55]. Generally, the higher the Debye temperature materials possess, the larger their microhardness. Clearly, the Debye temperature becomes greater with increasing the pressure. This indicates that the pressure is in favor of the improvement of the hardness for the C2/m-ZrB₃.

Due to the high bulk and shear moduli as well as the large Debye temperature for the predicted C2/m-ZrB₃, the hardness calculation is of great importance. The hardness of a material is the intrinsic resistance to deformation when a force is loaded, which depends upon the loading force and the sample quality. The Vickers hardness of a material is estimated by [39]:

$$H_v=2{(K^{2}G)}^{0.585}-3$$

(4)

where G is the shear modulus and K = G/B is the Pugh modulus ratio. Using this model, the estimated Vickers hardness of the C2/m-ZrB₃ is 42.2 GPa, which exceeds the lower limit of superhard materials (40 GPa) and is comparable to the known superhard materials B₆O (45 GPa) [3] and c-BN (48 GPa) [56]. Therefore, the C2/m-ZrB₃ is considered to be a superhard material.

The electronic structure and chemical bonding are key factors to further understand the origin of the mechanical properties of the C2/m-ZrB₃. For this purpose, the density of states (DOS) and bond characteristics are calculated, and the corresponding results are analyzed here. Figure 7a presents the total and partial density of states of the C2/m-ZrB₃ at 0 GPa, where the vertical dashed line denotes the Fermi level (E_F). As seen from this figure, the C2/m-ZrB₃ shows a metallic behavior because of the finite electronic DOS at the Fermi level. From the partial DOS, it exhibits that the peaks below −11 eV are mainly attributed to B-s and B-p states with slight contributions from Zr-s, Zr-p, and Zr-d states. The states above −11 eV largely come from Zr-d and B-p orbitals with small contributions of B-s, Zr-s, and Zr-p. Moreover, the partial DOS profiles of Zr-d and B-p have a very similar shape in the range from −11 to 0 eV, which indicates the significant hybridization between these two orbitals. This hybridization also reflects a strong covalent interaction between the Zr and B atoms. On the other hand, the DOS profile near E_F mainly originates from the 4d state of Zr. Another typical feature of DOS is that there is a deep valley, namely the pseudogap around E_F, which is regarded as the borderline between the bonding and antibonding states [57,58,59]. Significantly, the pseudogap appears below the E_F in the C2/m-ZrB₃, indicating the s-p and p-d bonding states started to be saturated. The nearly full occupation of the bonding states, and without filling the antibonding states, results in the high bulk modulus, large shear modulus, and small Poisson’s ratio, and also increases the structural stability of the C2/m-ZrB₃. To further understand the changes of the total and partial DOS under pressure, Figure 7b–d present the total and partial DOS of the C2/m-ZrB₃ at 20, 50, and 100 GPa, respectively. Compared with the case for 0 GPa (see Figure 7a), the shapes of the total and partial DOS are slightly changed with the increase of the pressure. This indicates that the C2/m-ZrB₃ structure is still stable under pressure up to 100 GPa. Traditionally, the stability of a solid is closely associated with the DOS at the Fermi energy. According to our calculations, the DOS values at the Fermi level are 0.294 eV for 0 GPa, 0.263 eV for 20 GPa, 0.227 eV for 50 GPa, and 0.186 eV for 100 GPa. Therefore, the DOS value at the Fermi level decreases with increasing the pressure. Based on the previous reference report [60], there is a principle stating that the DOS at the Fermi level would be hopefully smaller in an energetically-favored structure. This means that the applied pressure is in favor of the stability of ZrB₃. In addition, one can see from Figure 7 that the total DOS reveals a slight broadening with the enhancement of the pressure. All these results may be explained by the variations of the spacing between atoms.

To gain deeper insight into the bonding character of the predicted C2/m-ZrB₃, we calculate the electronic localization function (ELF), which is based upon a topological analysis related to the Pauli exclusion principle. The ELF is a contour plot in real space, in which different contours have values from 0 to 1. The upper limit ELF = 1 shows the perfect localization property of covalent bonds or lone pairs (filled core levels), whereas the lower limit ELF = 0 is typical for a vacuum (no electron density) or areas between atomic orbitals. ELF = 0.5 corresponds to the perfect free-electron gas, with values of this order meaning regions with bonding of a metallic character. Note that ELF is not a measure of electron density but a measure of the Pauli principle, and it is used to distinguish metallic, covalent, and ionic bonding. Figure 8 presents the calculated ELF contours of the predicted C2/m-ZrB₃ phase on the (001) and (102) planes. Because of the high ELF value between the two adjacent B and B atoms, as shown in Figure 8a, there is the existence of a strong covalent B–B bonding within the planar six-membered ring unit. Meanwhile, one can see from Figure 8b that the large ELF value between the Zr and B atoms implies the partially B–Zr covalent bonding interaction in the C2/m-ZrB₃ phase. Therefore, the strong covalent interaction between B–B bonds and B–Zr bonds is the major reason for its high hardness and stability.

The investigation on the thermodynamic properties of solids at high temperature and high pressure is an interesting topic in materials science. Therefore, we investigate the thermodynamic properties of the C2/m-ZrB₃ in the temperature range from 0 to 1500 K, where the quasi-harmonic model remains fully valid and the pressure effect is also investigated in the pressure range from 0 to 100 GPa. By applying this model, we can get the variations of the lattice heat capacity C_V or C_P, the thermal expansion coefficient α, and Grüneisen parameter γ with temperature and pressure. Figure 9 presents the temperature dependences of the heat capacity at constant volume C_V and the heat capacity at constant pressure C_P at different pressures, respectively. As seen from Figure 9, when the temperature is lower than 400 K, the difference of between C_V and C_P is very small. This can be explained by the relation between C_P and C_V, i.e., C_P = C_V + TVBα² [61], since the difference between the C_P and C_V is mainly determined by α² for the case of low temperature. Therefore, when the thermal expansion coefficient α becomes small, e.g., the temperature T decreases and the pressure P increases, the C_P-C_V difference also becomes small. On the other hand, the C_V and C_P increase rapidly with the temperature at low temperature. This rapid increase is simply due to the exponentially increased number of excited phonon modes. However, when the applied temperature is considerably high, due to the anharmonic effect, the C_P is different from the C_V. The former is proportional to T at high temperature, while the latter reaches a constant value which is the so-called Dulong-Petit limit (C_V(T) ~ 3R for monoatomic solids). Moreover, as the pressure keeps constant, the C_V and C_P increase with the temperature. As the temperature keeps constant, C_V and C_P decrease with the pressure. This shows that the temperature has a more significant influence on the heat capacity than the pressure.

Figure 10 plots the changes of the thermal expansion coefficient α for the C2/m-ZrB₃ with temperature and pressure. As shown in Figure 10a, the thermal expansion coefficient α increases sharply at low temperatures, especially for the case of 0 GPa, and then gradually approaches a linear increase at high temperature, finally becomes moderate at much higher temperature. In terms of α ~ C_V/B [61], the fast increase of α at low temperature is mainly due to that of C_V, and the linear increase of α at high temperature is attributed to that of bulk modulus B (C_V is saturated to the Dulong-Petit limit at high temperature). In addition, the effect of the pressure on α is very small for the case of low temperature, whereas the effect is clearly enhanced at sufficiently high temperature. When the pressure increases, α reduces drastically for a given temperature and the effect of the temperature on it becomes unnoticeable. From Figure 10b, one can see clearly that α is almost close to a constant value at high temperature and pressure. All these results are in good agreement with those of many kinds of materials by applying the Debye theory, such as OsB₄, FeB₄, and TcN [62,63,64].

The Grüneisen parameter γ can describe the anharmonic effects in the vibration of a crystal lattice, and it is commonly used to characterize and explore the thermodynamic behavior of a material at high temperature and high pressure, such as the thermal expansion coefficient and the temperature dependence of phonon frequencies and line-widths. For this purpose, we investigate the dependences of the Grüneisen parameter γ for the C2/m-ZrB₃ with the temperature and pressure, and the corresponding results are given in Figure 11. As seen from Figure 11a, for a given pressure, γ almost remains the same at the low temperature of ≤300 K. As the applied temperature is higher than 300 K, γ tonelessly increases with temperature for the given pressure. In addition, one can find from Figure 11b that γ decreases nearly exponentially with the pressure for a given temperature. As the applied pressure increases, the temperature has a much weaker influence on γ than the pressure.

4. Conclusions

In summary, a novel monoclinic C2/m structure is unraveled to be the ground-state structure for ZrB₃ via the PSO algorithm combined with first-principles calculations. The phonon dispersion and elastic constant calculations have verified that the C2/m-ZrB₃ is dynamically and mechanically stable. The formation enthalpy-pressure relationship shows that this new phase is energetically superior to the previously proposed FeB₃-, TcP₃-, MoB₃-, WB₃-, and OsB₃-type structures within the pressure range of 0–100 GPa. In addition, the high bulk modulus, large shear modulus, and small Poisson’s ratio of the predicted C2/m phase indicate that it is a promising low-compressibility material. Based on the calculated Vickers hardness (42.2 GPa) for the C2/m phase, it is a potential superhard material. The more detailed analyses of the electronic structure and electronic localization function have further demonstrated that the strong covalent B–B bonding and B–Zr bonding are the main driving force for the incompressibility and hardness for C2/m-ZrB₃. By using the quasi-harmonic Debye model, some fundamental thermodynamic properties, such as the heat capacity, thermal expansion coefficient, and Grüneisen parameter are estimated in the pressure range of 0–100 GPa and the temperature range of 0–1500 K. We hope that the present theoretical work will stimulate further experimental research on this material in the future.