The coupling of carbon non-stoichiometry and short-range order in governing mechanical properties of high-entropy ceramics

Introduction

Ultra-high temperature ceramics (UHTCs), typically defined as compounds with melting points exceeding 3000 °C, predominantly consist of carbides, nitrides, and borides derived from early transition metal elements in groups IVB and VB^1,2. Notably, the unique combination of covalent, ionic, and metallic characteristics in refractory transition metal carbides (TMCs) enables these ceramics to withstand extremely high temperatures, intense heat fluxes, elevated radiation levels, and significant mechanical stresses^3,4,5. However, as the demand for enhanced material performance grows, traditional binary carbide ceramics are increasingly inadequate for applications in extreme conditions due to their inherent brittleness, low toughness, and poor thermal shock resistance. Hence, developing advanced UHTCs is of urgent necessity, and a common approach is to mix cation elements to form carbide solid solutions^6,7,8,9.

Inspired by the entropy-stabilization effect, there has been significant advancement in high-entropy transition metal carbide ceramics (HECCs) to address the limitations of binary carbides. HECCs typically consist of four or more transition metal cations, often in equimolar or near-equimolar ratios, coexisting within the B1 crystal sublattice. Studies have demonstrated that these multicomponent transition metal carbides (MTMCs) exhibit superior properties, including high hardness^10,11,12,13, low thermal conductivity^14,15,16,17, outstanding oxidation¹⁸, irradiation resistance^19,20, and favorable high-temperature flexural behavior²¹. These characteristics position HECCs as promising candidates for emerging ultra-high temperature applications.

It is well established that carbides are susceptible to carbon deficiency due to the low formation energies of carbon vacancies. This vulnerability also extends to MTMCs, which typically contain carbon vacancies and can maintain stability across a wide range of carbon stoichiometries²². In contrast to binary TMCs, the formation of carbon non-stoichiometry in MTMCs with multiple transition metals depends on the local environments surrounding the carbon atoms and exhibits site preference. Additionally, the incorporation of multiple principal transition metals can lead to the formation of short-range order (SRO), a prominent characteristic in high-entropy systems that has garnered considerable interest, especially in high-entropy alloys^23,24. The presence of SRO can significantly alter the arrangement of elements, resulting in favorable and unfavorable local atomic environments around each site. However, the formation mechanism of SRO in HECCs and its effects have yet to be elucidated.

Both carbon non-stoichiometry and SRO induce atomic-level heterogeneity and profoundly influence the mechanical properties of MTMCs. For MTMCs containing Group IVB (Ti, Zr, Hf) and VB elements (V, Nb, Ta), experiments have demonstrated that carbon deficiency leads to lattice softening, which is crucial for promoting local plasticity, as evidenced by the reductions in elastic modulus and hardness^25,26,27. First-principles calculations for MTMCs such as (TiZrHfNbTa)C²⁸ have indicated that carbon deficiency-caused changes in hardness are primarily driven by variations in covalent bond density. Specifically, the metallic states exposed after carbon vacancy formation cannot compensate for the loss of covalent bonds, resulting in decreased nanohardness. In terms of SRO, findings in high-entropy alloys suggest that SRO can generally enhance mechanical strength; however, whether similar effects occur in HECCs remains unknown. Revealing the combined effects of carbon non-stoichiometry and SRO is an essential step toward a fundamental understanding of the structure-property relationship of HECCs.

The vast composition space and diverse local atomic environments (LAEs) associated with MTMCs make machine learning (ML) an ideal tool for investigating how local variations in atomic arrangments, including carbon stoichiometry and SRO, impact their properties^29,30,31,32. Nonetheless, the application of ML in HECCs is limited, with most previous models primarily focusing on macroscopic properties, such as single-phase formation ability. Recently, a random forest model was trained to predict the local environment dependence of the carbon vacancy formation energies in (HfNbTaZr)C³³, utilizing a simplified description of the local environment that emphasizes local chemistry. However, the effects of local geometry, especially the local relaxation effects and local lattice distortion, are not taken into account.

In this study, we investigate the coupling effects of carbon non-stoichiometry and SRO on the mechanical properties of HECCs by integrating density functional theory (DFT), Monte Carlo (MC) algorithm, and ML models. The (TiZrHfNb)C^34,35 system is selected because it has been experimentally demonstrated to be a single phase and exhibit superior mechanical strength. Our results reveal that Ti and Nb elements promote local carbon non-stoichiometry and lead to a local strengthening effect. Furthermore, these two elements are also the primary contributors to the formation of SRO. By capturing both the local chemical and geometric environment, we develop ML models capable of accurately predicting carbon vacancy formation and related mechanical properties. It is expected that this framework, which synergistically integrates DFT, MC, and ML, can effectively accelerate the exploration of HECCs and inspire the optimization of their properties.

Results

LAE dependence of carbon non-stoichiometry

