Subtle Variations of the Electronic Structure and Mechanical Properties of High Entropy Alloys With 50% Carbon Composites

Binary and ternary transition metal carbides are stable ceramic crystals with outstanding mechanical properties. In recent years, multicomponent single-phase high entropy alloys enjoyed explosive growth due to many of their outstanding physical properties with its large and flexible composition space. Hence the composite between them can be advantageous in forming a new class of ceramic materials with combined superiority in their properties for many applications. On the basis of a systematic large-scale ab initio simulations using density functional theory that are specifically designed for their compositional variations, subtle differences in their structures, electronic and mechanical properties are revealed and discussed in detail. Fifteen supercell models with 512 atoms with equal composition of C and high entropy alloys occupying the sub-lattices of the rock-slat structure were constructed. These models are fully optimized, and their properties carefully characterized, compared and contrasted. By applying the novel concept of total bond order density and its partial components, the partial bond order density, we revealed many subtle variations in their properties that have not been known before. This large database can play an important and valuable role in the design and synthesis of high entropy ceramic carbides.


INTRODUCTION
Following the pioneer work less than two decades ago (Senkov et al., 2013;Sun et al., 2017;Ding et al., 2018), high entropy alloys (HEAs) is now a flourishing research area in metal alloys with no insight. HEAs are complex multicomponent homogeneous single crystalline alloys (Gao et al., 2016;Lim, 2016;Ye et al., 2016;Miracle and Senkov, 2017) comprised of four, five or more principal components in near equal concentrations. They possess many outstanding physical properties including high strength and ductility (Senkov et al., 2011;Gali and George, 2013;Gludovatz et al., 2014;Lim, 2016;Ye et al., 2016;Li et al., 2016;Li et al., 2017;Tsao et al., 2017), hardness (Senkov et al., 2010;Senkov et al., 2013;Ye et al., 2016), superconductivity (Von Rohr et al., 2016) and much more. These novel alloys have high entropy of mixing thus favor the formation of single phase disordered solid solutions at higher temperatures Li et al., 2016). Although the enthalpy certainly plays a critical role in determining its composition and phase stability in the absence of long-range-order (Miracle and Senkov, 2017;Santodonato et al., 2018), uncertainty remains regarding the possible existence of short-range-orders and the nature of their origin . Disorder do exists in HEAs during experimental synthesis when they encounter precipitation, undergo elemental segregation, and variation in chemical ordering at different temperature. The fundamental theory for the formation of HEAs is still not fully established since their inception (Senkov et al., 2013;Sun et al., 2017;Ding et al., 2018) due to the complexity of their compositions, difficulties encountered in synthesis, precision in experimental measurements, and proper comparison with modeling efforts. Most of these efforts are based on different perspectives and techniques. They are quite scattered, which is one of the main reasons for HEAs to be a flourishing and highly crowded area of contemporary alloy research.
Transition metal (TM) carbides have been a significant portion of ceramic materials since early 1930s (Agte and Alterthum, 1930). They are binary or ternary crystals with well-defined crystalline structures and some of them are known for exceptionally high melting temperature e.g. Ti 3 C 2 Si and Ta 4 HfC 5 (Agte and Alterthum, 1930). The binary TM carbides have the FCC rock-salt structure with two sub-lattices and strong TM-C bonds. It is therefore obvious that the composites formed between multicomponent HEAs and carbon can be especially advantageous in forming a new class of ceramic materials with combined superiority in their properties such as increased thermal stability, enhanced mechanical properties, corrosion and oxidation resistance for special applications especially at ultrahigh temperature (Sarker et al., 2018;Harrington et al., 2019). The term for HEA-carbide is essentially to generalize the same TM-C bond as in the respective cases in the hope of additional advantages offered by multicomponent TMs. The first successful demonstration of the HEA-ceramics is by Rost et al. (2015). This is followed by the recent work on (Hf-Ta-Zr-Nb)C by Dusza et al. (2018) and Zhou et al. (2018) on (Ti-Zr-Hf-Nb-Ta)C.
It should be pointed out that there are many recent works in which a small amount of C, usually less than 5 atomic %, are added to mostly FCC HEAs such the as the standard Canton alloy CrMnFeCoNi (Stepanov et al., 2016;Cheng et al., 2017;Wang et al., 2017;Ikeda et al., 2019) to investigate the resulting local structural variations, formation of microstructures and precipitates and potential enhancement of their properties. These works are completely different from the HEAs-carbides composites in which C constitute 50% of the lattice sites.
Computational modelling has been an important component in HEA research and more recently in HEA-C (Sarker et al., 2018;Zhang et al., 2019) and other HEA-ceramics (Yang et al., 2018). They use different method and approaches follow the accepted concept of the random solid-solution-model (RSSM) for HEAs (Tian, 2017;Widom, 2018;LaRosa et al., 2019). These theoretical and computational methods for HEAs are widely spread out and include the use of CALPHAD (Senkov et al., 2013;Miracle and Senkov, 2017;Sanchez et al., 2019), special quasi-randomstructure (SQS) (Zunger et al., 1990;Curtarolo et al., 2013;Yang et al., 2018), coherent potential approximations (EMTO-CPA) (Tian et al., 2016), effective medium theory (Huang et al., 2018), and different levels of density functional theory (DFT) (Chen et al., 2018) on small-size supercells. Despite a plethora of methods used, few of them can provide a comprehensive view on their formation and the prediction of their properties. Recently, we have reported a systematic study of 13 biocompatible HEAs in BCC lattice consists of refractory TMs (Ching et al., 2020). We used large supercells of 250 atoms each and the random solidsolution-model protocol. The results include in-depth analysis of their electronic structures, interatomic bonding, and predict HEA properties based on the use of the novel quantum mechanical metric, the total bond order density (TBOD) and the partial bond order density (PBOD) resulting in many new insights. This work has special merits because of 1) using large supercells; 2) accurate and large-scale DFT calculations; and 3) detailed interpretation of the results on the subtle interactions between different metal atoms.
We now extend the above ab initio simulation on pure HEAs to HEA-carbide composites in the FCC rock salt (NaCl) structure that retain the short-range-orders in binary carbides. We specifically designed a total of 15 supercell models with 50% of C and 50% of HEA among nine different TM elements to test the data-driven concept of simulating a large number of models and identify subtle differences. In these models, TMs and C occupied different sub-lattice sites. There are two types of interatomic bonding. The nearest neighbors (NN) TM-C bonding and the second NN TM-TM bonding with larger distance separations.
It should be pointed out that we used a completely different philosophy and method from the recent work of Sarker et al. (2018) although the goals are the same. Our main conclusion is the total verification of the concept the HEA-carbides retain the same short range order as in the respectively binary carbides with new twist to adjust and improve the mechanical properties by extending to multi-component HEA of TMs. The main focus is on model m1 (Ti 10 Nb 10 Ta 10 Zr 10 Mo 10 C 50 ) shown in Figure 1. Other models essentially test the compositional variations from m1. Information on the structure and properties relationship at atomic scale of HEA carbide ceramics is much less available in comparison with the conventional HEAs. A large database on DFT calculations that can be used for future predictions is very valuable. In other words, computational modelling on HEA-C materials can help in exploring and developing novel HEA ceramics.

