Fundamental electronic structure and multiatomic bonding in 13 biocompatible high-entropy alloys

High-entropy alloys (HEAs) have attracted great attention due to their many unique properties and potential applications. The nature of interatomic interactions in this unique class of complex multicomponent alloys is not fully developed or understood. We report a theoretical modeling technique to enable in-depth analysis of their electronic structures and interatomic bonding, and predict HEA properties based on the use of the quantum mechanical metrics, the total bond order density (TBOD) and the partial bond order density (PBOD). Application to 13 biocompatible multicomponent HEAs yields many new and insightful results, including the inadequacy of using the valence electron count, quantification of large lattice distortion, validation of mechanical properties with experiment data, modeling porosity to reduce Young’s modulus. This work outlines a road map for the rational design of HEAs for biomedical applications.


INTRODUCTION
High-entropy alloys (HEAs) are complex metallic alloys 1-4 comprised of four, five, or more principal components of different concentrations. Such a structure leads to a high mixing entropy that favors the formation of single phase disordered solid solutions at higher temperatures 5,6 , although the enthalpy also plays a critical role in determining its composition and phase with no long-range-order (LRO) 1,7 . A significant amount of disorder exist in HEAs when they undergo elemental segregation, precipitation, and chemical ordering, but uncertainty remains regarding the existence and the nature of short-range-orders (SRO) 8 . The fundamental theory for the formation of HEAs is not fully established despite many theories and modeling efforts since their inception [9][10][11][12] . Due to the complexity of their compositions and the difficulties for accurate measurements for comparison, most of these efforts are based on different perspectives and quite scattered. They generally used the random solid-solution-model (RSSM) to evaluate their properties although the methods employed could be very different [13][14][15] . Current theoretical methods for HEAs are spread out and include but not limited to the use of CALPHAD 1,12,16 , quasi-random-structure (SQS) [17][18][19] usually in the framework of density functional theory (DFT) on small-size supercells 20 , coherent potential approximations (EMTO-CPA) or effective medium theory 21,22 . Despite a plethora of approaches and methods used, few of them can provide a comprehensive view on their formation and the prediction of their properties due to several obstacles. Firstly, due to relatively large chemical disorder in HEAs, the application of RSSM requires using large supercells. Secondly, the computational demand for accurate large-scale DFT calculations could be prohibitive. Finally, the lack of a non-empirical metric to interpret the results that involve subtle interactions between different metal atoms. By introducing the concept of total bond order density (TBOD) and partial bond order density (PBOD) based on the quantummechanical (QM) metrics (see the "Methods" section), we avoid the use of pure geometric parameters in describing the structures, compositions, and properties of HEAs. The new perspective is based on the understanding the nature of metallic bonding, crucially related to the theory of formation of HEAs. Although metallic bonding has been extensively discussed in the field of metallic glasses (MGs) 23,24 , it has not been thoroughly investigated for HEAs. Metallic bonding is multi-atomic in nature, different from the covalent or ionic bonding where the bond length (BL) is explicitly defined as the distance between the two atoms forming the bond. In MGs, and to a lesser extent in HEAs, the BL can be ambiguous, since all atoms within a certain distance of separation can contribute to metallic bonding. For a fixed distance of separation for a pair of atoms in the model, there could be many possible values of bond order (BO), a measure of the bond strength; and for a specific value of BO, there could be many possible pairs of atoms with the same distance of separation. A theory that is predominately depending on the geometric parameter of BL, or atomic size for interpretation could be oversimplified and problematic. On the other hand, the concept of the TBOD is still applicable as long as the BO values of all the contributing atoms are counted. This point has already been strongly argued for in reference 25 for MG. What differentiate HEAs from MGs is that HEAs have a basic crystal lattice (FCC, BCC, or HCP) as its structural backbone. However, based on the RSSM description, HEAs are not strictly crystalline. They possess no LRO and negligible or small SRO with different nearest neighbor (NN) and second nearest neighbor (SNN) pairs. This is the same predicament facing the vague description of the so-called lattice distortion (LD) that will be discussed later. The use of TBOD and its partial components, the PBOD, to explore the theory of formation in HEAs is clearly a novel approach.
As the populations of developed countries continue to age, there is an increasing demand for biocompatible implant materials. Traditionally, stainless steels and titanium alloys have been typically implemented as joint surrogate metals 26,27 . More recently HEAs 5,7 have been proposed as viable candidates due to their favorable properties, such as high strength and ductility, resistance to corrosion, wear, and fatigue. Importantly, HEAs may consist of refractory elements that are mostly non-toxic and hypoallergenic. They favor the structure of a single solid-solution phase in BCC lattice and have been proposed as a new class of metallic biomaterials 28 . With Mo present, it exhibited excellent biocompatibility compared to the pure Ti 28 . However, the Young's modulus of this alloy system is roughly 10 times greater than human bones 29 . Other challenges including the structural complexity involved at the interfaces with tissues and bones, and other soft matters in an aqueous environment involving body fluids.
In this article, we present the results on the electronic structure, interatomic bonding, and the application of TBOD and PBOD in addressing the challenges for fundamental understanding on the theory of formation of HEAs and its potential applications. We would like to point out the special merits of using TBOD and PBOD as key metrics for assessing the fundamental properties of multicomponent alloys (see "Methods" section for detail). They can be directly compared with each other irrespective of their atomic species, composition, or size. Moreover, they can be applied to other materials systems as long as all interatomic bonding between every pair of atoms are included and normalized by the volume of the cell. This characteristic is very different from other techniques based on ground state energies used in the enthalpy evaluation 19,30,31 . The total energies for different HEAs can be very different because of the difference in the reference energies, although the formation energy for each HEA can explicitly be calculated, which can be quite onerous and time consuming for multi-component HEAs with different compositions. Also, the PBOD resolved from the TBOD have different pairwise components is particularly useful in revealing the details of the interatomic interactions, since TBOD and PBOD are derived from quantum mechanical calculations, not from geometric parameters. Such information cannot be obtained easily based on the calculation of total energies.
In the present work, the electronic structures, interatomic bonding, and mechanical properties of the 13 bioinspired HEAs (Table 1) are investigated through advanced modeling using large supercells yielding many new and insightful results critical to the development and application of biocompatible HEAs.