Carbon deficiency is influenced by specific LAEs. Figure 1 illustrates the relationship between carbon vacancy formation energies and the number of first nearest neighbor (1NN) metallic elements. In general, regardless of whether the system consists of two, three, or four metallic elements, their influences on reducing the formation energy follows this order: Nb > Ti > Zr > Hf. For instance, when comparing Ti with Zr/Hf elements, vacancies are preferentially associated with Ti. Consequently, in (TiZr)C and (TiHf)C, vacancy formation energies decrease as the number of Ti atoms surrounding the vacancy increases. However, compared to Nb, vacancies do not exhibit a preference for Ti over Nb, resulting in an increase in vacancy formation energies as the number of Ti atoms increases in (TiNb)C.

The coupling of carbon non-stoichiometry and short-range order in governing mechanical properties of high-entropy ceramics — **Fig. 1: LAE dependence of carbon vacancy formation.**

The average and distribution range of vacancy formation energies are presented in Fig. 2a and compared to those of binary TMCs. In binary carbides, carbon vacancies are more likely to form in NbC, while they are less favored in TiC, ZrC, and HfC. This trend also extends to MTMCs, in which MTMCs containing Nb exhibit lower formation energies. Nonetheless, the average vacancy formation energies of MTMCs (depicted as blue triangles) are noticeably lower than the rule-of-mixture values of binary carbides (depicted as red triangles), except for (HfNb)C. The distribution range of formation energies in MTMCs is also influenced by the variations in their constituent elements. For example, the effects of Hf and Nb on vacancies differ significantly, resulting in a wider distribution of vacancy formation energies in (HfNb)C. In contrast, the effects of Zr and Hf on vacancies are similar, leading to a narrower distribution in (ZrHf)C.

**Fig. 2: Analysis of the vacancy formation mechanism.**

To highlight the critical role of local relaxation induced by carbon vacancies, we calculate the relative atomic displacement and local lattice strain for the 1NN (metallic atoms) and 2NN shells (carbon atoms):

$${varepsilon }_{d}=frac{{d}_{{rm{vac}}}-{d}_{{rm{perfect}}}}{{d}_{{rm{perfect}}}}times 100( % )$$

(1)

$${varepsilon }_{{rm{local}}}=frac{1}{m}mathop{sum }limits_{i=1}^{m}frac{{d}_{{rm{vac}}}-{d}_{{rm{perfect}}}}{{d}_{{rm{perfect}}}}times 100( % )$$

(2)

where m is the total number of atoms in the corresponding NN shell (m = 6 for 1NN and m = 12 for 2NN), d_perfect is the distance between an atom in the NN shell and the original carbon site in a perfect crystal, and d_vac is the distance between the same atom and the vacant carbon site after vacancy formation, assuming the carbon site remains in its original position. The distributions of atomic displacements for 1NN metallic atoms and 2NN carbon atoms are shown in Fig. 2c and Fig. S3a, respectively. The positive atomic displacements of 1NN metallic atoms indicate that they move away from the vacant site, while the negative atomic displacements of 2NN carbon atoms suggest a shift toward the vacancy. Notably, Fig. 2c reveals that Ti and Nb atoms experience more notable displacements, leading to greater local lattice strain that facilitates the formation of carbon vacancies. Conversely, Zr and Hf exhibit the opposite effect. Figure 2d further demonstrates that the vacancy formation energies decrease with increasing local lattice strain, which corresponds to a higher degree of local relaxation.

Figure 2e illustrates the charge density difference with and without carbon deficiency. The carbon vacancy induces the most significant changes in charge distribution around the metallic atoms in the 1NN shell, while the remaining carbon atoms at 2NN show minimal variation. Particularly, the charge that was originally located between the 1NN metallic atoms and the central carbon atom shifts towards the metallic atoms, away from the carbon vacancy. This movement strengthens the covalent bond between the 1NN metallic atoms and the outer carbon atoms within 2NN and 4NN. Moreover, there is a greater charge transfer for Ti/Nb compared to Zr/Hf, which aligns with the observation that carbon vacancies favor Ti/Nb but do not favor Zr/Hf. This combined analysis of local atomic relaxation and charge transfer explains the LAE dependence of carbon vacancy formation energies as observed in Fig. 1.

While 1NN metallic atoms exert the strongest influence on carbon vacancy formation energies, contributions from 3NN and beyond cannot be ignored. The dependence of 3NN and 5NN metallic elements on carbon vacancy formation energies is shown in Fig. S1 and S2. Interestingly, the influence of 3NN metallic elements exhibits an opposite trend to that of 1NN: Ti/Nb generally contributes to an increase in vacancy formation energies, while Zr/Hf tends to reduce them. Nevertheless, the overall influence of 3NN becomes less pronounced, and when extended to 5NN, it becomes negligible, showing no discernible correlation with vacancy formation energies.

