Structural, elastic, mechanical and thermodynamic properties of HfB4 under high pressure

The present work aims to study the structural, elastic, mechanical and thermodynamic properties of the newly discovered orthorhombic Cmcm structure HfB4 (denoted as Cmcm-HfB4 hereafter) under pressure by the first-principles calculations. The obtained equilibrium structure parameters and ground-state mechanical properties were in excellent agreement with the other theoretical results. The calculated elastic constants and phonon dispersion spectra show that Cmcm-HfB4 is mechanically and dynamically stable up to 100 GPa and no phase transition was observed. An analysis of the elastic modulus indicates that Cmcm-HfB4 possesses a large bulk modulus, shear modulus and Young's modulus. The superior mechanical properties identify this compound as a possible candidate for a superhard material. Further hardness calculation confirmed that this compound is a superhard material with high hardness (45.5 GPa for GGA); and the relatively strong B–B covalent bonds’ interaction and the planar six-membered ring boron network in Cmcm-HfB4 are crucial for the high hardness. Additionally, the pressure-induced elastic anisotropy behaviour has been analysed by several different anisotropic indexes. By calculating the B/G and Poisson's ratio, it is predicted that Cmcm-HfB4 possesses brittle behaviour in the range of pressure from 0 to 100 GPa, and higher pressures can reduce its brittleness. Finally, the thermodynamic properties, including enthalpy (ΔH), free energy (ΔG), entropy (ΔS), heat capacity (CV) and Debye temperature (ΘD) are obtained under pressure and temperature, and the results are also interpreted.