RESULTS
Electronic structure and interatomic bonding We start the "Results" section with the central part of this paper, the electronic structure, and interatomic bonding. The results for the 13 HEAs constitute a substantial amount of data. Model M3 (TiNbTaZrMo) is chosen as a representative one ( Fig. 1) for focused presentation of the results. They are sketch of the supercell structure of 250 atoms (Fig. 1a), total and partial density of states (DOS) (Fig. 1b), partial charge (PC) distribution (Fig. 1c), and results related to interatomic bonding (Fig. 1d, e). Results for other HEAs are presented in (Figs. S1, S2, and S4). The data on TDOS values at the Fermi level E F (N(E F )) and the minimum values at the locations of the deep valley in the conduction band (CB) above E F of the 13 HEAs are listed in Supplementary Table 1. The similarity in the DOS features amongst the 13 HEAs for energies ranging from −5 to +5 eV is because they are all derived from the 3d/4d/5d orbitals of the transition metals in the HEAs.
An important electronic structural property is the effective charge Q * associated with the so-called PC ΔQ = (Q 0 -Q * ) where Q 0 is the charge of the neutral atom also called the valence electron count (VEC) (see "Methods" section). Figure 1c displays the PC distribution of the 250 atoms in M3. Similar plots for other HEAs are shown in Supplementary Fig. 2. Table 2 lists the averaged PC and Q * for each atomic component in the 13 HEAs. We then average them over the HEAs containing that element (bottom of Table 2 and above the row for VEC). As can be seen, Ti, V, Zr, and Hf have a VEC of 4, Nb and Ta have a VEC of 5, and Mo has a VEC of 6. The calculated Q * values range from 4.10 (M3) to 4.37 (M9) for Ti; 4.09 (M1) to 4.25 (M4) for V; 3.67 (M3) to 3.87 (M9) for Zr, and 3.51 (M10, M13) to 3.72 (M9) for Hf, respectively. These Q* values depend on the compositions of the HEAs. The increases in Q * (Ti, Ta, V, and Mo) or decreases in Q * (Nb, Zr, and Hf) beyond their VEC depend on the atomic compositions and specific HEAs involved. The averaged Q * of each type of atoms for all 13 HEAs show significant deviations from VEC ( Supplementary Fig. 3). This trend strongly suggests that conventional theories using VEC as the key parameter in their formulation may be inadequate and is the main reason for the failure of using the rule of mixture (ROM) in the explanation of simulation results in other studies 32,33 . The trend is also particularly true in some more complicated HEAs with specific local structures, such as those involving defects or other microstructures. Since the entire HEA system should be electrically neutral with a zero net charge, a slight deviation in the atomic Table 1. Fully relaxed structures for the 13 biocompatible models. composition can offer opportunities for fine-tuning their targeted properties. Atom-specific effective charge Q * can also be used as a descriptor in the machine-learning (ML) technique but they are less important than the TBOD and PBOD. The distribution of the BO vs. the bond length ("BL") for all interatomic pairs up to a "BL" of 4.0 Å in M3 is shown in Fig. 1d. The "BL" (with quotation) is the distance of separation between a pair of atoms with a specific BO value. These BOs are calculated with contributions from all atoms near this particular pair since metallic bonding is multiatomic, not pairwise. The TBOD is determined by the strength of all interatomic bonds and their propensity normalized by the cell volume (see "Methods" section). The contributions from the specific PBOD that are decomposed from the TBOD is depicted in the pie chart (Fig. 1e). Similar BO vs. "BL" plots and pie charts for other HEAs are shown in Supplementary Fig. 4.
The TBOD and PBOD in the 13 biocompatible HEAs are presented in the form of histogram bars (Fig. 2a) and pie chart (Fig. 2b). This complex display requires the detailed analysis to appreciate the information that it contains. Some of the general trends are discussed in the figure caption. Still, we cannot be certain about whether there exists any trends related to their composition, atomic species, and their interatomic interactions in relation to PBOD. This critical issue will require a more penetrating analysis. An interesting issue frequently discussed in the HEA community is the possible existence of SRO and if we can quantify it. A commonsense rule is to consider that a single atom has no SRO, or the short-rang-order-parameter SROP = 0, and a diatomic molecule has the maximum SROP = 1. The SROP in HEAs can be obtained from the root-mean-square (RMS) deviation of [PBOD/ TBOD]. First, the ratio of PBOD/TBOD for each of the 13 models are calculated from the data listed in the  34 . It should also be pointed out that the above description is strictly based on the homomgenized RSSM used in the supercell and does not apply to cases where there may involve other factors under differernt assumptions 8,35 . Experimentally, the actual sample will be much larger and may involve multiphase structures, the presence of impurities or dislocations that may induce some form of SRO beyond the simple notion of number of atomic components. The present method of quantifying SRO could still be applied in such a senario as long as the supercell itself contains far more complicated nonhomogeneoius structures beyond the RSSM.
Lattice distortion LD plays a central role in controlling the properties of HEAs 36 . A large LD could be the harbinger of a possible phase transition or sample inhomogeneity. However, the degree of LD in HEAs is difficult to quantify beyond the geometric parameters (Table 1).
We can quantify LD from the BO vs. "BL" data ( Fig. 1d and Supplementary Fig. 4) by resolving them into interatomic pairs. For M3, the scattered plot shows that Ti-Zr, Ti-Ti, and Ti-Nb pairs from NN and SNN groups start to merge ( Supplementary Fig. 4c). This is a strong evidence of large LD. Similar observations for other HEAs are displayed in Supplementary Fig. 4a-h. The investigations on LD reported in the literature are mostly vague and sometimes misleading. The difficulty in providing a more precise description for LD arises from the lack of a basic reference for the undistorted lattice for HEAs due to the random distributions of the different metallic ions of varying sizes throughout the lattice. In other words, there is no such a thing as an undistorted lattice to be used as a reference. To quantify the notion of LD, the histogram distribution of "BL" of interatomic pairs in M3 is plotted (Fig. 1f). The data are fitted with two Gaussians exhibiting a bimodal distribution for the NNN and SNN pairs. The larger the two FWHM and their sum, the greater the LD. Similar analysis applied to all 13 HEAs are shown in Supplementary      3) There is no evidence that the size of the atom influences the PBOD. (4) There is some evidence that the Ti-Ta pair has a relatively large PBOD whereas the Nb-Nb pair has a relatively small PBOD. b Pie chart presentation for PBOD except M3 (in Fig. 1e).