To elucidate the contrasting behavior observed in 3NN, we analyze the atomic displacement and charge transfer patterns, as shown in Fig. S3. While 1NN metallic atoms move away from the vacant carbon site, 3NN atoms shift toward it. This opposing displacement is accompanied by a reverse charge transfer: 1NN atoms redistribute charge toward neighboring metallic atoms, away from the vacancy, whereas 3NN transfer charge toward the vacancy.

The distinct responses of 1NN and 3NN atoms to carbon vacancies account for their differing influences on vacancy formation energies. When a carbon vacancy forms, the disruption of metal-carbon bonds causes the 1NN metallic atoms to move outward. These metallic atoms redistribute charge and form bonds with outer carbon atoms, effectively accommodating the vacancy. In contrast, 3NN metallic atoms shift inward, transferring charge toward the vacancy. However, this inward displacement and charge transfer are insufficient to form new bonds with the original carbon site, making 3NN atoms less effective in adapting to the vacancy. Despite the distinct responses of 1NN and 3NN to the carbon vacancy, the magnitude of atomic displacement and charge transfer for 3NN is significantly smaller—approximately 15% of that observed for 1NN atoms. Consequently, the influence of 3NN atoms on vacancy formation energies is minimal, with 1NN metallic atoms remaining the primary contributors.

Influence of carbon non-stoichiometry on mechanical properties

The elastic constants of considered MTMCs with (96.9% carbon) and without (100% carbon) carbon non-stoichiometry are computed from DFT calculations, which allow for the derivation of the bulk modulus B, shear modulus G, elastic modulus E, and Poisson’s ratio υ. The overall changes in B, G, and E before and after the formation of carbon non-stoichiometry are plotted in Fig. 3. The presence of carbon vacancies results in a reduction in the modulus of MTMCs in general, with the observed trends showing consistent behavior. The distribution widths of the moduli are similar to those of vacancy formation energies, primarily influenced by the metallic elements constituting the system. Fig. S1 illustrates the trend between E, B, υ and the number of local metallic atoms surrounding the carbon vacancy. It is found that Ti and Zr in the local environment have a relatively minor impact on mechanical properties. In contrast, Hf and Nb exhibit a more pronounced influence, especially on the Poisson’s ratio. Consequently, Hf leads to a decrease in elastic modulus, while Nb causes an increase.

To better understand the changes in modulus caused by carbon non-stoichiometry, the projected density of states (PDOS) for local metallic atoms (d orbitals) and carbon atoms (p orbitals) are provided in Fig. 4. The introduction of carbon vacancies results in a decrease in the PDOS peak within the energy range of −4 eV to −2 eV (red circles), suggesting diminished hybridization between metal and carbon atoms. This reduction in hybridization accounts for the reduced modulus observed in non-stoichiometric systems. The PDOS of metallic atoms displays a typical pseudogap (blue circles), characterized by a valley between distinct peaks representing bonding and anti-bonding states. Notably, after the formation of vacancies, the PDOS in this valley increases, giving rise to a minor peak that indicates charge transfer.

**Fig. 4: Projected density of states (PDOS) analysis.**

The origin of this peak can be traced to the partial electronic charge density in the energy range from −1 eV to the Fermi level (E_F, set to 0 eV) before and after the formation of a carbon vacancy, as presented in Fig. 5. It can be found that at the vacant carbon site, the original carbon-metal bonds are replaced by notable d-d covalent bonds. This type of d-d bonding is particularly pronounced between Ti and Nb, which explains why the presence of these elements in the local environment mitigates the decline in elastic modulus.

**Fig. 5: Electronic structure analysis.**

Short-range order formation and its effects

The mixing properties of different transition metallic pairs in the B1 lattice of HECCs exhibit variability, leading to competing interactions that result in SRO formed by favored and disfavored atomic pairs. In HECCs, it is noteworthy that the metallic pairs are located at the 2NN positions, which may significantly expand the SRO domain compared to high-entropy alloys, where 1NN interactions are primarily involved. This elemental heterogeneity depends on the interplay of processing conditions and the thermodynamic affinities among different metallic elements. As a result, the arrangement of metallic elements is not random but instead exhibits specific patterns. To investigate the SRO formation in (TiZrHfNb)C, which has not been previously explored, we conduct Metropolis Monte Carlo (MMC) simulations using mixing energies obtained from DFT.

Table 1 presents DFT calculated mixing enthalpies between pairs of different metallic elements in (TiZrHfNb)C. Based on these energies, Fig. 6a, b illustrates the results of the MMC simulations, showing the evolution of the total mixing enthalpy and SRO parameters at 500 K over 3000 steps in a supercell with 512 atoms. After 1500 steps, both the total mixing enthalpy and SRO parameters reach a stable state. A noticeable presence of SRO is observed within the system. Ti-Ti and Zr-Zr pairs exhibit the most negative SRO parameters, indicating that Ti and Zr have a strong tendency to cluster with themselves. Additionally, Ti and Nb atoms show notable clustering tendencies, while Ti and Zr/Hf atoms tend to separate. This behavior is consistent with the most negative mixing enthalpies of Ti-Nb pairs and the positive mixing enthalpies of Ti-Zr/Hf pairs. Furthermore, the results align with the greatest charge transfer observed between the Ti-Nb pairs, as discussed previously.