MODEL CONSTRUCTION AND METHODS
The supercells are constructed from (4 × 4 × 4) of FCC rock salt (NaCl) lattice which has 512 lattice sites. The sublattice A or the 256 anion sites are fully populated by C atoms and the sub-lattice B or the 256 cation sites are randomly occupied by five TM elements in equal or nearly equal proportion. The refractory TM elements involved in these models are Ti, V, Cr, Zr, Nb, Mo, Hf, Ta, and W, They are distributed randomly in nearly equal percentage to construct two 4-component, seven 5-component, and one 6-component HEAs in the cation sub-lattice. The 4-TM component HEA-C contains 64 atoms each, the 5-TM components contains one TM with 52 and the rest with 51 atoms, and the 6-TM components contains two TMs with 42 atoms, and the other four with 43 atoms each. It is important to note that the TM components are distributed randomly in the cation sublattice sites of the supercell with periodic boundaries. The carbon atoms are in the other sub-lattice have different NN TMs. The supercell we constructed for the 15 HEA models differ from the quasi-random-structure SQS structure used in other studies) (Zunger et al., 1990;Curtarolo et al., 2013;Yang et al., 2018), and also claimed to be supercells but are usually much smaller in size.
In order to reveal the subtle differences due to different compositions and atomic components in these 15 models listed in Table 1, we consciously divide them into six groups (G1-G6) for specific comparison and focused discussion. G1: (TiNbTaZrMo) 50 C 50 . This consists of a single model m1 which is central to this paper since all other models are morphed from m1 with specific modifications.