Introduction
Ultra-incompressible and hard materials have attracted a great deal of attention owing to their outstanding physical and chemical properties in fundamental science and industry applications [1][2][3][4]. Previously, it was commonly accepted that superhard materials are those strong three-dimensional covalent compounds formed by light elements (B, C, N and O), such as diamond [5], c-BN [6], BC 2 N [7], B 6 O [8] and BC 5 [9] etc. However, these superhard materials are infrequent and expensive because they can be synthesized only under extreme high-temperature and high-pressure conditions. To overcome the shortcoming, many research activities attempt to synthesize new intrinsically superhard materials by combining small strongly covalent bonding light elements and large electron-rich heavy transition metals (TMs). The compounds formed by TMs and light atoms usually possess high valence electron density and directional covalent bonds, which contribute to improve their mechanical properties and hardness. Based on this design principle, a series of ultra-incompressible and hard materials were successfully synthesized [6,[10][11][12][13][14][15][16]. Recently, considerable attention has been paid to transition metal borides (TMBs), such as Re-B [3], Os-B [17], Ru-B [18] and W-B [19,20]. Especially the newly boron-rich TM compounds, such as WB 4 [21,22], CrB 4 [23,24], MnB 4 [25][26][27], FeB 4 [28,29], YB 4 [30] and VB 4 [31], have been extensively investigated. This may be due to boron having a strong ability to form covalent bonds with TMs. A lot of calculated results also show that all these materials exhibit excellent mechanical properties and high hardness. In addition, TMBs can be synthesized at ambient condition, which reduces the cost of synthesis and is beneficial for practical application. Therefore, in order to find new superhard materials, further studies are needed for the TMBs, especially the boron-rich TM compounds that have been proposed recently. The study will open up a new way for seeking new ultra-incompressible and hard materials.
In the boron-rich TM compounds, the studies of Hf-B compounds are surprisingly scarce, compared with the extensive research works on TMB 4 (TM = W, Cr, Mn and Fe) [21][22][23][24][25][26][27][28][29]. For the Hf-B system, three known binary phases have been synthesized experimentally [32], namely HfB, HfB 2 and HfB 12 . For HfB, earlier studies [33,34] suggested that it has an orthorhombic FeB structure and its thermal stabilization is in the range from 1200°C to 2000°C. A later experiment reported that HfB also has a NaCl-type phase [32]. With regard to HfB 2 , it possesses a hexagonal AlB 2 -type structure [32], exhibits many unique properties, such as high chemical and thermodynamic stability, high hardness and low electrical resistivity; these characteristics make HfB 2 become a promising candidate for various hightemperature structure materials [35]. For the HfB 12 , its cubic structure was experimentally synthesized at 6.5 GPa and 1600-2100°C [36]. Recently, an experimental study showed that the cubic structure of HfB 12 can be successfully stabilized in the Y 1−x Hf x B 12 system under ambient pressure [37]. Most recently, a novel promising phase, orthorhombic Cmcm structure HfB 4 , has been proposed by Zhang et al. [38] with the particle swarm optimization (PSO) algorithm. The calculated results show that the Cmcm-HfB 4 may be a potential promising superhard material and energetically much more stable than the previously proposed YB 4 -, ReP 4 -, FeB 4 -, CrB 4 -and MnB 4 -type structures at zero pressure. However, up to now, there are few works in the literature reporting on the physical properties of Cmcm-HfB 4 under pressure (to our knowledge). Moreover, many fundamental aspects are still not well understood because of its complex structural and chemical behaviour. Therefore, a detailed study of the physical properties of Cmcm-HfB 4 is undoubtedly necessary and may bring new insight for further understanding its unique mechanical properties.
In this study, we studied the structural, elastic and mechanical properties of the novel Cmcm-HfB 4 under ambient and high pressures by the first-principles approach based on density functional theory (DFT). The effects of pressure and temperature on the thermodynamic properties of HfB 4 are systematically obtained. This paper is organized as follows. In §2, the methods used are described. In §3, the main results of our calculations of structural and mechanical properties under pressure are given and compared with other TMBs. The computed thermodynamic properties are also reported in this section. Finally, in §4, the conclusions of this work are given.

Computational methods and details
In the current work, all theoretical calculations were carried out by the first-principles planewave pseudo-potential method based on DFT in Cambridge Serial Total Energy Package (CASTEP) code [39]. The contribution of core electrons is described using the Vanderbilt non-local ultrasoft pseudopotential [40]. Pseudo-atom calculations are performed for B: 2s 2 2p 1 and Hf: 5p 6 5d 2 6s 2 . To compare the performance of different approximations of exchange-correlation interaction, here we used the generalized gradient approximation (GGA) with the Perdew-Burke-Ernzerhof (PBE) functional [41] and the modified PBE functional (PBESOL) [42], which was proved to provide good results for solids of high density [43], as well as the local density approximation (LDA) with the form of Ceperley-Adler parametrized by Perdew & Zunger (CAPZ) [44], used to describe the exchange-correlation potentials. The structures were relaxed using the Broyden, Fletcher, Goldfarb and Shannon (BFGS) minimization method algorithm [45]. The lattice constants and atom coordinates were optimized by minimizing the total energy. Through a series of convergence studies with respect to cut-off energies and k-points, the cut-off energies were set at 500 eV, and the K-space integration over the Brillouin zone was carried out using a 5 × 8 × 2 k-point Monkhorst-Pack mesh for the compounds. The values of the kinetic energy cutoff and the k-point density were determined to ensure the convergence of total energies to within 0.01%. In geometrical relaxations, the residual force is less than 0.001 eV Å −1 , the maximum displacement of atoms to less than 5.0 × 10 −4 Å and the convergence threshold for the maximum stress to 0.02 GPa. All these parameters were tested to ensure that the self-consistent total energies converged to within 5.0 × 10 −7 eV atom −1 .
The elastic constants C ij , needed for the calculation of mechanical properties and to study the mechanical stability of Cmcm-HfB 4 , were calculated through C ij = σ i /ε j , where σ and ε are the elastic stress and strain, respectively. The subscripts i and j denote the Cartesian coordinates of the considered structures [46]. For the optimization of the internal atomic freedoms, the criteria for convergence were set as follows: the difference in total energy was within 1.0 × 10 −6 eV atom −1 , maximum force was 0.002 eV Å −1 , maximum displacement was 1.0 × 10 −4 Å and maximum stress within 100 GPa. Based on these elastic constants and Voigt-Reuss-Hill (VRH) approximations, the bulk modulus (B), shear modulus (G), Young's modulus (E) and Poisson's ratio (ν) were derived by the formulae: The effective bulk and shear moduli are also introduced by B = (B V + B R )/2 and G = (G V + G R )/2, where the subscripts V and R denote the Voigt and Reuss approximation, respectively.
The phonon dispersion curves and phonon density of states (DOS) were determined by using the linear response density functional perturbation theory [47,48]. The phonon frequencies in the Brillouin zone (BZ) centre are computed as second-order derivatives of the total energy with respect to atomic displacements. In addition, some phonon-related thermodynamic properties such as the enthalpy ( H), entropy ( S), free energy ( G), lattice heat capacity (C v ) and Debye temperature (Θ D ) are evaluated in a quasi-harmonic approximation. 3. Results and discussion 3.1

. Structural properties under pressure
The crystal structure of HfB 4 is orthorhombic with the symmetry of the space group Cmcm (no. 63, figure 1a-c). The Cmcm-HfB 4 contains four HfB 4 formula units (f.u.) in its unit cell, in which two inequivalent Hf and B atoms occupy the Wyckoff 4c (0, 0.4190, 0.75) and 16h (0.8312, 0.8767, 0.5795) position, respectively. Surprisingly, in Cmcm phase, each metal Hf atom is surrounded by 12 neighbouring B atoms, and B atoms form parallel hexagonal planes, which results in the formation of HfB 12 hexagonal columns (as shown in figure 1b). These structural units form intriguing B 6 -Hf-B 6 sandwiches stacking order along the crystallographic c-axis, which exhibit ultra-incompressible characterization. In the HfB 12 polyhedron, the three types of Nb-B bond distances are calculated to be 2.448 (× 4), 2.505 (× 4) and 2.612 (× 4) Å. In the quadrangle B ring, the shortest B-B bond distance is 1.7791 Å, which is smaller than that of the known hard materials ReB 2 (1.81 Å), OsB 2 (1.85 Å) and Os 2 B 3 (1.876 Å), indicating that the Cmcm-HfB 4 structure may have strong stability and good mechanical properties.
To obtain a stable geometry structure, we first make the structural relaxation and optimization with the two different exchange-correlation approximations (GGA and LDA), to determine the internal atomic coordinates and structure parameters. The obtained lattice parameters, density, unit cell volume and atomic positions within the GGA and LDA are listed in table 1 together with the available results and similar TMBs for comparison. From table 1, it can be seen that the obtained equilibrium lattice constants, density, cell volumes and atomic positions at the GGA/LDA level are also in excellent agreement with the theoretical values [31,49]. This good agreement shows the accuracy of our calculations. In addition, the results from GGA (PBE and PBESOL) have a better agreement with the previous data compared to the LDA (CAPZ), owing to the obtained lattice parameters by the two exchange correlation functions           (PBE and PBESOL) being almost same for the Cmcm-HfB 4 compound. Therefore, in this work, we only choose the GGA-PBE method in the following calculations.
As is well known, pressure can change the structural stability and mechanical behaviour of a compound during its synthesis process. To provide some insight into the pressure behaviour of HfB 4 , the normalized parameters a/a 0 , b/b 0 , c/c 0 and volume V/V 0 as a function of pressure are shown in figure 2, where a 0 , b 0 , c 0 and V 0 are the equilibrium structural parameters at zero pressure, respectively. From figure 1a, it is obvious that the values of the structural parameters decrease with increasing pressure, which means that the bond length of Cmcm-HfB 4 becomes shorter as the pressure increases. Meanwhile, the curves also present that the crystal cell of Cmcm-HfB 4 along the crystallographic b-axis is much more difficult to compress than in the other directions in the whole range of pressure. This may be due to the fact that there is very strong local covalent bonding between Hf atoms and B atoms. In addition, figure 2 shows that pressure effects will generally reduce the volume and increase the density of Cmcm-HfB 4 . Moreover, in the pressure range studied here, no transition between the different phases occurs.

Elastic and mechanical properties under pressure
The elastic properties (e.g. elastic constants, elastic moduli and elastic anisotropy etc.) under pressure are of key importance for us to understand the deformation behaviour of a solid in response to external forces, and can provide a deeper insight into macroscopic mechanical behaviour and help estimate the hardness of materials. Therefore, an explicit determination of the elastic properties of Cmcm-HfB 4 as a function of pressure is essential, and will have an important guidable significance to accelerate the synthesis of Cmcm-HfB 4 . The predicted zero pressure elastic constants C ij , nine for the orthorhombic phase (C 11 , C 22 , C 33 , C 44 , C 55 , C 66 , C 12 , C 13 and C 23 ) of Cmcm-HfB 4 are determined by the strain-stress method. These results are listed in table 2 along with other typical TMBs (TMB 4 , TM = Y, Fe, Cr, Mn, Ta, Zr, Mo, Re, V and Os) for comparison. From table 2, our results are in reasonable agreement with those obtained by Zhang et al. [38]. The good agreement between the current and previous results supports the accuracy and reliability of our elastic calculations.
Next, we consider the mechanical stability, which is a necessary condition for the stability of crystal. For an orthorhombic crystal, the requirement for mechanical stability obeys the following restrictions on the elastic constants: [55] 23 ) > 0 and [C 11 + C 22 + C 33 + 2(C 12 + C 13 + C 23 )] > 0. Our                      This suggests that Cmcm-HfB 4 may be a potential superhard material. In addition, the Young's modulus is also an important parameter for measuring the stiffness of a material in addition to the bulk modulus and shear modulus. The larger the Young's modulus a material has, the harder it is to deform. From table 3, the Young's modulus (535 GPa for GGA and 577 GPa for LDA) of HfB 4 is larger than the other compounds [31,[49][50][51][52], indicating that HfB 4 has a higher hardness than the other TM compounds. All of these excellent mechanical properties strongly suggest that m-HfB 4 is a potential candidate for a superhard material.
Owing to the large elastic constants, and high bulk and shear moduli for Cmcm-HfB 4 , its hardness calculation is of great importance. It denotes the resistance to elastic and plastic deformations when a force is loaded. Here, we predicted the hardness of Cmcm-HfB 4 , which can be correlated with the product of the squared Pugh's modulus ratio and the shear modulus. The Vickers hardness is estimated by [56] H where K represents the G/B ratio and G is the shear modulus. Using the model, the estimated Vickers hardness H v of Cmcm-HfB 4 is in good agreement with a previous study [38]. The high hardness of Cmcm-HfB 4 probably stems from the relatively strong B-B covalent bonding and the presence of a planar boron six-membered ring unit in Cmcm-HfB 4 structure.
In addition, one should be aware that the values presented in tables 2 and 3 were obtained at P = 0 GPa, and that pressure effects will generally increase the elastic constants and mechanical modulus. To further explore the influence of pressure on the elastic and mechanical properties of Cmcm-HfB 4 , one can use the relationship between the elastic constants and pressure shown in figure 3. Unfortunately, there are no literature data to compare with our predicted results for the pressure derivative of elastic properties of Cmcm-HfB 4 . Therefore, our results can serve as a prediction for future studies. From figure 3a,b, clearly, the elastic constants and the elastic modulus monotonically increase under pressure. The elastic constants are positive and satisfy the well-known Born criteria [57], indicating that Cmcm-HfB 4 is a mechanically stable phase up to P= 100 GPa. From figure 3a, it is clear that the calculated C 22 value is higher than those of C 11 and C 33 , implying that the resistance to deformation along the b-axis is stronger than that along the a-and c-axis. This just corresponds with the above-determined fact that it is much more difficult to be compressed along the b-axis than along the a-and c-axis under pressure. In addition, C 44 is an important indicator for the hardness of a material. The large C 44 value (247 GPa for GGA and 270 for LDA) indicates its relatively strong strength against shear deformation. From figure 3b, it can be seen that the elastic moduli monotonically increase under pressure, suggesting that Cmcm-HfB 4 is more difficult to compress with increase in pressure.
To further probe into the physical properties of Cmcm-HfB 4 , we analysed its ductility and brittle nature under pressure according to the Pugh criterion [58]. A material with a high G/B ratio value (greater than 0.57) is associated with brittle nature; otherwise, materials should have ductility. Our results show that the G/B ratio changes from 0.99 (P = 0 GPa) to 0.67 (P = 100 GPa). These values indicate that Cmcm-HfB 4 shows behaviour as a brittle material in the whole pressure range from 0 to 100 GPa (figure 3d). The relative directionality of the bonding in the material also has an important effect on its hardness and can be determined from the G/B ratio. HfB 4 at 0 GPa has the largest G/B ratio, which suggests the strongest directional bonding between the Hf and B atoms, and this bonding decreases as the pressure increases.
The calculated Poisson ratio (ν) is an additional argument used to describe the ductility and brittle nature of materials, with 0.33 as the critical reference point. For brittle metallic materials, these values are small, whereas for ductile metallic materials, the value is bigger than 0. 33. In this study, ν changes from 0.13 (P = 0 GPa) to 0.22 (P = 100 GPa); therefore, Cmcm-HfB 4 shows behaviour as a brittle material and then it tends to become ductile (figure 3d). In addition, the ν value is far below 0.33, suggesting a high degree of covalent bonding, which contributes to the material's hardness. The Poisson's ratio of Cmcm-HfB 4 (0.128 for GGA) is smaller than that of the previously reported hard material ReB 2 (0.171), indicating that the degree of covalent bonding of Cmcm-HfB 4 is very good because of the presence of a planar six-membered ring boron network. The Debye temperature Θ D as a function of pressure is shown in figure 3c. Moreover, the calculated Θ D value shows a linear increase with increase in pressure due to stiffening of the crystal structure.
To further understand the mechanical properties of Cmcm-HfB 4 under pressure, we calculated the DOS at 0 and 100 GPa, and the corresponding results are shown in figure 4. The Fermi energy level was taken as the origin of the energy. As is clearly shown in figure 4, Cmcm-HfB 4 presents a metallic character. There is a deep valley and significant pseudo-gap at E F , which indicates the covalent interaction between Hf and B of Cmcm-HfB 4 . In addition, the pseudo-gap is widened with the increase in pressure, revealing that the bond length is reduced under pressure, while the covalent interaction between Hf and B is increased with the increase in pressure. Therefore, this is also the reason for the high elastic modulus and the high hardness of Cmcm-HfB 4 .