Table 1 Mixing enthalpies (in meV/atom) between pairs of metallic elements in (TiZrHfNb)C

Full size table

**Fig. 6: Short-range order evolution.**

To quantify the degree of SRO, we calculate the sum of absolute values of SRO parameters for all atomic pairs:

$${rm{SRO}},{rm{degree}}=mathop{sum }limits_{i=1}^{n}mathop{sum }limits_{j=i}^{n}|{alpha }_{ij}|$$

(3)

where n = 4 is the number of metallic elements, α_ij is the SRO parameter between metallic pair consisting of i and j. As shown in Fig. 6c, the degree of SRO increases rapidly within the range of 0-600 steps, followed by a slower ascent and stabilization after 1500 steps.

The impact of SRO on the mechanical properties of the considered MTMC is further studied through DFT calculations. For this purpose, SRO is created in a supercell containing 64 atoms following the MMC procedure. Structures at 0, 300, 600, 900, 1200, and 1500 steps are selected for calculating elastic constants, with the results presented in Fig. 6c. It can be observed that as the degree of SRO increases, the elastic modulus of (TiZrHfNb)C gradually rises. Therefore, a certain degree of SRO is beneficial for improving the mechanical strength of HECCs, consistent with the findings in high-entropy alloys³⁶.

Machine learning

Both carbon non-stoichiometry and SRO alter the local atomic arrangement and significantly influence the mechanical properties of MTMCs. While carbon stoichiometry may soften MTMCs, SRO can enhance their strength. To elucidate these two competing effects, we develop ML models based on the established DFT database spanning from binary to quinary MTMCs to predict the LAE dependence of mechanical properties for HECCs. We have selected five categories of features:

(1) Compositions. Given that the dataset consists of several different MTMCs, one-hot encoding is employed to represent the composition of each system.

(2) Nearest neighbors. For each carbon atom, the 2NN and 4NN positions are occupied by carbon atoms. Hence, the local environment depends only on the arrangement of metallic atoms in the 1NN, 3NN, and 5NN shells. We count the number of metallic elements within the 1NN, 3NN, and 5NN shells as input features.

(3) Mean and deviation of local properties. We calculate the mean and standard deviation of local properties in 1NN, 3NN, and 5NN shells as follows:

$$overline{{p}^{k}}=mathop{sum }limits_{i=1}^{N}{x}_{i}^{k}{p}_{i}^{k}$$

(4)

$${sigma }_{p}=sqrt{mathop{sum }limits_{i=1}^{N}{x}_{i}^{k}{[{p}_{i}^{k}-overline{{p}^{k}}]}^{2}}$$

(5)

where ({p}_{i}^{k}) denotes either the properties of the element i or the properties of binary carbides constituted from element i in the k^th neighbor shell, N is the total number of species, and x_i is the concentration of species i in the k^th neighbor shell. For predicting vacancy formation energies, the properties p considered encompass the atomic number, period number, group number, atomic weight, electronegativity, covalent radius of the element, along with the lattice constant, melting point, valence electron concentration, and vacancy formation energy of the corresponding binary carbide. For the prediction of mechanical properties, in addition to the aforementioned properties, shear modulus, bulk modulus, elastic modulus, and hardness of each binary carbide are also considered.

(4) AGNI fingerprints³⁷. AGNI values are used to capture the local atomic arrangement within a 5 Å radius of the vacant carbon site. Specifically, we utilize the initial unrelaxed structures as input to ensure that the ML model can make predictions without requiring additional DFT calculations. For example, when considering Ti, we exclude other metallic elements and compute the AGNI values solely for Ti atoms. This approach ensures that the AGNI fingerprints represent the geometric configuration formed by Ti atoms surrounding the vacancy site, effectively reflecting the local atomic arrangement of Ti near the vacancy.

(5) Gaussian symmetry functions G² and G⁴ ³⁸. Gaussian symmetry functions are calculated for each species within a 5 Å radius of the vacant carbon site. Similar to AGNI fingerprints, the initial unrelaxed structures are used as input to capture the local atomic arrangement. Different from AGNI, Gaussian symmetry functions use radial (G²) and angular (G⁴) terms to encode atomic positions and angles within the local environment, providing a comprehensive description of both pairwise distances and angular distributions.

In this set of five feature categories, the first describes the overall composition, while the second and third primarily capture the local chemical information. AGNI and Gaussian symmetry functions, on the other hand, capture the local geometric features. These features are straightforward to obtain and do not require additional calculations, allowing for deployment without extra computational overhead.