Structural Relaxation
The Vienna ab initio simulation package (VASP) (Kresse and Furthmüller, 1996;Kresse and Furthmuller, 2002), which is highly effective for structure optimization, was employed to relax the initial model described above. Accurate optimization of the supercells for the TM-C composites is very important since each C atom will have different TM as NN with different interatomic bond lengths and interactions. We have used the projector augmented wave method with the Perdew-Burke-Ernzerhof exchange correlation functional (Perdew et al., 1996). Perdew-Burke-Ernzerhof is one of the most popular generalized gradient approximation (GGA) potentials and is reasonably accurate for alloy systems. A relatively high energy cutoff of 500 eV and the stringent electronic convergence criterion of 10 -5 eV were adopted. The force convergence criteria for ionic relaxation was set at −10 -3 eV/Å. We used a single k-point calculation since our models are in the form of large supercells and a single k-point calculation at the zone center is sufficient. This has been successfully demonstrated in many of our recent studies in bio-molecular systems, organic, inorganic and metallic crystals and glasses (Poudel et al., 2015;Adhikari et al., 2016;Poudel et al., 2016;Poudel et al., 2017;San et al., 2018).

Electronic Structure and Interatomic Bonding
The electronic structure calculations use the in-house developed orthogonalized linear combination of atomic orbitals (OLCAO) method (Ching and Rulis, 2012) with the VASP-relaxed structure as input. The all-electron OLCAO method uses the atomic orbitals in the basis expansion which is far more economical than the plane-wave basis and enable us to apply it to complex HEAs using supercells. The merits of the OLCAO method is well documented. In particular, the OLCAO method can provide the effective charge (Q*) or partial charge (PC) on each atom as well as the bond order (BO) values ραβ between any pairs of atoms. They are defined as: In the above equations, S iα . jβ are the overlap integrals between the i th orbital in α th atom in the j th orbital in β th atom. C m jβ are the eigenvector coefficients of the mth occupied band. The PC is the deviation from the neutral charge Q 0 α from the effective charge Q p α on the same atom α. This BO value is affected by the presence of all the nearby atoms which contribute to the BO. The BO, which is basis-dependent and only for short-ranged atomic orbitals, defines the relative strength of the bond. Comparisons of BO calculation using different basis set or methods should thus be treated with caution. The atomic-scale interactions based on DFT calculations are critical for providing the accurate information necessary for their fundamental understanding. More details on the OLCAO method can be found in Ref. (Ching and Rulis, 2012). An important concept critical to the HEAs and its composites with carbon is the TBOD and PBOD. 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 or a supercell model (Dharmawardhana et al., 2014). The TBOD can be resolved into partial components or the PBOD for the different types of atomic pairs (TM-C, TM-TM), or other structural units or compositions in the supercell. The use of TBOD and PBOD as descriptors has special merits. They are based on the quantummechanical calculations and avoid the use of pure geometric parameters in describing their structures and properties. It can be applied to models with different numbers of atoms and compositions, or variations in the structures. For pure HEAs, this new concept is based on understanding the nature of metallic bonding critically related to the theory of formation of HEAs. Although metallic bonding has been extensively discussed in the field of metallic glasses (MGs) (Miller and Liaw, 2007;Pelletier and Qiao, 2019), it has not been thoroughly investigated for HEAs. Metallic bonding is multi-atomic in nature, different from the covalent or ionic bonding. In comparison with the use of the total energy or enthalpy of the system, it is far much simpler and avoid the use of large samplings necessary in the Boltzmann distribution such as in the descriptor called entropy-forming-ability (EFA) (Sarker et al., 2018). Entropy-forming-ability is based on a different theoretical method called automatic flow-partial occupation which utilizes a large number of samples of different configurations.