Mechanical properties
One of the most important properties of HEAs is the mechanical properties. They are intimately related to the electronic structure and bonding (see "Methods" section). The calculated results of the 13 HEAs include the elastic coefficients, compliance tensor, bulk modulus (K), shear modulus (G), Young's modulus (E), Poisson's ratio (η), G/K ratio and estimated Vickers harness (Hv) together with the available experimental data [37][38][39][40][41][42][43][44][45][46] is presented in  Supplementary Fig. 6. The Pugh's ratio G/K is well within the range of 0.30 and 0.20 with the average around 0.25 which is on the ductile side for metallic systems, somewhere between bulk metallic glass (BMG) and polycrystalline metals 48 . There is a good correlation between the calculated bulk modulus K and TBOD (Fig. 3), since the TBOD is a single metric that quantifies the internal cohesion. Correlations of TBOD with G and E are less obvious since they involve directional dependence in the strain.

Porosity in TiNbTaZrMo
For the biocompatible HEAs proposed for biomedical applications, the most relevant property is Young's modulus E, which range from 83.7 GPa for M5 (NbTaTiVZr) to 113.5 GPa for M13 (TiNbMoHfTa). To be compatible with the strengths of bones and joints in a human body, the E needs to be reduced by at least 50%. It is well known that there is a correlation between Young's modulus and the porosity in porous materials 49 . We provide some preliminary test calculation by introducing porosity in the simulation for M3 to reduce E. Seven cases of porous models (p1-p6 and px) are constructed from the non-porous model or M3 (p0) by deleting a portion of connected atoms from p0 making sure that the atomic composition is still remain equal partition as in the original model. This is not an easy task but is necessary. The porosity values for each model from p0 to p6 and px are determined by using the PLATON software to be 0.0%, 2.90%, 4.00%, 7.50%, 11.2%, 18.2%, 23.5%, and 30.0%, respectively. These models are shown in Fig. 4a using Crystalexplore17 software. The porosity can reduce Young's modulus significantly in M3 (Fig. 4b). Figure 4a shows that p6 with porosity of 23.5%, the E was reduced to 60.83 GPa from 102.65 GPa in p0, closer to E values exhibited by human cortical bones typically ranging from 7 to 30 GPa 50 , but still way too large. The other porosity-dependent  (1) resonant ultrasound spectroscopy; (2) in-situ neutron; (3) ultrasound spectroscopy; (4) compression; (5) nano-indentation; (6) use Vickers hardness testing. mechanical properties also show similar decrease in bulk and shear modulus with the Poisson's ratio decreases from 0.39 at p0 to 0.30 at p6. This implies that the HEA represented by M3 becomes more brittle as the porosity is increased. However, this point is less clear and requires detailed modeling for biocompatible HEAs in an aqueous environment. Another positive fact is that the introduction of porosity greatly reduce the weight of the materials, an important consideration for biomedical applications.  4 Porosity modeling for M3. a Ball and stick models that exhibit different levels of porosity from p0 to px are plotted using the CrystalExplorer17 software. The different level of porosity is achieved by removing a portion of connected metal atoms in an equal proportion starting with p0 (no porosity) up to p6, which has a large single pore. These models are fully relaxed with the cell volume fixed and their mechanical properties are calculated in the same manner as for the 13 HEA models. b Young's modulus E vs. porosity. The data for porosity values from p1 to p3 contain some fluctuations due to the larger error/volume ratio of the small pores. When the modeling was applied to the model with largest porosity of 34.2% (px), the simulation failed to converge because the ratio of the size of the pore to p0 is too large.
Currently, most biomedical applications involve the use of Ti alloys which is the lightest refractory element in the HEAs.
The current study on porosity simulation is limited to the starting supercell size of 250 atoms. A larger supercell containing 432 atoms (see Supplementary Materials) or even larger and with multiple smaller pores would be necessary for more realistic simulations. Young's modulus for the 432-atom supercell is 100.77 GPa, slightly <102.65 GPa for the 250-atom supercell (Supplementary Table 2). This exploratory study is to provide the proof-of-concept that supercell modeling can provide the guideline to reduce Young's modulus in biocompatible HEAs by introducing porosity of varying shapes that can accommodate bone structures in the human body. In doing so, our simulation method may facilitate the discovery of viable means to fabricate HEAs that are good candidates for biomedical applications.

