Computational exploration of biomedical HfNbTaTiZr and Hf0.5Nb0.5Ta0.5Ti1.5Zr refractory high-entropy alloys

Refractory high entropy alloys (RHEAs) have been proven to be a potential candidate in the biomedical field due to their balanced mechanical properties and biocompatible composition. Recent experimental findings show that RHEAs like HfNbTaTiZr and Hf0.5Nb0.5Ta0.5Ti1.5Zr have good mechanical properties such as high polarization and wear resistance than others which establish them as potential materials for biomedical application. In this work, we performed first-principles density functional theory calculations on the mechanical and thermal properties of HfNbTaTiZr and Hf0.5Nb0.5Ta0.5Ti1.5Zr. The predicted lattice constant, density, Young’s modulus, and Vickers hardness are consistent with the available experimental report, which verifies the accuracy of the applied model. The thermal coefficient of linear expansion of both RHEAs has been investigated by utilizing the Debye theory. The present methods could be applied to study other future RHEAs on exploration of their physical properties.


Introduction
The use of biomaterials in implantation has attracted the interest of researchers in the biomedical field due to their attractive physical properties such as high strength, biocompatibility, corrosion resistance, and wear resistance [1][2][3][4]. The commonly used biomedical alloys includes 316L, stainless steel, CoCrMo, Ti6Al4V, and titanium-based alloys [5][6][7][8][9]. However, these materials are very sensitive to body fluids, which can easily corrode and release toxic metallic ions (Ni, Cr, Co and Al) when used for a long term [7,[10][11][12][13][14]. For example, titanium alloy, Ti6Al4V is promising and is widely used for implant applications, but it releases ions of V and Al which affects the organs and tissues of patient's body. Recent study shows the aluminum ions is a possible cause for Alzheimer disease, while vanadium oxide is a toxic chemical for human body [15][16][17]. This results in early failure of the implants system which is not desirable. The demand for implants is increasing in developed countries to maintain the quality of life of aged people. Therefore, it is necessary to design metallic alloys with balanced mechanical properties, high hardness, high corrosion resistance, and biocompatibility [18,19].
Refractory high entropy alloys (RHEAs) are advanced novel class metallic alloys that are comprised of five or more different elements, each element has a concentration ranging from 5 to 35 at.% [20,21]. As compared to traditional metallic alloys, RHEAs consist of alloying refractory elements such as titanium (Ti), zirconium (Zr), molybdenum (Mo), niobium (Nb), iron (Fe), vanadium (V), and tantalum (Ta). RHEAs are considered to be favorable refractory alloys with non-allergic and non-toxic properties [22][23][24][25]. Therefore, compared to traditional metallic alloys, they have great potential application in biomedical field [26].
HfNbTaTiZr was extensively studied for high temperature application in the past [27][28][29][30][31][32]. Recently Amir et al [33] synthesized equimolar HfNbTaTiZr and non-equimolar Hf 0.5 Nb 0.5 Ta 0.5 Ti 1.5 Zr RHEAs by arc melting method and investigated their microstructure, mechanical, corrosion and wettability properties for possible biomedical application. They compared the properties of two RHEAs, HfNbTaTiZr and Hf 0.5 Nb 0.5 Ta 0. 5 [34][35][36][37][38][39][40]. However, most experimental studies focus on the structural and mechanical properties with limited study on other thermodynamics properties. For the real application of RHEAs, detailed investigation of electronic, structural, mechanical, and thermodynamic properties is very important. Moreover, computational methods are economical and easier to investigate the alloys' structural, mechanical, and thermal properties as compared to experimental techniques.
In this report, the comparative study of electronic, mechanical, and thermodynamics properties of two RHEAs, HfNbTaTiZr and Hf 0.5 Nb 0.5 Ta 0.5 Ti 1.5 Zr, are performed by using first-principles density functional theory (DFT) calculations. The accuracy of the predicted properties of both RHEAs was further validated by comparing with available experimental data. It has been observed that the calculated mechanical properties are in good agreement with the available experiments. The quasi-harmonic Debye-Grüneisen theory was adopted to compute the thermodynamics properties.

Computational method
All mechanical properties were calculated using Vienna Ab-Initio Simulation Package (VASP) 5.4 [41] in MedeA software [42]. The generalized gradient approximation (GGA) [43] of Perdew-Burke-Ernzerhof (PBE) [44] is used as an exchange-correlation function. The plane wave cut off 500 eV was used in all calculations. The convergence criterion is 0.02 eV/Ang. The electronic iterations convergence is 10 −5 eV using the blocked Davidson algorithm. The applied k-spacing of 0.2 per Angstrom, leads to a 4×3×3 mesh in our calculation. This gives a k-spacings of 0.185×0.197×0.197 per Angstrom. The first order Methfessel-Paxton smearing with a width of 0.2 eV was used. A first-principles calculation is performed to predict the mechanical properties of RHEAs. Previously, this method has been successfully applied for studying the alloying effects of element on structural, mechanical, and thermal properties of RHEAs [45][46][47][48]. Hu et al [45] studied the hydrogen behaviors and microstructural development in the HEA TiZrHfMoNb by applying first-principles method.
The random structure of BCC supercell is constructed by utilizing the python code and Knuth shuffle model [49] consisting of 100 atoms and is shown in figure 1. Virtual crystal approximation (VCA) [50] and special quasi-random structure (SQS) [51] are often used methods for simulating random distribution of RHEAs components on a given crystal lattice. However, VCA methods do not consider the effect of lattice distortion that exists in complex RHEAs whereas the SQS method considers the lattice distortion effect by combining the idea of cluster expansion with large supercells [52]. For complex RHEAs, large supercells are required to resemble the random distribution of elements. As an example, Tian et al [53] have constructed a large supercell of 250 and 128 atoms using the Alloy Theoretic Automated Toolkit to study the elastic properties of various RHEAs. The major limiting factor of the SQS model is the high computational cost associated with a large supercell. The proposed 100 atoms random supercell model is capable of distributing the components of the complex alloy randomly and its convergence test has been succinctly described in a previous publication [54]. This 100 atoms supercell is appropriate for providing highly accurate results with low computational costs. The stress-strain method [55,56] which is implemented in MedeA software is used to calculate the elastic constants C , 11 C , 12 and C . 44 VASP computes the stress tensor by applying analytic expressions. Several strains are used to get more points for the fitting procedure involved in the calculation of the elastic constants. The equilibrium supercell is optimized with high accuracy to avoid zero strain. The Voight-Reuss-Hill approximation [57] was used to calculate bulk modulus (B), Shear modulus (G), Young's modulus (E), and Poisson ratio (n).
The formula to compute B, G, E, and n are Voight 11 12 44 and S is compliance matrix, The Grüneisen-Debye approximation was applied to compute the thermal expansion coefficient and Grüneisen constant G [58]. The internal energy-volume equation of state provided by Mayer et al [59] was used to find the Grüneisen constant G where V 0 is equilibrium volume and B is bulk modulus. After calculating the γ G , Debye temperature is calculated as, where V 0 is volume, q represents atoms in the unit cell, , and B K are Planks and Boltzmann constant, respectively. The specific heat capacity, C , V as a function of temperature, T, is estimated as, where R is ideal gas constant; x and x i j are the atomic percentages of the ith and jth element, respectively; r i is the radius of the ith element;r is averaged atomic radius; and ( ) VEC i is the valence electron concentration of ith element.
From the report of Zhang et al [61] and Sheng et al [64], the required conditions for HEAs to form a solid stable solution are: 11ΔSmix 19.5 (J/K.mol), −22ΔHmix7 (kJ mol −1 ), and Ω>1.1 and δ<6.6%. The calculated parameters are shown in table 1, which suggest that two current RHEAs will form a solid solution structure as it satisfies the conditions proposed by references [61,62]. According to Guo et al [63], if the calculated VEC of alloy is lower than 6.67, it will form a BCC crystal structure and both RHEAs have VEC<6.67   indicating the existence of BCC phase in them. The experimental reports [27,29,31,33] show only BCC phase exists for both RHEAs, which confirms the VEC could be the indicative parameter to identify the phase in alloys. The lattice constant of both RHEAs was calculated using the relaxed super cell volume and results were compared with available experiment data which is shown in figure 2. The calculated lattice constant of HfNbTaTiZr is 3.403 Å which is equal to the experimental values of 3.40 Å [27, 29] and 3.405 Å [31]. Similarly, the calculated lattice constant of Hf 0.5 Nb 0.5 Ta 0.5 Ti 1.5 Zr was equal to 3.391 and its experimental value is 3.405 [33]. It is seen that the changes in the atomic concentration of alloys from equimolar to non-equimolar lowers the atomic size of alloys as shown in table 1. This lowering of atomic size is responsible for lattice contraction which slightly decrease the lattice constant of Hf 0.5 Nb 0.5 Ta 0.5 Ti 1.5 Zr [65].

Mechanical properties
To explore the mechanical properties of two RHEAs, we obtained the elastic constants which are important parameters in determining the stability of materials and calculated values are shown in figure 3. According to the Born Rule, the mechanical stability criteria required for materials are 44 0, and | | > C C . 11 12 The values of the elastic constants C , 11 C , 12 and C 44 meet the Born stability criteria for both RHEAs, showing that their crystal structures are mechanically stable. The experiment conducted by Dirras et al [32] found the C 11 and C 44 of HfNbTaTiZr are 172±6 GPa and 28±1.5 GPa, respectively, which is closer to our simulation value. Similarly, Fazakas et al [66] calculated the value of C 11 and C 44 of HfNbTaTiZr by firstprinciples method, which are 160 GPa and 62.4 GPa, respectively. The first-principles calculations overestimate the C 44 value and underestimate the value of C . 11 The elastic constants C 11 and C 12 of HfNbTaTiZr are larger than Hf 0.5 Nb 0.5 Ta 0.5 Ti 1.5 Zr. The reason for higher elastic constant of HfNbTaTiZr could be the lattice distortion brought by size difference in atomic radii of elements [67]. However, there is no significant changes in the C 44 value for both alloys and the existing first-principles method is not able to provide a conclusive explanation. The Cauchy pressure (C 11 -C 44 ) of materials is useful in predicting the nature of bonding which, if positive, is likely to form a metallic bond with ductile behaviors [68]. It is found that the value of Cauchy pressure for both alloys are positive which indicates the formation of metallic bonds exists in them and they are ductile in nature. The anisotropy factor (A) [69] is an excellent parameter that reveals how much anisotropy a cubic crystal possesses which is calculated by the following equation 44 11 12 For isotropic crystals, the anisotropy factor must equal one, and any value of A less or more than one relates to the crystal's degree of elastic anisotropy. The calculated A from elastic constants is 0.87 and 0.91 for HfNbTaTiZr and Hf 0.5 Nb 0.5 Ta 0.5 Ti 1.5 Zr RHEAs, respectively. As a result, one can conclude that Hf 0.5 Nb 0.5 Ta 0.5 Ti 1.5 Zr is somewhat more anisotropic than HfNbTaTiZr. Furthermore, the higher value of elastic constants C 11 and C 12 increases the other mechanical properties like bulk modulus (B), shear modulus (G), Young's modulus (E) and Vickers's hardness (H v ) in HfNbTaTiZr. The theoretically calculated density (ρ), B, G, E and H v for both RHEAs are shown in table 2. The Tian's et al [70] model was used to calculate the Vickers hardness (H v ) as The hardness of the HfNbTaTiZr was found to be 2.92 GPa, which is higher than the Hf 0.5 Nb 0.5 Ta 0.5 Ti 1.5 Zr. The increased hardness of HfNbTaTiZr could be due to the cocktail effect [72] which can result in unexpected properties changes in metallic alloys. For metallic alloys, the cocktail effect can be achieved by changing the composition and alloying of materials. According to Pugh's rule [73], ductile material will have a B/G ratio >1.75 and brittle material will have B/G<1.75. The calculated B/G ratio for both RHEAs is greater than 1.75, which suggests that these RHEAs will have a ductile nature. The calculated Young's modulus of Hf 0.5 Nb 0.5 Ta 0.5 Ti 1.5 Zr is 86 GPa which is lower than the HfNbTaTiZr. Due to the lower elastic modulus, Hf 0.5 Nb 0.5 Ta 0.5 Ti 1.5 Zr could be a potential alloy for bone implants [71]. Furthermore, slightly higher value of poison ratio of Hf 0.5 Nb 0.5 Ta 0.5 Ti 1.5 Zr indicates its better plastic property [74]. Higher plastic material is easy to handle than brittle material during the manufacturing process. There is a good agreement between the calculated elastic constants, Young's modulus, shear modulus, Poison's ratio (ν), Vickers hardness with the available experimental reports which provides validation of the theoretical analysis. This demonstrates that the random supercell model merged with the DFT-based calculation is a feasible approach to predict the elastic constants and other mechanical properties correctly.

Thermal properties
It is very necessary to study the thermal stability of alloys at different temperatures in order to use them as a high temperature structural material. Therefore, the thermal properties of both alloys have been studied. Debye temperature (θ D ) has been calculated which helps to identify the natural covalent bond of solid materials. It is found that the material with larger θ D has a higher chance of forming a strong covalent bond [75]. The calculated θ D , Grüneisen parameter, and V m of HfNbTaTiZr are 235.7 K, 2.25, and 2136 m s −1 , respectively and those of Hf 0.5 Nb 0.5 Ta 0.5 Ti 1.5 Zr are 242.9 K, 2.32, and 2205 m s −1 . The thermal coefficient of linear expansion (a L ) of two RHEAs is shown in figure 4. It is a very useful parameter for explaining how temperature can change the size of materials and its value has wide application in determining the structural as well as mechanical applications. It can be seen from the diagram that for Hf 0.5 Nb 0.5 Ta 0.5 Ti 1.5 Zr, the a L value is almost the same as HfNbTaTiZr up to 70 K and becomes higher after 100 K. The higher value of a L shows its sensitivity nature towards high temperature. Since a L is inversely proportional to the melting temperature of materials, Hf 0.5 Nb 0.5 Ta 0.5 Ti 1.5 Zr has a higher melting temperature than HfNbTaTiZr which results in increased a L as temperature increased. The θ D and a L of Hf 0.5 Nb 0.5 Ta 0.5 Ti 1.5 Zr is greater than the HfNbTaTiZr which implies it has better thermal properties. There is no experimental and theoretical study of thermal properties related to both RHEAs as per the authors' findings. Therefore, simulations of the thermal properties of current RHEAs could provide important guidance for future experimental and ab-initio investigation. It is clear from this computational study that changing the composition of elements of RHEAs could affect the mechanical and thermodynamic properties, often described by four core effects [76]. In this study, various structural, mechanical, and thermal properties of RHEAs have been explored and compared with available experiments. More experimental study is needed to validate some of our computational findings, and it is expected to be confirmed in the near future.

Conclusion
The first-principles DFT calculations have been performed to study the structural, mechanical, and thermal properties of two RHEAs, HfNbTaTiZr and Hf 0.5 Nb 0.5 Ta 0.5 Ti 1.5 Zr, with BCC structure. The calculated elastic constants indicate that both RHEAs are mechanically stable. The calculated lattice constant, density, Young's modulus, and hardness for both RHEAs are consistent with the experiments. It is found that the Hf 0.5 Nb 0.5 Ta 0.5 Ti 1.5 Zr RHEA have low Young's modulus, high plasticity, and better thermal properties which is likely to be a potential candidate for future biomedical application. The agreement of our computational findings with the available experimental data confirms the efficiency and accuracy of 100 atoms supercell model which is capable of predicting the mechanical properties precisely. Therefore, DFT computational technique could be an alternative cost-effective and efficient method, besides related experiments, to explore the structures, mechanical, and thermodynamical properties of the novel RHEAs that have potential application in biomedical fields.