Mechanical Properties
For the elastic and mechanical properties of the HEAs, we used the stress (σ j ) vs strain (ε j ) response analysis scheme (Nielsen and Martin, 1983;Yao et al., 2007) 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, 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 (Voigt, 1928;Reuss, 1929;Hill, 1952). In the Voight-Reuss-Hill approximation, the elastic constantans are obtained from the average of the two calculations, one based on using maximum stress and the other based on maximum strain. The mechanical parameters K, G, E and Poisson's ratio η are obtained using the standard equations from C ij , Vickers hardness H v is estimated by using the method proposed by Tian et al. (Tian et al., 2012). The G/K ratio also called Pugh ratio (Pugh, 1954) is a useful parameter based on simple empirical arguments on poly crystalline samples , it was claimed that for G/K ratio greater (less) than 0.571, the material is more brittle (ductile).

RESULTS AND DISCUSSION
Before we specifically present and discuss the results for the six groups. We will firstly collectively present the equilibrium structure for each of the 15 models in Table 1 and their elastic and mechanical properties in Table 2 and the PC distribution on each atom in Table 3. Specific details and discussion on the electronic structures and interatomic bonding are presented below for the six groups, starting from m1, in order to facilitate the inter-comparison identify their salient features.
(TiNbTaZrMo) 50 C 50 , m1 Group 1 or model m1 is the main focus of HEA-carbide of the present work since other models are the derivatives of this model. It contains 50% of five TM elements Ti, Nb, Ta, Zr, Mo, and 50% of C. These five TMs in pure HEAs usually form in either FCC or BCC lattices (Todai et al., 2017;Ching et al., 2020). However, each  of them, with very different atomic sizes form binary and ternary carbide with different specific properties. We choose these five refractory TMs as the most representative TM group for the study of HEA carbides (Sarker et al., 2018). As a matter of fact, the same 5-component in the pure HEAs in the BCC structure was the focus of the recently published work on 13 bio-compatible HEAs (Ching et al., 2020). In presenting the calculated results for m1, we also present the results of the 15 models collectively and discuss those results specific for each group later. Figure 1 Shows the ball and stick illustration of model m1 in the FCC rock-slat structure of the 512-atom supercell. Similar illustrative figures of other models are not shown. Table 1  It is interesting to point out that the five refractory TM elements in m1 form pure HEA in BCC lattice [m3 in ref. (Ching et al., 2020)]. The number of TM in HEA-C is 256, very close to the 250 TM in pure HEA m3 of ref. (Ching et al., 2020). The main differences between the model TiNbTaZrMo in BCC and FCC-C lattice are all listed in Supplementary Table S1. It can be seen that the in HEA BCC lattice there is a significant lattice distortion, whereas the FCC HEA-C rock salt lattice, lattice distortion is negligible since all carbon atoms are in the anion   sub-lattice. There is no TM-TM NN boning in the FCC-C model. This fundamental difference in the two models results in drastic difference in their mechanical properties. The FCC-C are significantly stronger than BCC HEA. The Vickers hardness of FCC-C is almost 10 times larger than BCC model. In addition to the comparison between FCC rock salt and BCC, the PC and effective charge of each element in these two models is quite different as shown in Supplementary Table S2.
The Calculated total DOS and atom-resolved partial DOS (PDOS) for m1 is shown in Figure 2. And those for all 15 models are shown in Supplementary Figure S1. As can been seem, the Fermi level at 0.0 eV locates very close to a local minimum in the TDOS, a signature of the stability of these TM carbides. The values of the DOS at the Fermi level or N (EF) and their atomic decomposition are listed in Table 4. For m1, N (EF) 154.059 (States/eV) and their components are respectively 36.221, 20.405, 23.975, 18.290, 12.870, 42.299 (States/eV) for C, Ti, Nb, Ta, Zr, Mo, respectively.
The calculated PC distribution in every atom in m1 is shown in Figure 3. Similar plots for all other models are shown in Supplementary Figure S2. Clearly, the C atoms lost charges to all TM elements. The amount of charge transfer depends on the specific TM and their local environment in the supercell. In m1, Mo receives less electrons from C, Nb and Zr receive the most with Ti and Ta in the middle. We now refer to Table 3 which list the PC, the effective charge Q* and the valence electron count (VEC) of each element in each in each model. In Supplementary Table S2, we compare the effective charges Q* and PC of same 2 models FCC-C-512 and BCC-250 in Supplementary Table S1. As can be seen, the same five TM elements Ti, Nb, Ta, Zr, Mo have very different values and trends although they have same VEC. This clearly demonstrate that VEC should never be used in any theory or explanations for HEAs or their composites.
In Figure 4, we display BO-BL for NN TM-C bonds and SNN TM-TM bonds. Similar plots for all other models are shown in Supplementary Figure S3. It can be seen that for the TM-C bonds, the data are somewhat clustered depending on the TM with the general trend of longer the BL, less the BO. The same trends also shown in the TM-TM they are more scattered. In particular, Mo-Mo bonds tend to have smaller BL and higher BO values than other pairs. Such trend can be seen much easier in the Figure 5, where the data in Figure 4 are further break down into specific atomic pairs.
From the sum of all the BO values, we can obtain the TBOD and the PBOD when normalized by the cell volume. These are shown in Figure 5 for the 15 models in six groups including m1. The highest TBOD is in m8 [0.0614 e/( A) ] and the lowest TBOD is in m5 [0.0503 e/( A) ]. It is obvious that the TM-C bonds contribute far more than the TM-TM bonds to the TBOD. For m1, the contribution to TBOD follow the trend of C- Ta  Apparently, the above listed numbers show that 88% of bonding in HEA-C for m1 comes from very strong TM-C bonding, and only 12% from TM-TM bonding. The fact that carbon bonding has dominated in HEA-C means it will have better and much stronger mechanical properties.
In Table 2, we listed the calculated mechanical parameters for all 15 models in six groups. These consist of the elastic coefficient C 11 , C 12 , and C 14 for a cubic structure, bulk modules K, shear modulus G, and Young's modulus E and Poisson's ratio η. Also listed are the values for the Vickers hardness Hv, TBOD, and energy per atom since they are all intimately related to the mechanical properties. Figure 6 shows the correlation matrix between the mechanical properties among themselves and TBOD in all 15 models. The largest correlation is between Young's and Shear modulus, and the second largest will be Bulk with Young's and Shear moduli. These strong correlations reflect the fact that all models are related or derived from m1 with specific variations. The figure also shows that there is a negative correlation of Poison's ratio with any other properties although Poisson's ratio is not correlated with TBOD. It is vital to make a correlation between physical properties in models with different chemical compositions in order to reveal the subtle differences that could impact meaningful predictions in the synthesis of HEA ceramic composites with carbon.