DISCUSSION
In this study, several new results with lasting impacts are revealed. The method and approach we used here can be readily extended to other multi-phase HEAs, or composites with lighter nonmetallic elements, and the presence of defects or other special refractory elements, such as tungsten (W) 51 greatly expanding its applicability. There is now a new trend to expand HEAs to composite materials, which contain light elements, such as C, B, N, and O (especially C) to form HEA ceramics 32,33 . Although the random solid solution model may still hold, the supercell construction for the composites will be much more challenging because of the involvement of non-metallic elements. Other possible extensions of HEAs include the deviation from equalatomic compositions or even the mixed phases to optimize certain properties. Last but not the least, more realistic modeling for biocompatible HEAs in aqueous solution or body fluid is possible. With the ever-increasing computing power, such large-scale computational modeling in complex materials is very realistic in the near future.

Supercell construction
The cubic supercells for biocompatible high entropy-alloys (HEAs) in the BCC lattice are constructed based on the random solid solution model (RSSM). The size of the supercell, or the number of atoms N it contains BCC supercell is determined by the simple formula N = 2 × (n) 3 . In the present work, n is 5 so the supercell has 250 atoms in the cubic cell of length equal to na, where a is the lattice constant of the single atom BCC crystal of a typical transition metal. We believe that a 250-atom supercell may be the minimal size to justify the use the RSSM with sufficient confidence. It should be pointed out that our supercell is different from those used in the SQS structure in the simulations by other researchers which is usually much smaller. To account for the different possible structural configurations, many such so-called "supercells" have to be used. In the present study, we assert that the statistical distribution of random distribution of metals is sufficient, since the supercell is sufficiently large and with periodical boundary conditions and can account for the random distribution of the NN, second NN and also the third NN and beyond for each atom in the model. These two approaches are similar in spirit but different in strategy. We believe the use of large supercells for many HEAs is more efficient, tractable, and conducive for detailed bond analysis and calculation of TBOD. Four or five atomic species of equal percentage are chosen from the following seven refractory elements: Ti, V, Zr, Nb, Mo, Hf, and Ta, and are distributed randomly in the lattice sites of the supercell with periodic boundaries. For the 5-component HEAs, there are 50 atoms each. For the 4-component HEAs, two of them have 62 atoms, and the other two have 63 atoms each. The initial lattice constant for the supercell is obtained from the scaled average of the crystal lattice constant for each atom. To fully comply with the spirit of RSSM, a script is written such that the atomic occupation of each site is completely random with no restriction to their NN atoms and beyond. So for the 4 (5)-component HEA models, the supercell will have 10 (15) possible NN pairs.

