Anisotropies in Elasticity, Sound Velocity, and Minimum Thermal Conductivity of Low Borides V x B y Compounds

: Anisotropies in the elasticity, sound velocity, and minimum thermal conductivity of low borides VB, V 5 B 6 , V 3 B 4 , and V 2 B 3 are discussed using the ﬁrst-principles calculations. The various elastic anisotropic indexes ( A U , A comp , and A shear ), three-dimensional (3D) surface contours, and their planar projections among different crystallographic planes of bulk modulus, shear modulus, and Young’s modulus are used to characterize elastic anisotropy. The bulk, shear, and Young’s moduli all show relatively strong degrees of anisotropy. With increased B content, the degree of anisotropy of the bulk modulus increases while those of the shear modulus and Young’s modulus decrease. The anisotropies of the sound velocity in the different planes show obvious differences. Meanwhile, the minimum thermal conductivity shows little dependence on crystallographic direction.


Introduction
Superhard materials have long attracted tremendous attention for a wide range of applications, such as high-temperature applications, surface protection, and abrasive materials [1][2][3][4][5][6][7][8][9][10][11][12][13][14][15]. Diamond is the hardest substance, which maintains records of bulk modulus, shear modulus, and high hardness of 432 GPa, 535 GPa, and 60-150 GPa, respectively. Meanwhile, it is the material with the highest thermal conductivity,~2000 W/m·K. However, the weak thermal stability of diamond at temperatures higher than 800 • C restricts the application in industry [3]. Therefore, the search for and design of new superhard materials with other unique properties is an important topic in the field of materials research. Previous studies of superhard materials have usually focused on short covalent materials [4][5][6], however, pure covalent materials show poor metallicity. The balance of high hardness and excellent metallicity is a central challenge to design superhard metals. However, great success has been achieved with the discovery of transition-metal (TM) light-element (LE) compounds [7][8][9][10][11][12][13][14][15]. It is generally believed that their high hardness stems from a combination of high valence-electron concentration to resist volume compression and strong covalent bonding to counteract shape deformation [7,8]. In addition, the TM-TM bonds of TM-LE compounds provide electron transporting paths in the crystals. Therefore, it is possible for TM-LE compounds to achieve high hardness and electrical (thermal) conductivity at the same time.
Transition-metal borides (TM x B y ) with low boron content reduce the level of covalent bonding and significantly improve the degree of metallic bonding, which improves metallicity of materials, but not for high hardness due to the low resistance against the formation and propagation of dislocations. Surprisingly, the hardness over 40 GPa had 2 of 13 been firstly found in low borides W 0.5 Ta 0.5 B by substituting W with Ta of WB [12]. Subsequently, Liang et al. had reported that a series of transition-metal monoborides exhibited surprisingly anomalous hardening behavior using the first-principles calculations [13]. Li et al. also reported that increasing boron content did not lead to higher hardness in boron-rich tungsten borides [11]. These exciting materials have raised great expectations to design superhard metals.
As candidates for hard materials, vanadium borides show favorable mechanical properties and are quite easily synthesized at normal temperatures and pressures [14][15][16][17]. The studies of low borides are surprisingly scarce, contrary to studies of other boron-rich vanadium borides [15,16]. Liang et al. analyzed the mechanical properties and mechanism of anomalous hardness of VB, V 5 B 6 , V 3 B 4 , and V 2 B 3 [14]. Pan et al. had pointed out that the high elastic modulus and hardness of vanadium borides with low B content are determined by the bond strength and direction [17]. The anisotropy of a crystal reflects the periodicity of the atom arrangement along the different crystallographic orientations in the lattices, which results in the different physical properties of crystals along the different directions. A more detailed exploration of the anisotropic features of V x B y compounds with low boron content is helpful to understand superhard metal. In this paper, the anisotropies of VB, V 5 B 6 , V 3 B 4 , and V 2 B 3 in mechanical properties (bulk modulus, shear modulus, and Young's modulus, sound velocity) and minimum thermal conductivity are discussed.

Method and Computation Details
Four different vanadium borides, including VB (space group: Cmcm), V 5 B 6 (space group: Cmmm), V 3 B 4 (space group: Immm), and V 2 B 3 (space group: Cmcm), were considered in this paper. The primitive cell was used to complete the calculations in the present work. The first-principles method based on density functional theory was used with the Perdew-Burke-Ernzerhof form of generalized gradient approximation for the electronic exchange-correlation interaction implemented in MedeA VASP [18,19]. All calculations were completed with a cutoff energy of plane wave basis~550 eV and Monkhorst-Pack k point meshes ≤0.02 Å −1 [20]. The tolerances of geometry optimization were set as the difference in the total energy being within 5.0 × 10 −6 eV/atom with a maximum force per atom 0.01 eV/Å, a maximum stress of 0.02 GPa, and a maximum displacement of 5.0 × 10 −4 Å. The elastic constants were obtained using the stress-strain method [21] with six strain steps for each optimized structure. Based on the elastic constants, the polycrystalline bulk modulus B and shear modulus G are obtained through the Voigt-Reusse-Hill (VRH) approximation [22]. Young's modulus E is given by Based on the assumption that the indentation size is correlated with the shear modulus G and the width of a formed imprint is proportional to the bulk modulus B, Chen et al. proposed a semi-empirical model to describe the Vickers hardness of materials [23]. This model has been widely used to predict the Vickers hardness of many materials and is given by

Equilibrium Structure and Stability
To estimate the structural stability, the cohesive energy (E coh ) and formation enthalpy (E eff ) are investigated.
where x (y) is the number of V (B) atoms in the unit cell, E total is the total energy of the V x B y compound, E atom V and E atom B are the energies of the free V and B atoms, E g V and E g B are the energy per atom in its ground state. The optimized lattice constants, cohesive energy, and formation enthalpy are listed in Table 1, together with the available theoretical and experimental results. The negative cohesive energy and formation enthalpy indicate that these structures are thermodynamically stable. The stability follows the order VB > V 5 B 6 > V 3 B 4 > V 2 B 3 ; the same result can also be obtained from the position of Fermi level (E F ) of the total density of states (see Figure 1). The structural stability of materials can be judged from the position of Fermi surface in the density of states curve [24]. The higher the value of E F , the lower the stability of the structure. E are the energy per atom in its ground state. The optimized lattice constants, energy, and formation enthalpy are listed in Table 1, together with the available cal and experimental results. The negative cohesive energy and formation entha cate that these structures are thermodynamically stable. The stability follows the > V5B6 > V3B4 > V2B3; the same result can also be obtained from the position of Fe (EF) of the total density of states (see Figure 1). The structural stability of materia judged from the position of Fermi surface in the density of states curve [24]. Th the value of EF, the lower the stability of the structure.

Elastic Properties
The present elastic moduli and Vickers hardness H V data are shown in Table 2. Generally, the compressibility of solids under pressure is expressed by the bulk modulus. The greater the bulk modulus of solids, the less compressible they are. The shear modulus reflects the resistance to shape change under shear force. Young's modulus is related to the stiffness of materials. It can be seen that the bulk, shear, and Young's moduli of V x B y all follow the order V 2 B 3 > V 3 B 4 > V 5 B 6 > VB. Furthermore, the values are very close to each other. Table 2 shows the theoretical Vickers hardness value of V x B y compounds is close to 40 GPa, indicating these V x B y compounds are potential superhard materials [25].

Elastic Anisotropy
Elastic anisotropy has great research value for ceramics due to its correlation with the generation of microcracks and lattice deformation [26,27]. The orientation dependence of three-dimensional (3D) elastic moduli are typically used to describe the elastic anisotropy of a crystal. Considering the effects of shear surface and direction on shear modulus, the average shear modulus over all possible directions is calculated. The formulae of bulk modulus, shear modulus, and Young's modulus of an orthorhombic structure can be obtained from the direction cosines (l 1 , l 2 , and l 3 ) and elastic compliance constant S ij [28], and are given by Figure 2 shows the 3D surface contours and their projections in different planes of the bulk modulus. In this context, a 3D plot is spherical for an isotropic material, and the deviation from a sphere describes the degree of anisotropy. The ellipsoids of the 3D bulk modulus of four structures indicate relatively strong anisotropy in linear compressibility. Furthermore, the 3D contours of VB and V 2 B 3 are similar. The shrinking of the 3D plot depicted in blue indicates the small linear incompressibility along that direction. Based on the ratio of maximum to minimum of directional bulk modulus B max /B min , the order of anisotropy in the bulk modulus is VB < V 5 B 6 < V 3 B 4 < V 2 B 3 . In their planar projections, the nearly circular shape of VB and V 2 B 3 in (100) plane is due to the approximately equal values of C 22 and C 33 . The larger distorted shapes in the (010) and (001) planes indicate its strong anisotropy in linear compressibility. For V 5 B 6 , the similar anisotropy is observed due to the coincident profiles in the (100) and (010) planes, while the anisotropy of the linear compressibility in the (001) plane is very weak due to the small difference of C 11 and C 22 . The anisotropies of the linear compressibility in the (100) and (001) planes are also similar for V 3 B 4 , while that in the (010) plane is also very weak.    Figure 3 shows the 3D surface contours and their projections of shear modulus in different planes. The strongly curved surface shapes suggest strong anisotropy in the shear modulus. Based on the ratio of the maximum to minimum value of the directional shear modulus G max /G min , the order of anisotropy in the shear modulus is VB > V 5 B 6 > V 3 B 4 > V 2 B 3 . For VB, the shear moduli in all planes show relatively strong anisotropy. However, the weakest anisotropy is found in the (100) plane for V 5 B 6 with G max /G min~1 .096, in the (100) and (010) planes for V 3 B 4 with G max /G min~1 .095 and~1.097, and in the (100) and (001) planes for V 2 B 3 with G max /G min~1 .056 and~1.044, respectively. Similarly, the 3D surface contours of the Young's modulus shown in Figure 4 are also distorted, suggesting the strongly anisotropic feature. According to the ratio of maximum to minimum of directional Young's modulus Emax/Emin, the order of anisotropy in the Young's modulus is VB > V5B6 > V3B4 > V2B3. The projections of the Young's modulus in the different planes indicate relatively strong anisotropy, except in the (010) plane for V3B4 with Emax/Emin~1.084 and in the (001) plane for V2B3 with Emax/Emin~1.048. The anisotropy of the Young's modulus reflects a change of strength of the chemical bonds. Similarly, the 3D surface contours of the Young's modulus shown in Figure 4 are also distorted, suggesting the strongly anisotropic feature. According to the ratio of maximum to minimum of directional Young's modulus E max /E min , the order of anisotropy in the Young's modulus is VB > V 5 B 6 > V 3 B 4 > V 2 B 3 . The projections of the Young's modulus in the different planes indicate relatively strong anisotropy, except in the (010) plane for V 3 B 4 with E max /E min~1 .084 and in the (001) plane for V 2 B 3 with E max /E min~1 .048. The anisotropy of the Young's modulus reflects a change of strength of the chemical bonds.  The elastic anisotropy can be estimated by the universal elastic anisotropic index ) and B R (G R ) are the bulk (shear) modulus in Voigt and Reuss approximations) [29] and percent anisotropy in compression [30], respectively. For an elastic isotropic crystal, A U = A comp = A shear = 0 [29,30]. A larger deviation degree from zero means a stronger anisotropy. Figure 5 gives the variation trends of various elastic anisotropic indexes and the ratios of B max /B min , G max /G min , and E max /E min of V x B y . It is obvious that the variation trends of the anisotropy in A comp , A shear , and A U agree well with the results obtained from 3D surface contours of the bulk, shear, and Young's modulus. The compression anisotropy increases, while the shear anisotropy decreases with the increased content of B atoms. The universal elastic anisotropic index includes the contributions of polycrystalline shear and bulk modulus [31]. The elastic anisotropy is mainly determined by the shear anisotropy due to the similar variation tendency between A shear and A U . According to the values A U and E max /E min , the order of elastic anisotropy is VB > V 5 B 6 > V 3 B 4 > V 2 B 3 .
Metals 2021, 11, x FOR PEER REVIEW crystal, AU = Acomp = Ashear = 0 [29,30]. A larger deviation degree from zero means a st anisotropy. Figure 5 gives the variation trends of various elastic anisotropic indexes and tios of Bmax/Bmin, Gmax/Gmin, and Emax/Emin of VxBy. It is obvious that the variation tre the anisotropy in Acomp, Ashear, and AU agree well with the results obtained from 3D contours of the bulk, shear, and Young's modulus. The compression anisotropy inc while the shear anisotropy decreases with the increased content of B atoms. The un elastic anisotropic index includes the contributions of polycrystalline shear and bul ulus [31]. The elastic anisotropy is mainly determined by the shear anisotropy due similar variation tendency between Ashear and AU. According to the values AU and Em the order of elastic anisotropy is VB > V5B6 > V3B4 > V2B3.

Anisotropy in Sound Velocity
Based on the classic theory of transportation of phonon gas solid, the sound v of a crystal structure can be employed to estimate the thermal conductivity [32]. In cific crystallographic plane, the sound velocity consists of three orthogonal modes, transverse mode vt and two mixed modes v±, by solving the Christoffel equation [3 three modes of velocity are given by

Anisotropy in Sound Velocity
Based on the classic theory of transportation of phonon gas solid, the sound velocity of a crystal structure can be employed to estimate the thermal conductivity [32]. In a specific crystallographic plane, the sound velocity consists of three orthogonal modes, a Metals 2021, 11, 577 9 of 13 pure transverse mode v t and two mixed modes v ± , by solving the Christoffel equation [33]. The three modes of velocity are given by ρv 2 t = a 1 sin 2 θ k + a 2 cos 2 θ k (8) 2ρv 2 ± = a 3 sin 2 θ k + a 4 cos 2 θ k ± a 5 sin 2 θ k − a 6 cos 2 θ k 2 + (2a 7 sin θ k cos θ k ) 2 1/2 (9) where the values of a i can be obtained using the elastic constants, and θ k is the angle of wave vector k with respect to one of the symmetry axes. The projections of sound velocity (v t and v ± ) in the different crystallographic planes for V x B y are plotted in Figure 6. Table 3 gives the ratio of maximum to minimum of sound velocity of V x B y . It can be seen that the faster mode v + is larger than pure transverse mode v t and slow mode v -, and the values of v t and vare comparable. The profiles of three acoustic branches (v t and v ± ) of four structures are different in (100), (010), and (001) planes. The ratio of maximum to minimum of sound velocity in different planes (Table 3) indicates the degree of anisotropy. The anisotropies in v t of V 5 B 6 , V 3 B 4 , and V 2 B 3 are relatively weak at (100), (010), and (001) planes, which is mainly determined by the small difference of elastic constants C 44 , C 55 , C 66 . The degree of anisotropy of two mixed modes is determined by the diagonal and off-diagonal elements of the elastic constant matrix.

Anisotropy in Minimum Thermal Conductivity
The minimum thermal conductivity is an important parameter for high-temperature applications of materials, which can be evaluated via the Cahill model [34]. The obtained three sound modes (v t and v ± ) by solving the Christoffel equation are also used.
Here, k B is Boltzmann's constant, n is the number of atoms per molecule, and V is the volume per molecule. The unit of κ min is W/(m·K). Figure 7 gives the planar projections of minimum thermal conductivity in (100), (010), and (001) crystallographic planes for V x B y . It can be found that the differences in profiles of minimum thermal conductivities in different planes are small. Furthermore, these curves show little dependence on crystallographic direction, indicating the anisotropy of minimum thermal conductivity is very weak, which also can be obtained from the ratio of maximum to minimum of minimum thermal conductivity in the different crystallographic planes (Table 4). Based on the Callaway-Debye theory, the minimum thermal conductivity is correlated to the Debye temperature [35]. Thus, the anisotropic feature of minimum thermal conductivity also reflects the anisotropy of the Debye temperature.

Conclusions
The first-principles method is used to investigate the stability, anisotropies in elasticity, sound velocity, and minimum thermal conductivity of low borides VB, V5B6, V3B4, and V2B3. All compounds are thermodynamically stable. The stability follows the order VB > V5B6 > V3B4 > V2B3. The calculated Vickers hardness of VB, V5B6, V3B4, and V2B3 are close to 40 GPa, indicating that these VxBy compounds are potential superhard materials. The elastic anisotropic indexes (AU, Acomp, and Ashear), 3D surface contours, and their planar projections in different crystallographic planes of bulk, shear, and Young's moduli are used to characterize the elastic anisotropy. All VxBy compounds show a relatively strong elastic anisotropic feature, which is mainly determined by the shear anisotropy. The order of elastic anisotropy is VB > V5B6 > V3B4 > V2B3. The anisotropies of sound velocity in different crystallographic planes show obvious differences. The minimum thermal conductivities for all VxBy compounds in different planes are comparable and show very weak anisotropy.

Conclusions
The first-principles method is used to investigate the stability, anisotropies in elasticity, sound velocity, and minimum thermal conductivity of low borides VB, V 5 B 6 , V 3 B 4 , and V 2 B 3 . All compounds are thermodynamically stable. The stability follows the order VB > V 5 B 6 > V 3 B 4 > V 2 B 3 . The calculated Vickers hardness of VB, V 5 B 6 , V 3 B 4 , and V 2 B 3 are close to 40 GPa, indicating that these V x B y compounds are potential superhard materials. The elastic anisotropic indexes (A U , A comp , and A shear ), 3D surface contours, and their planar projections in different crystallographic planes of bulk, shear, and Young's moduli are used to characterize the elastic anisotropy. All V x B y compounds show a relatively strong elastic anisotropic feature, which is mainly determined by the shear anisotropy. The order of elastic anisotropy is VB > V 5 B 6 > V 3 B 4 > V 2 B 3 . The anisotropies of sound velocity in different crystallographic planes show obvious differences. The minimum thermal conductivities for all V x B y compounds in different planes are comparable and show very weak anisotropy.