Reduction of Carbon Contents in m1: (m1A-m1E)
The Group 2 has five models with decreasing contents of C from 256 C atoms in m1 with we designate to be 100% of C. The percentage of C atoms in the five models from m1A to m1E are reduced respectively to 96%, 94%, 90%, 85%, and 80%. In these Group 2 models the number of C atoms are randomly removed for each % of reduction. This certainly affects their structure and Our discussion in this subsection focuses on the trend associated with the steady C reduction and in comparison with experimental data (Conclusion). Briefly, these effects can be summarized as follows: 1) Deducting carbon content causes volume decreasing since there are less atoms in the supercell.
2) The mechanical properties also decrease because of decreasing in the strong TM-C bonds. This is reflected in the reduced Vickers hardness Hv and tends to be more ductile with a decreasing G/K ratio.
3) The energy per atom is getting slightly larger since there are less strong TM-C bonds. 4) It can be observed that the PC distributions are more scattered. From m1A to m1E.
To correlate the change in mechanical properties due to C reduction, we show in Figure 7A, the plot of TBOD and the bulk modulus K against the C percentage for models m1 and m1A-m1E. Within this range, TBOD decreases in a linear fashion as expected since the reduction in C content mirrors the reduction of the strong TM-C bonds. It can also be seen that this trend is closely correlated with K as expected. There is a subtle difference at the 96% C data point where the data point for bulk modulus falls slightly below the straight line. This subtotal difference is related to the specific TM-C bonds missing due to the removal of 4% of C atoms. In Figures 7B,C, we show the similar plots for shear modulus G and Young's modulus E vs. C%, and Poisson's ratio η and G/K vs C% respectively. It shows that both G and E decrease with C% but not in a linear fashion as in Figure 7A but have a slight change in the slope at 96% of C content. Thus, the properties change due to reduction in C could be more complicated and will definitely affect the interpretation of experiments related to C deficiency in HEA-C composites.
Variation in Ti/Mo Ratio From m1: (m1F, m1G) The two TM elements Ti and Mo in m1 always play some special roles in HEAs. To investigate their opposite influence on the structures and properties of TM carbides, the Group 3 compares the two models (m1F, m1G) with composition formula of Ti 5 (NbTaZr) 10 Mo 15 C 50 and Ti 15 (NbTaZr) 10 Mo 5 C 50 with different Ti/Mo ratios of 1-3 and 3-1 respectively. m1F has a larger cell volume than m1G because Mo is a much larger atom than Ti. This is partly reflected into the larger TBOD of m1G over m1F as well as their difference mechanical properties and Poisson's ratio. From the PC table for all 15 models in Table 3, it can be seen that both Ti and Mo have their effective charges bellow the VEC of 4 and 6 respectively, which also make the PC in C to be less effected by different Ti/Mo ratios in m1F and m1G. On the other hand, the N (EF) values listed in Table 4 are quite different, 173.36 States/eV for m1F and 128.044 States/eV for m1G because of significantly different contribution from Ti and Mo atoms at the Fermi level. It suggests that m1F is more stable than m1G. This is quite obvious by looking into Supplementary Figure S1 that the PDOS of Ti and Mo in m1G has a deep local minimum at the Fermi level. We also observe from Supplementary Figure S4 on the BO vs BL plots  This group is similar to m1, but with different ratio of Ti/Mo components. The key difference between group 1 and group 3 is that the mechanical properties of G3 is lower than G1. Moreover, changing the Ti/Mo ratio from 1:3 (Ti 5 /Mo 15 ), to 1:1(Ti 10 /Mo 10 ), and to 3:1(Ti 15 /Mo 5 ), respectively, makes the TBOD slightly higher. Thus, we can see that increasing (decreasing) number of atoms of Ti in systems, gives higher (lower) TBOD. Not only Ti is essential for the system but increasing the number of atoms of Mo is also important with a constant number of Ti (See A Six-Components High Entropy Alloy-Carbide m6: (TiNbTaZrVMo) 50 C 50 ).