The accuracy of predicted carbon vacancy formation energies, bulk modulus, shear modulus, and elastic modulus by different ML models is given in Fig. 7a–d. The determination coefficients R² on the test set are slightly lower than those on the training set, suggesting that the models are not overfitted. Besides, it can be observed that even simple linear models, such as logistic regression (LR), least absolute shrinkage and selection operator (LASSO), and ridge regression (RR), demonstrate good performance. This suggests that the selected features effectively capture the local environment, leading to highly effective model training. For predictions, we adopt the models with the best performance. Specifically, kernel ridge regression (KRR), random forest (RF), support vector machine (SVM), and KRR are chosen for the predictions of vacancy formation energies, bulk modulus, shear modulus, and elastic modulus, respectively. The predicted properties from these four models are depicted in Fig. 7e–h. Notably, some of the mechanical moduli are discretely distributed, a phenomenon attributed to the significant differences in modulus across various MTMC systems, resulting in a segmented distribution of modulus values. Overall, the ML models yield R² scores exceeding 95%, demonstrating exceptional model performance.

To gain a deeper understanding of the model performance, it is crucial to clarify the underlying physical reasons behind why some models outperform others. For vacancy formation energy predictions, linear models such as LR, LASSO, and RR perform well, suggesting a strong linear relationship between vacancy formation energies and input features. In contrast, ensemble methods like RF, gradient boosting regression (GBR), and extreme gradient boosting (XGB) show signs of overfitting, likely due to their complexity. KRR achieves the best performance by effectively capturing subtle non-linear patterns through kernel methods while avoiding excessive model complexity. Its ability to map data into higher-dimensional spaces makes it particularly suitable for modeling the intricate relationships between LAEs and vacancy formation energies.

For modulus predictions, ensemble models such as RF, GBR, and XGB outperform linear models, indicating that the relationship between input features and modulus is more complex and non-linear. Ensemble models excel at capturing higher-order interactions and non-linear dependencies, resulting in improved accuracy. Models like SVM and KRR also perform well, as they are capable of handling non-linear relationships through their kernel functions. SVM efficiently captures patterns in high-dimensional spaces, while KRR extends ridge regression by incorporating kernel methods, enabling it to model more complex relationships.

With the well-trained ML models, we further analyze the importance of different features on prediction accuracy. To explore the roles of the five categories of features, each category is individually used to train ML models; the results are shown in Fig. 8a, b. In predicting vacancy formation, it is evident that the nearest neighbors and local properties contribute the most to the accuracy of most models. In contrast, the overall compositions, AGNI fingerprints, and Gaussian symmetry functions contribute less. Therefore, the importance of local chemistry is greater than that of local atomic arrangement for predicting carbon vacancy formation.

**Fig. 8: Feature importance analysis.**

In predicting elastic modulus, apart from the nearest neighbors and local properties, the overall composition also plays an important role. The influence of different nearest neighbor shells on the model performance is further analyzed, as shown in Fig. 8c–f. As the nearest neighbor shells expand from 1NN to 5NN, the determination coefficient R² of vacancy formation energy and elastic modulus prediction results increases while the mean absolute error (MAE) decreases, indicating improved model performance.

Finally, focusing on the two feature categories that contribute most to the model — nearest neighbors and local properties — we separately train models using each category and analyze the feature importance using the Shapley additive explanation (SHAP). The SHAP values are illustrated in Fig. 8g–j. Among the nearest neighbor atoms, Nb plays the most crucial role. An increase of Nb in the nearest neighbor leads to a decrease in vacancy formation energy and an increase in elastic modulus, consistent with earlier DFT results. Regarding local properties, the mean melting point of binary carbides and the standard deviation of vacancy formation energy rank first in predicting vacancy formation energy and elastic modulus, respectively. The melting point, carbon vacancy formation energy, and elastic modulus together reflect the intrinsic physical properties of carbon-metal (C-TM) bonds. The relationship between SHAP and feature values indicates that higher melting points correspond to higher vacancy formation energies and, thus, higher moduli in the system. The analysis provides a solid explanation for the physical significance of our ML models.

Coupling effects of SRO and carbon non-stoichiometry

The emergence of SRO leads to the rearrangement of elements and alters carbon vacancy formation, creating a coupling between these two processes. Given that their impacts on mechanical properties are contrary, this coupling exerts a complex influence on the mechanical response of HECCs. Based on the trained ML model, we can predict vacancy formation and mechanical properties across different SRO states. Specifically, for each SRO structure, vacancy formation energies and elastic moduli are predicted for each possible vacant site, assuming an equal probability of vacancy formation at any carbon atomic position. The predicted results are shown in Fig. 9a, b. It is observed that as the SRO degree increases, there is no significant change in the average value of carbon vacancy formation energies. Additionally, the elastic moduli are minimally affected by the coexistence of SRO and carbon non-stoichiometry in this context.

**Fig. 9: The coupling effects of carbon non-stoichiometry and SRO on mechanical properties.**