Structural relaxation
The initial BCC supercells for HEAs are fully optimized using the Vienna Ab initio Simulation Package (VASP) 52,53 . VASP is a plane-wave-based DFT method using the pseudopotential. It is very efficient for the structural optimization and elastic properties calculations. In the present work, the PAW-PBE potential for the exchange correlation potential within the generalized gradient approximation (GGA) was used 54 . We adopt a relatively high energy cutoff of 500 eV. The electronic and ionic force convergences were set at 10 −5 eV and 10 −3 eV/Å, respectively, with a 2 × 2 × 2 k-point mesh. Additional testing with selected sample for higher kpoint showed no discernable difference. The final relaxed structures of the 13 biocompatible HEAs (Table 1) based on which all the properties are evaluated. The sketch of the supercell for M3 is shown in Fig. 1a.

Electronic structure and interatomic bonding
For electronic structure and bonding calculations, the in-house developed package, the orthogonalized linear combination of the atomic orbitals (OLCAO) method is used 55 with the VASP-relaxed structures as the input. The OLCAO method is another DFT-based method using atomic orbitals for the basis expansion, which are expressed as a linear combination of Gaussian type of orbitals (GTO). The use of GTO enables the efficient evaluation of three center integrals in the analytic form and thus makes the method highly efficient, especially for large complex systems as demonstrated in many recent publications. In the present calculations for HEAs, a more localized minimal basis (MB) set is used, which consists of the core orbitals and the open shell of valence orbitals. Using Ti as an example, the MB has core orbitals of (1s, 2s, 2p, 3s, and 3p) and valence orbitals of (4s, 4p, and 3d).
In addition to the usual electronic structures such as DOS and partial DOS (PDOS), the most important part of the calculation is the effective charges Q Ã α and the BO values, ρ αβ , between each pair of atoms (α, β) based on the Mulliken population analysis scheme 56 : In Eqs. (1) and (2), S iα,jβ are the overlap integrals between the ith orbital in the αth atom and jth orbital in the βth atom. C m jβ are the eigenvector coefficients of the mth band for the jth orbital in the βth atom. The PC, ΔQ, or the charge transfer for each atom, is the deviation from neutral atomic charge (Q 0 ) from the effective charge (Q * ) or ΔQ = Q 0 -Q * . The accurate determination of PC is extremely important in interpreting many of the material properties and their functionality, especially in HEAs. This function arises from the multi-component nature of the HEAs, which consists of transition metal elements with different d-electron occupations. The BO value, ρ αβ , between a pair of atoms in Eq. (2) provides a quantitative measure of the contribution to the metallic bond from atoms α and β with a specific distance of separation. This BO value is affected by the presence of all the nearby atoms which contribute to the BO. The summation of all BO values normalized by the cell volume gives us the TBOD, which is a single metric to assess the internal cohesion in the crystal 57 . The TBOD can be resolved into partial components or the partial BO density (PBOD) for the different types of atomic pairs, or different groups of atoms in the structural units in the crystal or the supercell.

Mechanical properties
For the elastic and mechanical properties of the HEAs, we used the stress (σ i ) vs. strain (ε j ) response analysis scheme 58,59 on the fully relaxed structure from VASP. A small strain ε (±0.5%) is applied to the supercell to obtain the elastic coefficients C ij and compliance tensor S ij (i, j = 1, 2, 3, 4, 5, W.-Y. Ching et al. 6) by solving the following set of linear equations: From the calculated C ij and S ij , other mechanical properties such as the bulk modulus (K), shear modulus (G), Young's modulus (E), and Poisson's ratio (η) are obtained using the Voight-Reuss-Hill (VRH) polycrystalline approximation [60][61][62] . More detailed information on the methods used for mechanical properties are described in the Supplementary Materials S1.

DATA AVAILABILITY
All the data in this paper including those in the supplementary materials are freely available by contacting the corresponding author (chingw@umkc.edu).

CODE AVAILABILITY
There are several software and packages were used for this study; the analyses were done using Crystalexplore17, Platon, and Origin software. Besides the Vienna ab initio structural Simulation Package (VASP) 53 , the orthogonalized linear combination of the atomic orbitals (OLCAO) 55 method was used.