Comparison of Five Models With Five Transition Metals by Replacing Mo in m1
With V, Cr, W, Hf: (m1-m5) The next group, Group 4, for us to discuss is to compare the five models with five TM components including m1 to m5 specifically identified as group 1 for in-depth study. m2-m5 replaces Mo in m1 with V, Cr, W, Hf respectively. This enables us to identify the subtle difference of changing only one TM element in TM carbides. In the nine TM involved in the 15 models for HEAs-C, Cr in m3 is the only non-refectory TM element. We observe that the PC of Cr in m3 is both gaining and losing electrons, and one C atom actually gaining an electron. Another observation is that replacing Mo with Hf makes m5 much weaker since it also does not contribute much in the bonding with C. This reveals that Mo is essential for HEAs. Also, in all 15 models Nb does not show the contribute of bonding with itself and with other elements.
We also noted that Zhang et al. (2019) also investigated the electronic structure and mechanical properties of a model similar to our m5 but with a much smaller elongated supercell of 40 atoms using a different method. Their results show a much higher total energy per atom, lower Young's modulus and Shear modulus, and lower G/K ratio. Our results are much accurate and substantial and reflect the importance of using large supercell in the modeling of HEA-C composites.
A Six-Components High Entropy Alloy-Carbide m6: (TiNbTaZrVMo) 50 C 50 The group 5 has only one model m6, but it is the only one with six TM components and hence a larger intrinsic entropy. The difference between this model with m1 and m2 is adding one more element, V, in m1 and replacing Mo with V in m2. Since V is a light atom, the change in properties will not be significant. For instance, the comparison of m6 with m1 shows m6 has higher TBOD. M6 also possesses lower mechanical properties than m1. On the other hand, in comparing m2 without Mo to m6 with Mo, m2 has a larger volume and lower TBOD and mechanical properties. In addition, m2 is more ductile than m6. As has been noted, Mo tends to strengthen the mechanical properties for the HEAs. This effect can be clearly seen also in pure HEAs system in BCC lattice for models m2, m3, and m5 in Ref. (Ching et al., 2020).
Two Models by Doubling the Percentage of Zr and Ta Over the Other Four Transition Metals in m1: (m7, m8) The group 6 is the last group with 2 models m7 and m8. These two models have double percentage of Zr and Ta over the reaming four TMs again explore the subtle difference due to compositional changed concentrating on Zr and Ta, both are present in m1. It can be observed that m7 has much larger volume than m8 that related to the atomic size of Zr and Ta. This is somewhat counter intuitive since Ta atom is at the lower row in the Periodic Table with more electrons. This simply illustrate that in complex HEA-C composites, inter atomic interaction and local bonding controls the structure and properties, particularly the TBOD. m7 has the smallest TBOD and m8 has the largest TBOD among all 15 models signifying the contrasting difference between Zr and Ta. This contrasting difference between m7 and m8 can also be noticed in their PC distributions of Zr in m7 and Ta in m8 as shown in Table 3. Their deviations are much larger than the same elements in the other 13 models. The DOS value at the Fermi level N(E F ) between these two models are naturally to be very different due to different electronic structure (140.05 States/eV in m7 and 183.70 States/eV in m8. Since Ta is at a higher row in the Periodic Table with 5 days electrons than Zr with 4 days electrons, it seems to bond more strongly with C than Zr-C bond. Ta-C bonds contribute 42.57% to the total bonding in m8, whereas m7 has double number of Zr-C bonds atoms contribute only 28.74% to the total bond in m7 even though the number of Zr atoms in m7 and Ta atoms m8 are the same, both twice that of the other TMs. Frontiers in Materials | www.frontiersin.org November 2020 | Volume 7 | Article 575262

CONCLUSION
We have presented a detailed study of 15 models of TM HEA carbides using large supercells of 512 atoms each. The results are presented and discussed by dividing them into six groups for specific comparison and analysis. Detailed electronic structure, interatomic bonding, PC distributions, provide the fundamental understanding of these relatively new HEA-ceramic carbides and their mechanical properties. This unprecedented large database enables us to delineate the subtle differences in their fundamental properties that have not been revealed before or have been long ignored. Specifically, we have reached the following solid conclusions. The method and approach we developed in this paper can be extended to other more complex HEA-ceramics including oxides, nitrides, borides or silicides. The opportunity is unlimited.
(1) The TM carbides with 50/50 composition provide an ideal structure for HEA ceramics in which the NN consist of TM-C bonds and the SNN consist of TM-TM pairs. There are no C-C bonds.
(2) The rock-salt structure mimics the ideal TM-C binary alloys. However, we now have the additional advantages of having HEAs in the TM distribution, greatly expand the composition space and variety. We are able to analyse the results at a deeper level quantitatively especially in the BO vs BL plots than just using the vague argument of atomic size misfit prevailed in many published literatures. (3) We utilize the novel concept of TBOD and PBOD to characterize the complex alloy systems in HEAcarbides which is more straightforward and easier to understand. (4) Specific modelling for C-deficiency provide another possibility of extending HEA-ceramics with different C contents, further expand the composition space. (5) The design of different composition of the TM elements in this work reveal many of the subtle difference due to atomic origin. These include the special insights obtained that Ta and Ti as well as Mo components make a large contribution in bonding with carbon, which enhance their strength. (6) Our ab initio calculation in this work demonstrate that the routine use of the valence electron count for TM elements as key parameters in characterizing HEA components in untenable. (7) The methods developed and demonstrated in this paper can be easily extended into other HEA-ceramics such as nitrides or borides. The insights gained can facilitate the design new radiation resistant HEA ceramic composites at high temperature. Frontiers in Materials | www.frontiersin.org November 2020 | Volume 7 | Article 575262