However, in practice, only carbon atomic sites with lower vacancy formation energies are likely to form vacancies, suggesting that preferential carbon vacancy formation should be considered. For this reason, we focus on vacant sites with vacancy formation energies below the average, as these sites are more prone to develop vacancies. Figure 9c illustrates the evolution of elastic modulus during SRO development when accounting for such site preferences. The results indicate that SRO enhances the modulus of (TiZrHfNb)C, even when carbon vacancies are present. This enhancement can be attributed to the clustering of Ti and Nb atoms due to their negative mixing enthalpy, as depicted in Fig. S3. This favorable SRO clustering increases the likelihood of encountering Ti and Nb atoms in the vicinity of carbon vacancies. According to our previous analysis, the presence of Ti and Nb elements in the local environment contributes to the enhancement of elastic modulus. Hence, an increase in elastic modulus is observed when preferential carbon vacancy formation is taken into account.

In Fig. 9b, c, we also display the average elastic modulus and its distribution range for pristine (TiZrHfNb)C. The elastic modulus of structures with carbon vacancies generally falls below the average value (dashed line) of perfect structures despite the strengthening effects of SRO. Therefore, while SRO contributes to the enhancement of mechanical strength, it cannot fully counterbalance the softening effects caused by carbon non-stoichiometry.

Coupling effects of SRO and varying carbon content

The above results based on a single carbon vacancy indicate that carbon vacancies tend to soften HECCs, even when potential chemical SRO is considered. To further investigate the effects of decreasing carbon content, we calculate the elastic modulus of (TiZrHfNb)C_x for values of x ranging from 1 to 0.8125 using DFT calculations. To account for different arrangements of carbon vacancies, ten SQS structures are constructed for each carbon content. The resulting average elastic moduli are presented in Fig. 10a. As x decreases from 1 to 0.8125, the elastic modulus decreases from 431 GPa to 344 GPa.

**Fig. 10: Coupling effects of SRO and varying carbon concentration.**

To explain the mechanism behind the decline in elastic modulus, we show the total density of states (TDOS) at varying carbon vacancy concentrations (Fig. 10b). As the carbon content decreases, the primary peak corresponding to C-p and M-d hybridization in the range of −4 eV to −2 eV gradually diminishes, which accounts for the reduction in mechanical properties. Concurrently, the minor peak in the range of −1 eV to 0 eV, representing d-d bonds between metals, gradually increases, consistent with the trend observed for a single carbon vacancy.

Intriguingly, as the concentration of carbon vacancies rises from 0 to 20%, the shape of the TDOS remains largely unchanged compared to that of a single vacancy, while the elastic modulus decreases linearly as a function of carbon content. This result suggests that interactions between carbon vacancies have minimal impact on the mechanical properties of HECCs when the carbon vacancy concentration is below 20%. Consequently, the variation in the mechanical response can be analyzed in terms of single vacancies, and our ML model can, therefore, predict the elastic modulus of HECCs across different carbon concentrations. Specifically, assuming that the reduction in elastic modulus caused by each carbon vacancy remains consistent, the decline for structures with two or more carbon vacancies can be predicted using the percentage reduction in modulus induced by monovacancies and calculating their cumulative softening effect. Figure 10c shows the ML-obtained coupling effects of SRO and varying carbon concentrations. Generally, regardless of whether the HECCs are carbon-deficient, the presence of SRO enhances the elastic modulus. However, this enhancement is insufficient to compensate for the reduction caused by carbon deficiency. These findings underscore the dominant role of carbon non-stoichiometry in softening the mechanical strength of MTMCs, despite the possible strengthening effects of SRO.

Discussion

TMCs composed of metallic elements from the same group usually exhibit similar properties, such as comparable hardness³⁹, stacking fault energies⁴⁰, and fracture toughness⁴¹ due to their similar electronic structures. However, in our study of (TiZrHfNb)C, we find that Ti (group IVB) and Nb (group VB) respond similarly to carbon vacancy formation, despite belonging to different groups.

To further investigate this phenomenon, we compare the carbon vacancy formation energies of binary TMCs from groups IVB and VB, as reported in the literature (Table S2). Generally, the carbon vacancy formation energies follow the trend: VC < NbC < TaC < TiC < ZrC < HfC, with TMCs from group IVB exhibiting lower vacancy formation energies compared to those from group VB. Notably, this trend remains consistent in HECCs. For example, in (ZrHfNbTa)C^33,42, the order of increasing vacancy formation energies is Nb < Ta < Zr < Hf. Similarly, in (TiZrHfNb)C, the carbon vacancy formation energies follow the sequence Nb < Ti < Zr < Hf, aligning with the trends in binary TMCs. This consistency indicates that the influence of high-entropy mixing on vacancy formation energies is approximately linear, as evidenced by the good performance of linear ML models in capturing these trends. Thus, the similarity between Ti and Nb in (TiZrHfNb)C is not a result of group similarity but rather a consequence of their relatively low carbon vacancy formation energies compared to Zr and Hf.