Elastic anisotropy under pressure
Considering that the microcracks are easily induced in materials owing to significant elastic anisotropy especially for hard materials, it is worthwhile to investigate the elastic anisotropy of materials. The shear anisotropic factors provide a measure of the degree of anisotropy in the bonding between atoms in different planes. The shear anisotropic factor for the {100} shear planes between the 011 and 010 .   Table 4. The pressure (in GPa) dependence of the shear anisotropic factors A 1 , A 2 , A 3 , anisotropy of bulk modulus A B and A G (in %), and compressibility anisotropy factors A Ba and A Bc of Cmcm-HfB 4 .      13 .

directions is
For the {010} shear planes between the 101 and 001 directions, it is 23 .
And for the {001} shear planes between the 110 and 010 directions, it is 12 . (3.4) For orthorhombic crystals, the elastic anisotropy arises from the anisotropy of the linear bulk modulus in addition to the shear anisotropy. The anisotropies of the bulk modulus along the a-axis and the c-axis with respect to the b-axis can be written as where B a , B b and B c are the bulk moduli along different crystal axes, defined as For an isotropic crystal, these factors must be 1.0. Any departure from 1.0 is a measure of the elastic anisotropy of the crystal. A more practical measure of elastic anisotropy is defined as follows: where B and G denote the bulk and the shear modulus, and the subscripts V and R represent the Voigt and Reuss approximations, respectively. A value of zero is associated with elastic isotropy, while a value of one indicates the largest possible elastic anisotropy. The calculated anisotropic factors under different pressure are collected in table 4. For Cmcm-HfB 4, A B = 0.18%, A G = 0.8% at 0 GPa, indicating that it is nearly isotropic. The percentage of bulk modulus anisotropy is smaller than the percentage of shear modulus anisotropy under pressure, suggesting that there is less anisotropic incompressibility than in shear. In addition, the compressibility anisotropy factor A B b is nearly equal to 1.0  become closer, and their interactions become stronger to play an important role to form a superhard material. These mechanical properties also imply the anisotropy of HfB 4 .
To better understand the anisotropic behaviour, we obtained the 3D surface representations of the direction-dependent Young's modulus and its plane projections at 0 and 100 GPa, as shown in figure 5ad). From figure 5, one can discover that the (001) plane exhibits isotropic property but the (100) and (010) planes exhibit anisotropic behaviour at 0 GPa; these results are consistent with those obtained by Zhang et al. [38] In addition, from figure 5c,d, it is not difficult to find that the anisotropy is significant with increase in pressure, which is in accordance with the above analysis from the anisotropy factor.

