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

Using the first-principles particle swarm optimization algorithm for crystal structural prediction, we have predicted a novel monoclinic C2/m structure for ZrB3, which is more energetically favorable than the previously proposed FeB3-, TcP3-, MoB3-, WB3-, and OsB3-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 ZrB3 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 ZrB3 are also systemically investigated.


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 6 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 2 [4], OsB 2 [5], WB 4 [6,7], CrB 4 [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 2 , and ZrB 12 , have been synthesized experimentally. Among them, ZrB with the NaCl-type cubic structure has a high melting

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 2 2p 1 and 4s 2 4p 6 4d 2 5s 2 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 3 structure and considered FeB 3 -type (No, 11, P2 1 /m), TcP 3 -type (No, 62, Pnma), MoB 3 -type (No, 166, R-3m), WB 3 -type (No,194, P6 3 /mmc), and OsB 3 -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 3 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 3 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 12

Results and Discussion
Through the PSO technique, we perform a variable-cell structure prediction simulation for ZrB3 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 ZrB3 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   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-ZrB3 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 FeB3, TcP3, MoB3, WB3, and OsB3 are also considered for ZrB3. As seen from Figure 3a, our predicted C2/m structure for ZrB3 has a lower energy minimum than the FeB3-, TcP3-, MoB3-, WB3-, and OsB3-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 ZrB3 with six different structures relative to the FeB3-ZrB3 phase. One can see clearly that the predicted C2/m phase is much more energetically favorable than the previously proposed FeB3-, TcP3-, MoB3-, WB3-, and OsB3-type structures in the whole pressure range. 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 3 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 3 , TcP 3 , MoB 3 , WB 3 , and OsB 3 are also considered for ZrB 3 . As seen from Figure 3a, our predicted C2/m structure for ZrB 3 has a lower energy minimum than the FeB 3 -, TcP 3 -, MoB 3 -, WB 3 -, and OsB 3 -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 3 with six different structures relative to the FeB 3 -ZrB 3 phase. One can see clearly that the predicted C2/m phase is much more energetically favorable than the previously proposed FeB 3 -, TcP 3 -, MoB 3 -, WB 3 -, and OsB 3 -type structures in the whole pressure range.  We also calculate the formation enthalpies of the considered structural candidates of ZrB3 in the pressure range of 0-100 GPa. The formation enthalpy of ZrB3 with respect to the separate phases can be examined by following the reaction route 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 ZrB3 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-ZrB3 and its respective competing phases with respect to elemental Zr and B. Here we choose ZrB, ZrB2, and ZrB12 as the competing phases because they have been synthesized in experiments and are thermodynamically stable. From Figure 4b, the predicted C2/m-ZrB3 is much more energetically stable for decomposing into elements (Zr + B), compounds ZrB + B and ZrB12 + Zr in the pressure range from 0 to 100 GPa, whereas the competing phase ZrB2 + B is the most stable phase within the given pressure range. It is noteworthy that the relative enthalpy difference between the C2/m-ZrB3 and ZrB2 + B becomes smaller and smaller with the increasing pressure. This indicates that the ZrB3 may be the most stable phase against which to decompose into compounds ZrB2 + B when the pressure exceeds 100 GPa.  We also calculate the formation enthalpies of the considered structural candidates of ZrB3 in the pressure range of 0-100 GPa. The formation enthalpy of ZrB3 with respect to the separate phases can be examined by following the reaction route 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 ZrB3 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-ZrB3 and its respective competing phases with respect to elemental Zr and B. Here we choose ZrB, ZrB2, and ZrB12 as the competing phases because they have been synthesized in experiments and are thermodynamically stable. From Figure 4b, the predicted C2/m-ZrB3 is much more energetically stable for decomposing into elements (Zr + B), compounds ZrB + B and ZrB12 + Zr in the pressure range from 0 to 100 GPa, whereas the competing phase ZrB2 + B is the most stable phase within the given pressure range. It is noteworthy that the relative enthalpy difference between the C2/m-ZrB3 and ZrB2 + B becomes smaller and smaller with the increasing pressure. This indicates that the ZrB3 may be the most stable phase against which to decompose into compounds ZrB2 + B when the pressure exceeds 100 GPa. We also calculate the formation enthalpies of the considered structural candidates of ZrB 3 in the pressure range of 0-100 GPa. The formation enthalpy of ZrB 3 with respect to the separate phases can be examined by following the reaction route ∆H = H ZrB 3 − H Zr − 3H 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 3 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 3 and its respective competing phases with respect to elemental Zr and B. Here we choose ZrB, ZrB 2 , and ZrB 12 as the competing phases because they have been synthesized in experiments and are thermodynamically stable. From Figure 4b, the predicted C2/m-ZrB 3 is much more energetically stable for decomposing into elements (Zr + B), compounds ZrB + B and ZrB 12 + Zr in the pressure range from 0 to 100 GPa, whereas the competing phase ZrB 2 + B is the most stable phase within the given pressure range. It is noteworthy that the relative enthalpy difference between the C2/m-ZrB 3 and ZrB 2 + B becomes smaller and smaller with the increasing pressure. This indicates that the ZrB 3 may be the most stable phase against which to decompose into compounds ZrB 2 + 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 (Cij) of the C2/m-ZrB3. The calculated elastic constants are listed in Table 1 along with the theoretical values and available experimental data of other zirconium borides, MoB3, WB3, RuB3, and OsB3. For a stable monoclinic crystal, the independent elastic stiffness tensor consists of 13 components, C11, C22, C33, C44, C55, C66, C12, C13, C23, C15, C25, C35, and C46, and its mechanical stability criteria can be given by: (1) 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 3 . The calculated elastic constants are listed in Table 1 along with the theoretical values and available experimental data of other zirconium borides, MoB 3 , WB 3 , RuB 3 , and OsB 3 . For a stable monoclinic crystal, the independent elastic stiffness tensor consists of 13 components, C 11 , C 22 , C 33 , C 44 , C 55 , C 66 , C 12 , C 13 , C 23 , C 15 , C 25 , C 35 , and C 46 , and its mechanical stability criteria can be given by: (1)

Structure
Work As shown in Table 1, the elastic constants of the predicted C2/m-ZrB 3 completely fulfill the elastic stability criteria for a monoclinic crystal, indicating that it is mechanically stable at ambient conditions. The values of C 11 , C 22 , and C 33 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 11 and C 22 of the predicted C2/m phase are much higher than that of C 33 , 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 33 with regard to C 11 and C 12 is understandable. C 44 is an important indicator for the hardness of a material. The large C 44 value (207 GPa) of the predicted C2/m-ZrB 3 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 12 [26,28,43] but comparable to those of ZrB 2 [17,20,21,28,42], ZrB 4 [27], and other TMB 3 (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 3 agrees well with that directly obtained from the fitting results (B 0 = 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 3 , 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 3 is almost the same as those of ZrB 2 [17,20,21,28,42] and ZrB 4 [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 3 is similar to those of the ZrB 2 [17,20,21,28,42] and ZrB 4 [27], it is conceivable that the C2/m-ZrB 3 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 3 (0.14) is a little larger than that of ZrB 2 [42], the same as ZrB 4 [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 3 is a potential candidate for an ultra-incompressible and superhard material.
Debye temperatures of the C2/m-ZrB3 on the pressure. At zero pressure and zero temperature, our calculated Debye temperature is 998 K for the C2/m-ZrB3, which is much larger than the values of the known ultra-incompressible RuB2 (780 K) [52], ReN2 (735 K) [53], ReO2 (850 K) [54], and ReB2 (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-ZrB3.  Due to the high bulk and shear moduli as well as the large Debye temperature for the predicted C2/m-ZrB3, 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]: 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 B T; 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 3 is estimated from the average sound velocity (v m ) by the following expression: 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: 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 3 on the pressure. At zero pressure and zero temperature, our calculated Debye temperature is 998 K for the C2/m-ZrB 3 , which is much larger than the values of the known ultra-incompressible RuB 2 (780 K) [52], ReN 2 (735 K) [53], ReO 2 (850 K) [54], and ReB 2 (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 3 . Due to the high bulk and shear moduli as well as the large Debye temperature for the predicted C2/m-ZrB 3 , 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]: 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 3 is 42.2 GPa, which exceeds the lower limit of superhard materials (40 GPa) and is comparable to the known superhard materials B 6 O (45 GPa) [3] and c-BN (48 GPa) [56]. Therefore, the C2/m-ZrB 3 is considered to be a superhard material.  Due to the high bulk and shear moduli as well as the large Debye temperature for the predicted C2/m-ZrB3, 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]: The electronic structure and chemical bonding are key factors to further understand the origin of the mechanical properties of the C2/m-ZrB 3 . 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 3 at 0 GPa, where the vertical dashed line denotes the Fermi level (E F ). As seen from this figure, the C2/m-ZrB 3 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 3 , 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 3 . 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 3 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 3 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 3 . 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.  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-ZrB3 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 CV or CP, 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 CV and the heat capacity at constant pressure CP at different pressures, respectively. As seen from Figure 9, when the temperature is lower than 400 K, the difference of between CV and CP is very small. This can be explained by the relation between CP and CV, i.e., CP = CV + TVBα 2 [61], since the difference between the CP and CV is mainly determined by α 2 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 CP-CV difference also becomes small. On the other hand, the CV and CP 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 To gain deeper insight into the bonding character of the predicted C2/m-ZrB 3 , 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 3 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 3 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 3 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α 2 [61], since the difference between the C P and C V is mainly determined by α 2 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.  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-ZrB3 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 CV or CP, 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 CV and the heat capacity at constant pressure CP at different pressures, respectively. As seen from Figure 9, when the temperature is lower than 400 K, the difference of between CV and CP is very small. This can be explained by the relation between CP and CV, i.e., CP = CV + TVBα 2 [61], since the difference between the CP and CV is mainly determined by α 2 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 CP-CV difference also becomes small. On the other hand, the CV and CP 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 CP is different from the CV. The former is proportional to T at high temperature, while the latter reaches a constant value which is the so-called Dulong-Petit limit (CV(T) ~ 3R for monoatomic solids). Moreover, as the pressure keeps constant, the CV and CP increase with the temperature. As the temperature keeps constant, CV and CP 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-ZrB3 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 α ~ CV/B [61], the fast increase of α at low temperature is mainly due to that of CV, and the linear increase of α at high temperature is attributed to that of bulk modulus B (CV is saturated to the Dulong-Petit limit at high temperature). In addition, the effect of the pressure on α is very small for  Figure 10 plots the changes of the thermal expansion coefficient α for the C2/m-ZrB 3 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 4 , FeB 4 , and TcN [62][63][64].

Conclusions
In summary, a novel monoclinic C2/m structure is unraveled to be the ground-state structure for ZrB3 via the PSO algorithm combined with first-principles calculations. The phonon dispersion and elastic constant calculations have verified that the C2/m-ZrB3 is dynamically and mechanically stable. The formation enthalpy-pressure relationship shows that this new phase is energetically superior to the previously proposed FeB3-, TcP3-, MoB3-, WB3-, and OsB3-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-ZrB3. 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. 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 3 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.

Conclusions
In summary, a novel monoclinic C2/m structure is unraveled to be the ground-state structure for ZrB3 via the PSO algorithm combined with first-principles calculations. The phonon dispersion and elastic constant calculations have verified that the C2/m-ZrB3 is dynamically and mechanically stable. The formation enthalpy-pressure relationship shows that this new phase is energetically superior to

Conclusions
In summary, a novel monoclinic C2/m structure is unraveled to be the ground-state structure for ZrB 3 via the PSO algorithm combined with first-principles calculations. The phonon dispersion and elastic constant calculations have verified that the C2/m-ZrB 3 is dynamically and mechanically stable. The formation enthalpy-pressure relationship shows that this new phase is energetically superior to the previously proposed FeB 3 -, TcP 3 -, MoB 3 -, WB 3 -, and OsB 3 -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 3 . 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.