To explain the underlying reason behind the observed trends in vacancy formation energies, we calculate the local electronegativity difference between carbon and neighboring metallic elements around a vacancy for (TiZrHfNb)C:

$$Delta {chi }_{{rm{local}}}=mathop{sum }limits_{i=1}^{4}{c}_{i}|{chi }_{i}-{chi }_{{rm{C}}}|$$

(6)

where c_i is the concentration of the i^th element in the 1NN shell, and χ_i and χ_C denote the electronegativity of the i^th element and carbon, respectively. The electronegativity values for Ti, Zr, Hf, Nb, and C are 1.54, 1.33, 1.30, 1.60, and 2.55, respectively. As shown in Fig. S9a, vacancy formation energies generally increase with increasing local electronegativity difference in (TiZrHfNb)C.

This correlation can be attributed to the role of the electronegativity difference between metal and carbon in determining the nature of metal-carbon bonding. Metals with lower electronegativity exhibit larger electronegativity differences with carbon, promoting greater charge transfer. This enhanced bond stability makes vacancy formation less favorable, resulting in higher formation energies. In contrast, metals with higher electronegativity have smaller electronegativity differences with carbon, leading to reduced charge transfer, which facilitates vacancy formation and lower formation energies.

A similar trend is observed in the reduction of elastic moduli, as shown in Fig. S9b. A larger electronegativity difference between metal and carbon enhances charge transfer, leading to more stable metal-carbon bonds. The disruption of these stable bonds results in a greater decrease in elastic modulus. In contrast, a smaller electronegativity difference results in less stable bonds, making the modulus less sensitive to vacancy formation.

In binary TMCs, it is well established that carbon non-stoichiometry significantly influences their performances. Specifically, the carbon content, or carbon-to-metal (C/M) ratio, affects microstructure, densification, phase evolution, and mechanical properties. An optimal C/M ratio is typically identified for achieving the best mechanical strength; deviations from this optimal ratio can lead to undesirable microstructural changes and mechanical degradation. Therefore, controlling carbon stoichiometry is crucial for optimizing the microstructure and properties of carbides.

This principle also applies to HECCs where carbon stoichiometry strongly impacts mechanical properties. However, the inherent chemical complexity of HECCs introduces additional considerations. The diverse LAEs created by multiple metallic constituents allow for selective carbon deficiency to manifest in certain LAEs. Coupled with SRO formation that changes local atomic arrangement, the role of carbon non-stoichiometry in HECCs has another level of complexity in governing their mechanical properties. Our results reveal that while carbon deficiency may soften HECCs, the presence of SRO can enhance strength by promoting more favorable atomic arrangements. This interplay highlights the need for a nuanced understanding of how these factors collectively influence the mechanical behavior of HECCs. Our ML model, built upon local atomic features, effectively captures this interplay, enabling rapid evaluation of mechanical performance under both carbon non-stoichiometry and SRO simultaneously.

In summary, we have elucidated how changes in local atomic arrangements, induced either by carbon non-stoichiometry or SRO, can affect the mechanical properties of HECCs based on DFT, MC, and ML. We demonstrate that macroscopic mechanical properties, such as the elastic modulus, reflect atomic-level local interactions and can be viewed as a statistical average of the local properties. By predicting these local properties through ML, we can develop strategies for optimizing macroscopic mechanical behavior, such as adjusting the overall metallic composition or introducing SRO. In this way, our study bridges the gap between microscopic local environments and macroscopic mechanical performance, providing a theoretical foundation for the design and optimization of HECCs through strategic element selection and precise control of local chemistry and structure.

Methods

First-principles calculations

First-principles calculations were carried out using the Vienna ab initio simulation package (VASP) based on DFT^43,44. The Kohn-Sham equation was solved by the projector augmented wave (PAW) method⁴⁵, with standard PAW potentials selected for the constituent elements Ti, Zr, Hf, Nb, and C. The exchange-correlation interactions were modeled using the Perdew-Burke-Ernzerhof (PBE) form⁴⁶. The “mcsqs” code from the Alloy Theoretic Automated Toolkit (ATAT)⁴⁷ was employed to generate SQS structures for the considered MTMCs⁴⁸, each containing 64 atoms. For all structural optimization calculations, Brillouin-zone integrations were performed on a Γ-centered 5 × 5 × 5 k-point grid. The kinetic energy cutoff for plane waves was set to 600 eV, and the convergence criterion for electronic self-consistency was 10⁻⁵eV. Both lattice vectors and atoms were fully relaxed until the force components on atoms were smaller than 10⁻² eV/Å.

The vacancy formation energies and elastic constants of 420 carbon vacancy sites were obtained by DFT calculations. Specifically, the vacancy formation energies were calculated by:

$${E}_{V}^{f}={E}_{vac}-{E}_{per}+{E}_{carbon}$$

(7)