Thermodynamic properties under pressure
At present, most of the studies on the relative stability of polymorphs mainly consider the minimum lattice energies at 0 K and 0 bar condition. [59,60] However, considering the influence of the pressure and temperature on the structure, it is necessary to investigate the thermodynamic properties, such as the Gibbs free energy ( G), which can be used to determine the relative stability of polymorphs at desired thermodynamic conditions [43]. In addition, theoretical studies of the thermal properties are complementary to experiments and are useful for extending the regions of phase space which cannot be reached experimentally. Moreover, considering the pressure can change the mechanical stability and the reaction kinetics of a compound during its synthesis process, an explicit determination of its thermodynamic properties under pressure is essential. Figure 7a shows the enthalpy, free energy and entropy under pressure. The Helmholtz free energy, F = E -T*S (where F is the free energy) is the relevant potential in an ensemble, where the volume and temperature are independent variables. As shown in figure 7a, below 100 K, the values of the enthalpy, free energy and entropy are almost zero. Above 100 K, the free energy at 0.0 GPa decreases   faster with increasing temperature with respect to 100 GPa, and the TS at 0.0 GPa increases rapidly as the temperature increases with respect to 100 GPa, resulting in a linearly increasing relationship between the variations of enthalpy. Figure 7b shows the variations of the lattice heat capacity C V with temperature. We observe that the C V increases with the applied temperature but decreases under pressure. Below