where E_vac is the energy of the supercell with one carbon vacancy, E_per is the energy of the perfect supercell, and E_carbon is the chemical potential of the carbon atom, assumed to be the energy per atom in the graphite phase here. The elastic constants with and without considering carbon vacancies were calculated using the strain-energy method. Specifically, a set of strains was applied to the B1 structure of MTMCs, and the corresponding changes in energy were recorded. The elastic tensors were determined by fitting a second-order polynomial to the strain-energy relations. The bulk modulus B, shear modulus G, elastic modulus E, and Poisson’s ratio ν were estimated based on the Voigt-Reuss-Hill relation⁴⁹.

Dataset and ML models

A total of 420 data were obtained through DFT calculations to construct the dataset for ML. Among them, there were six two-metal MTMCs: (TiZr)C, (TiHf)C, (TiNb)C, (ZrHf)C, (ZrNb)C, and (HfNb)C. For each system, two SQS supercells were generated, and within each supercell, fourteen different carbon vacancy sites were considered. The three-metal MTMCs encompassed four systems: (TiZrHf)C, (TiZrNb)C, (TiHfNb)C, and (ZrHfNb)C. Three SQS supercells were generated for each system, with fourteen carbon vacancy sites chosen in each supercell. In the case of the four-metal MTMC (TiZrHfNb)C, six SQS supercells were generated, and fourteen carbon vacancy sites were selected for each supercell. The increasing number of SQS supercells from 2-metal to 4-metal MTMCs was due to the growing complexity of local structural configurations along with the increase in the number of metallic components. Thus, for systems incorporating a larger variety of cations, an increased number of structures helps capture a broader range of local environments.

In this work, we have tested the performance of various ML algorithm for our purpose. ML models, including logistic regression (LR), least absolute shrinkage and selection operator (LASSO), ridge regression (RR), kernel ridge regression (KRR), support vector machine (SVM), random forest (RF), gradient boosting regression (GBR), extreme gradient boosting (XGB), and artificial neural network (ANN) were implemented with Scikit-learn⁵⁰. The hyperparameters in the models were optimized by the grid search method. The considered hyperparameters and details about selected ranges are listed in Table S1. The dataset was divided into an 80% training set and a 20% validation set to assess the accuracy of the models.

Metropolis Monte-Carlo (MMC) simulations

MMC simulations were conducted at 300 K using a 4 × 4 × 4 supercell containing 512 atoms. To explore the evolution of SRO in (TiZrHfNb)C, the mixing enthalpies of different cations were calculated by:

$$Delta {H}_{mix}[({A}^{i}{A}^{j})C]={E}_{DFT}(({A}^{i}{A}^{j})C)-{E}_{DFT}({A}^{i}C,G)-{E}_{DFT}({A}^{j}C,G)$$

(8)

where E_DFT((AⁱA^j)C) is the DFT total energy of (AⁱA^j)C with a rock-salt structure, E_DFT (AⁱC, G) and E_DFT (A^jC, G) are the DFT total energies of the AⁱC and A^jC carbides in their ground-state phases. The nearest neighbor model (NNM) was employed to derive the local mixing enthalpy of a cation within the supercell. Specifically, the mixing enthalpy of a cation depends on its interactions with surrounding atoms. The NNM assumes that only the 1^st nearest neighbor (NN) cations interact with the central cation and treats all the 1NN cations as independent entities. Hence, the local mixing enthalpy around a central cation was calculated by ref. ⁵¹:

$$Delta {H}_{local}=frac{1}{{N}_{local}}mathop{sum }limits_{j=1}^{{N}_{local}}Delta {H}_{mix}[({A}^{i}{A}^{j})C]$$

(9)

where N_local is the coordination number (N_local = 12). The total mixing enthalpy for the entire supercell was determined by summing the local mixing enthalpies of all N_total = 256 cations:

$$Delta {H}_{total}=mathop{sum }limits_{i=1}^{{N}_{total}}Delta {H}_{local}^{i}$$

(10)

The change in total mixing enthalpy between i and i−1 steps serves as an energy criterion for conducting MMC simulations in the (TiZrHfNb)C supercell:

$$Delta H=Delta {H}_{total}^{i}-Delta {H}_{total}^{i-1}$$

(11)

During the swap process of MMC simulations, if ΔH < 0, the trial configuration was directly accepted. However, if ΔH > 0, the trial configuration was accepted by a probability p_B determined according to the Boltzmann factor:

$${p}_{B}=Min{exp (frac{-Delta H}{{k}_{B}T}),1}$$

(12)

where k_B and T are the Boltzmann constant and temperature, respectively. The simulations were started from a quasi-random configuration and continued until the enthalpy reached equilibrium, which is observed around 3000 steps. The SRO parameter of the 1NN shell was calculated by ref. ^52,53:

$${alpha }_{ij}=1-frac{{p}_{ij}}{{c}_{j}}$$

(13)

where p_ij represents the conditional probability of species j appearing in the 1NN shell of species i, while c_j is the concentration of species j in the entire system.