A systematic study of oxide passivation induced by chemical composition and microstructure regulation in high entropy alloys

Abstract

The oxidation behavior of high entropy alloys (HEAs) has garnered growing attention owing to their potential applications in extremely harsh environments. Although some efforts have been made in newly developed data science methods, bottlenecks have arisen regarding the complex effect of multi-component alloying on the oxidation behavior. To extend the prior knowledge of the physical process of oxidation, we systematically analyze the elemental composition, local chemical environment, microstructure, and thermal oxidation behavior from a theoretical perspective, along with their relationships. We reveal that the chemical short-range order (CSRO) of Cu-containing HEAs repels oxygen occupation regardless of Cr's advantageous oxygen affinity and the octahedral interstitial sites where oxygen is prone to reside, which is validly attributed to the impact of the d orbital of the principal elements. Microstructures, such as grain boundary (GB), serve as rapid pathways for oxygen ingress into the matrix and tune the chemical composition in a refined manner. We thus summarize the series of factors that play essential roles in passive quality, including microstructure design, the consequent elemental regulation, and the diversity of elements in chemical affinity and diffusion rate. These insights will provide a fundamental understanding of oxidation resistance optimization strategies in HEAs.

---

Introduction

With the development of Accident Tolerant Fuels (ATF), such as U₃Si₂¹, the urgent need for a new material to replace Zircaloy on the outer surface of the cladding, with enhanced corrosion resistance at high temperatures, has become apparent. High entropy alloys (HEAs),^2,3 as an emerging class of advanced structural materials, have garnered widespread attention due to their advantages such as high strength,^4,5 high-temperature oxidation resistance,^6,7 corrosion resistance,^8,9 and wear resistance.¹⁰

For alloys with potential applications in high temperature scenarios, such as coating for reactor containment under extreme conditions, the oxidation resistance at operating temperatures is crucial. Holcomb et al.¹¹ studied the oxidation behavior of the Cantor alloy (CoCrFeMnNi) at intermediate temperatures and found that the long-term oxidation behavior at 650 °C and 750 °C in air was predominantly governed by the formation of Mn/Cr oxide layers. However, it has been observed that the presence of Mn can render alloy passivation unstable due to the fast-outward diffusion of Mn,^12,13 such as the formation of a porous oxide. Therefore, Mn was replaced with Cu, which has commonly been used in traditional alloys, such as steels and Mg–Al alloys, to enhance solid solubility and improve corrosion resistance.^14,15 Currently, the CoCrCuFeNi alloy¹⁶ has been successfully fabricated utilizing 3D printing lattice technology through a mixture of Co₃O₄, CuO, Fe₂O₃, NiO, and Cr₂O₃ oxide powders. Kai et al.¹⁷ compared the oxidation behavior of the CoCrCuFeNi alloy and the ternary alloy FeCoNi at temperatures ranging from 800 °C to 1000 °C, noting an improvement in oxidation resistance attributed to the presence of Cr and Cu. Mittireddi et al.¹⁸ found that at temperatures above 573 K, Sm₂(Co, Fe, Cu, and Zr)₁₇ mainly formed CuO, hindering the penetration of oxygen into the underlying layers and significantly reducing the oxidation process due to the formation of semi-passivating oxide. Besides, Cu-rich precipitates near the subsurface were observed after 3.0% Cu was added to an 18% Cr–8% Ni alloy.¹⁹

However, there is a consensus that the precipitation of Cu in HEAs can result in severe intergranular fracture.^20–22 The associated explanation of Cu precipitation in 3d transition HEAs has been made from a thermodynamic perspective, considering the positive mixing enthalpy between Cu and the other three elements (Co, Cr and Fe), and the relatively small negative mixing enthalpy with Ni.^23,24 However, this explanation concerning the positive mixing enthalpy resulting in element precipitation is deemed unreliable because the mixing enthalpy between Ni and the other 3d transition elements is negative,²⁵ yet experimental and theoretical observations indicate Ni segregation.^26–29 Thus, the thermodynamic interpretation based on enthalpy is not consistent with the experimental data. The concept of segregation enthalpy has guided extensive research on binary and ternary nanoscale alloys,³⁰ but it has not been applied to multi-component alloys. Some researchers have attempted to explain Cu intergranular clustering or sub-grain clustering from the perspective of melting point differences, but specific studies and analyses have not been conducted.³¹

Additionally, the abnormal segregation of Cr at grain boundaries (GBs) has been observed in CoCrCuFeNi, allowing the alloy to maintain nano-scale grains even after annealing at temperatures up to 1073 K.³² More importantly, Cr has been found to be one of the most significant elements in the capability of forming compact, dense and stable passive films.^16,33–35 Thermodynamic parameters, such as the enthalpy of mixing (ΔH_mix) and entropy (ΔS_mix), satisfactorily assess the phase stability of MPEAs on a global scale; however, it is challenging from a mathematical perspective to accurately represent the atomic-level details of elemental precipitation. The complex chemical environment represented by elemental ordering/clustering has been a focal point in the study of HEAs, necessitating a more reasonable mechanistic explanation.

Numerous studies have reported the thermodynamics and kinetics of the oxide process in HEAs. Experimentally, Huang et al.³⁶ investigated the oxidation kinetics, characteristics of oxide layers and formation of oxide layers in CoCrCu_1.2FeNi in solid and semi-solid states, attributing the formation of outer oxide layers to the pronounced outward diffusion of Cu ions and the formation of inner oxide layers to the internal diffusion of oxygen. Bai et al.³⁷ suggested that elemental diffusion during the oxidation of FeCoNiTiCu is mainly related to the Kirkendall effect, with the oxidation sequence influenced by atomic diffusion rates and the oxide formation energy of constituent elements. Theoretical calculations by Liu et al.³⁸ using first principles calculations found that octahedral positions with high local Cr concentrations or low local Ni concentrations in FCC CrCoFeNi high entropy alloys are favorable for the insertion of oxygen, indicating that chemical ordering plays a crucial role in oxidation in high entropy alloys. These results suggest that the complex chemical environment plays a significant role in the solution behavior of oxygen, and tuning the solution behavior can effectively impact the oxidation and corrosion processes of alloys.

This study primarily investigates the mechanisms behind the formation of local short-range order (SRO) induced by Cr and Cu, and how the local environment affects the occupation behavior of oxygen in different interstitial sites. Three derivative systems of Cantor alloys (CoCrFeNi, CoCuFeNi and CoCrCuFeNi) are selected owing to their distinct SRO. Subsequently, we elucidate the underlying mechanisms of SRO from the perspective of atomic stress, distinguishing them from the thermodynamic approach of the mixing enthalpy. For materials exposed to an oxygen atmosphere, Cr is a well-known element that forms an oxide barrier for protecting materials from corrosion. Meanwhile, Cu inhibits the catalytic activities of Fe, Ni, and Co. Based on the characteristics of these two elements and the local chemical environment inherent in the three alloys, we investigate the solution behavior of oxygen atoms. Finally, our research indicates that, despite significant differences between the oxygen affinities of Cr and Cu, the introduction of GBs with varied excess free volumes of lattice enhances the solution of oxygen as interstitials rather than the effect of oxygen affinity from element tuning. Overall, this study sheds light on the decisive roles of chemical affinity, microstructure and elemental diffusion rate in the passive quality during the oxidation process of HEAs, synergistically offering guidance for future alloy design.

Methods

First principles calculations

This study employs first principles calculations based on density functional theory (DFT) and the local self-consistent multiple scattering (LSMS) method, with atomic stress evaluations conducted on structures following thorough relaxation. DFT was performed to probe atomic physical properties on an electronic scale utilizing the Vienna ab initio simulation package (VASP).³⁹ The projector augmented wave (PAW) potential method was used to describe the interactions between ions and electrons, and the exchange correlation function was selected as the generalized gradient approximation (GGA) determined by Perdew–Burke–Ernzerhof (PBE).⁴⁰ For the bulk structures, a gamma-centered k-point grid of 4 × 4 × 4 was used, and for the Σ5 GB structure, the k-point grid was 11 × 4 × 1. To achieve a convergence finer than 0.01 meV per atom, we used a plane-wave basis energy set at 350 eV. To enhance system convergence, Methfessel–Paxton smearing with a Sigma value of 0.1 eV was employed. Proper sets of Γ-centered k-point meshes for each supercell, as outlined in Table S1, were utilized for Brillouin-zone integrations, which ensured total energy convergence within 1 × 10⁻⁵ eV and Hellmann–Feynman forces on all atoms below 0.02 eV Å⁻¹. Atomic positions and cell shapes underwent full relaxation in all calculations.

The crystal models were constructed as follows: for bulk structures, the 3 × 3 × 3 supercell used for HEA bulk stacking contained 108 atoms. Additionally, common Σ5(310)[001] GB (hereafter referred to as Σ5) was established using the Coincident Site Lattice (CSL) method.⁴¹ To initiate alloy structure simulations, the Special Quasi-random Structure (SQS) method⁴² was employed to achieve a fully random distribution of atoms in the FCC lattice, approximating the early ideal solid solution state of HEAs.⁴³ Because the SQS method is merely a mathematical approximation of the fully random state, the atomic size, interatomic preferential combination, and other factors are ignored in the alloy design of the SQS. Therefore, to obtain stable structures, our self-developed Metropolis Monte Carlo–Density functional Theory (MMC-DFT) algorithm⁴⁴ was employed to generate an equilibrium state in energy, experiencing 5000 MC steps. As shown in Fig. S1, all the samples reach the convergence of potential energy. Structural convergence can be reflected by the SRO parameter evolution with MC steps, as depicted in Fig. 1 and Fig. S2, S3.

Fig. 1 Elemental segregation/ordering in the bulk structures of CoCrFeNi and CoCuFeNi via MC-DFT calculations.

The temperature selected for MC calculations in this study was 800 K, which was considered in the range of practical application temperature,⁴⁵ and variations in material properties were observed at elevated temperatures.⁴⁶

MC-DFT calculations

We constructed a hybrid method that combines Monte Carlo (MC) and DFT. The hybrid MC-DFT code is self-compiled,⁴⁴ which works in the Gibbs semi-grand ensemble, allowing it to predict the equilibrium phases accurately, regardless of the number of components and the crystallographic structures. The MC method rooted in the Metropolis algorithm⁴⁷ is widely used to search for equilibrium states of liquids and solids along a Markov chain. Because DFT is included in this methodology, the advantage of MC-DFT is rooted in its comprehensive consideration of atomic physical properties on an electronic scale. The MC-DFT method bypasses the sluggish physical dynamics encountered in the system, such as diffusion phenomena while yielding average compositional distributions in thermodynamically equilibrated systems. Its core methodology accounts for atomic relaxations, vibrations, and mutual diffusion in the model. The method operates under the assumption that, before the exchange of positions between two atoms within the crystal, the total energy of crystal A, denoted as E₁, signifies the “old structure”. Subsequently, following the exchange of positions, the total energy of crystal B, denoted as E₂, represents the “new structure”. If E₂ is less than E₁, indicating a lower energy state and greater stability for the “new structure”, this optimized structure is accepted, with probability P_AB being 1. Conversely, if E₂ exceeds E₁, the new optimized structure is accepted with a probability, P_AB, defined as follows:where k_B is the Boltzmann constant, T is the thermodynamic temperature, and ΔE = E₂ − E₁. This iterative process persists, with the “new structure” supplanting the “old structure”.

Atomic stress calculations

LSMS calculations were based on the analysis of single-particle Green's function characteristics and the static properties of the local spin density approximation (LSDA) free energy function.⁴⁸ We calculated an approximate value of the local spin density of the electronic Green's function diagonal part. In this algorithm, multiple electron scattering processes were applied within each atom and in the Local Interaction Zone (LIZ) centered on clusters of M atoms. The fundamental assumption of this method is that an accurate approximation of the electron density, and state density at the atom site in the LIZ can only be obtained by performing multiple scattering calculations within the LIZ, thus providing an accurate estimate of the total system energy E. Atomic potentials were approximated using a spherical symmetry, and the von Barth and Hedin function was employed for the local approximation of the exchange–correlation energy. Subsequently, integrations over energy and electron density within Voronoi polyhedra were performed to obtain atomic-level information.

To accurately assess atomic stress, a series of volume static tensile strain experiments were conducted. For each structure subjected to volume stretching, the LSMS algorithm was employed to investigate atomic stresses. Subsequent relaxation optimization calculations of the LSMS algorithm were no longer performed on the relaxed structures from VASP. To achieve static stretching, perturbations are applied to the given system, and the energy of each atom is calculated at two volumes (with ±0.5% strain in three directions, and ±1.5% volume strain). The stress was then calculated by dividing the energy difference between the two strains by the volume change. Positive stress values indicate compressive (repulsive) pressure, while negative values indicate tensile (attractive) tension.

Warren–Cowley short-range order

The Warren–Cowley SRO parameters^49,50 were employed to characterize the degree of chemical ordering. The SRO parameters can effectively characterize the mixture of alloying elements and capture the transition from a solid solution to segregation. Describing the SRO in HEAs with n components is significant and more complex than in binary cases. Compared with binary alloys, where a single Warren–Cowley SRO parameter can describe the situation for each adjacent shell, an n-component alloy requires n(n − 1)/2 linearly independent SRO parameters for description. The expression for the order parameters is as follows:where p_ij represents the probability of finding an element j in the vicinity and c_j is the atomic concentration of element j. For a random solid solution, the SRO parameter is zero. Positive and negative values correspond, respectively, to a decrease and an increase in the tendency for pairs.

Solution energy

The formula for the solubility energy of oxygen atoms in the alloy is given bywhere E_matrix is the relaxed energy of the matrix without oxygen, E_matrix+O is the energy of the matrix containing oxygen, and E_O2 is the chemical potential of an oxygen molecule in the gas phase. Negative values indicate a preference for oxygen occupation at the specific interstitial site, while positive values signify the opposite.

First principles calculations indicated that both octahedron (OCT) and tetrahedron (TET) sites acted as stable sites for oxygen interstitials in the BCC MoNbTaW HEA.⁵¹ In this study, it was found that OCT and TET sites were also available for occupancy by oxygen interstitials in the FCC alloy, CoFeNi–X (X = Cr or Cu). In Σ5(310)[001], there are additional interstitial sites, such as pentagonal bipyramid (PBP), bitetrahedron (BTE) and capped trigonal prism (CTP).⁵² Detailed schematic diagrams are provided in Fig. S4.

Results and discussion

Local chemical environmental characteristics

The hybrid MC-DFT method was performed to search for the lattice equilibrium state of the systems studied here (Fig. S1 illustrates the energy vs. MC step for all the systems). The simulation started with an SQS structure. When the energy deviation between sequential acceptance states was smaller than 5 meV per atom, we observed that the Markov chain converged. Actually, when this convergence criterion was met, we still extended the iterative calculations to 5000 steps for further verification. We initially investigated bulk structures without GBs, in which the dispersion of the free volumes for interstitial polyhedrons was relatively smaller than that in the GB structure. The SRO parameters for Cr/Cu in CoFeNi–X (X = Cr or Cu) exhibited contrasting trends. Fig. 1 illustrates the trends of Cr–Cr and Cu–Cu pairing during MC iterations. It is evident from Fig. 1 that Cr–Cr atomic pairs are less prone to aggregation in the bulk structure, while Cr atoms prefer to occupy sub-lattice lattice sites. This observation aligns with the computational results by Mizuno et al.^53,54 In contrast, Cu atoms in CoCuFeNi exhibit a significant tendency towards clustering, which is consistent with numerous experimental findings.^55,56 Compared to the quaternary CoCuFeNi alloy, the quinary CoCrCuFeNi alloy exhibits a more pronounced trend of Cu segregation. Fig. S2 illustrates the evolution of SRO parameters in the CoCrCuFeNi bulk structure as a function of MC steps. Out of anticipation, Cr atoms no longer occupy the sub-lattice sites, which differs from the phenomenon occurring in CoCrFeNi. That is to say, the local chemical environment of the CoCrCuFeNi HEA is significantly affected by Cu atoms, reflected by the SRO of the Cr–Cr pair, which is consistent with a fully random distribution within the quinary HEA. The Cu–Cu SRO parameter approaches −2, still retaining the tendency of Cu clustering. This is demonstrated in the analysis of atomic stress.

Atomic stress associated with chemical short-range order

This atomic stress mechanism, integrating atomic size, structural characteristics, and electronic effects, can also influence the occupancy behavior of oxygen atoms in different interstitial sites. Our study involves a comparative analysis of CoFeNi–X (X = Cr or Cu) high entropy alloys. Both Cr and Cu belong to the 3d transition metal group and exhibit differences in electronegativity and atomic size. Utilizing Bader charge analysis,⁵⁷ we discovered that charge transfer of atoms in HEAs is strongly correlated with the inherent electronegativity difference of principal elements, as shown in Fig. 2(a) and (b). Herein, the electronegativity difference was calculated as the disparity between the Allen electronegativity⁵⁸ of each principal element and the average electronegativity in HEA. For instance, the electronegativities of the four principal elements within the CoCuFeNi HEA are closely aligned, resulting in minimal charge transfer between the atomic constituents. In the case of CoCrFeNi, the situation is disparate. Then, we obtained an average atomic radius of 1.39 Å for Cu in CoCuFeNi, while the average atomic radius for Cr was 1.34 Å in CoCrFeNi. This indicates that the atomic sizes in the alloy differ from the metallic volumes observed in ref. 59 (Cu: 1.27 Å; Cr: 1.30 Å).

Fig. 2 Charge transfer (histogram) and electronegativity difference (squared line) in the three studied HEAs: (a) CoCuFeNi, (b) CoCrFeN, and (c) CoCrCuFeNi.

Fig. 3 illustrates the atomic Voronoi volumes and atomic stress diagrams for CoCuFeNi and CoCrFeNi for one of the configurations at the equilibrium state, where the bar of the deviation from the average value marks the upper and lower limits in the fluctuation of the atomic volume. The variances in atomic volume and stress exactly reflect the feature of the equilibrium structure owing to the lattice distortion of HEA. Fig. S5 shows the variance of atomic volume and stress for another 4 equilibrium configurations for SRO and SQS, respectively. The 4 configurations are also obtained by randomly sampling the equilibrium state of MC steps. We can observe that Fig. S5 has the distribution of stress vs. volume and the order of elements is similar to those of the configuration depicted in Fig. 3, indicating that the inherent feature of atomic volume and stress dispersion is substantially reflected in our calculations. As shown in Fig. 3(a), the SQS and SRO bulk structures of CoCuFeNi exhibit a linear relationship between the average Voronoi volume and average atomic stress. In particular, the SRO structure acquired from the MC iteration calculation exhibits a linear fit, with a correlation coefficient very close to 1. This may be a consequence of close alignment for the four elements in the 3d periodic table with small charge transfer between elements. This phenomenon suggests that atoms with smaller volumes in the CoCuFeNi bulk are more prone to negative atomic stress (tension), while larger volumes are more susceptible to positive atomic stress (compression). Besides, the fluctuation in atomic Voronoi volume generally stems from the local lattice distortion caused by the inherent atomic size variations in HEAs.^60,61 The existence of local chemical ordering generally intensifies the lattice distortion in CoCuFeNi, as reflected in the broader distribution of Voronoi volumes for the four elements compared to the SQS structures. Simultaneously, the fluctuation in atomic stresses is shrunk for the SRO structures and the atomic stresses of different elements approach the average value of the whole system (referring to the dashed line of zero depicted in Fig. 3). Extremely high tensile and compressive stresses are alleviated in the SRO structure.

Fig. 3 Atomic Voronoi volumes and atomic stress diagrams in the bulk structures of CoCuFeNi and CoCrFeNi. (a) CoCuFeNi; (b) CoCrFeNi. Left corresponds to SQS structures, and the right corresponds to SRO structures.

We revealed the phenomenon of random/ordering of elements in multi-principal element alloys based on the relationship between atomic volume and atomic stress. This contributes to understanding element segregation in structures containing larger free space volumes, such as those occurring in grain boundaries. Fluctuations of atomic stress and Voronoi volume for individual elements in the Σ5 GB structures of CoCuFeNi and CoCrFeNi HEAs are different from those of the bulk structures owing to the presence of large excess free volume in GB, as depicted in Fig. 4 and Fig. S6, with another 4 samples at the equilibrium state. The features of these samples are similar. In particular, the centroid of the Cr stress distribution shifts upward toward zero in the SRO CoCrFeNi alloy containing GB compared with the one without GB (Fig. 3b). This is attributed to the alleviation of high atomic stress, including Cr negative and Ni positive stress. Therefore, the concentration of Cr is increased to 33% at the Σ5 GB. The Cr enrichment in the region of large free volume has been corroborated by experimental results under ion irradiation⁶² and computational simulations.^63,64 Xu et al. achieved a stable state of defect concentration and chemical ordering in a CoCrFeNi high entropy alloy through 2.4 MeV Cu ion irradiation. Significantly, Cr segregation occurred owing to the rearrangement of the chemical environment in CoCrFeNi caused by the abundance of vacancy defects generated by irradiation. Theoretically, simulations of ion irradiations on a CoCuFeNi alloy model, which was generated using a hybrid MC/MD method by Koch et al.,⁶³ supported these findings. The simulation results indicated that the post-irradiation structure approached a state of nearly random solid solution but traces of local precipitation remained. Specifically, Cu with the largest size and resulting high stress fills in vacancy-type defects, which provide a large excess volume created under irradiation. Besides, Fig. 5(a) and Fig. S7(a) (another 4 samples) depict the relationship between atomic stress and atomic Voronoi volume in the random state of the CoCrCuFeNi bulk structure. Compared with the SQS structure, Cu clustering in the SRO structure reduces the averaged atomic stress of Cu, whose average value approaches zero, as shown in Fig. 5(b) and Fig. S7(b). Moreover, Fig. 5(c) and Fig. S7(c) illustrate that the positive stress values for Cu atoms generally decline when introducing Σ5 GB and the enlargement of the Voronoi volume owing to lattice distortion. The summation of all atoms for each structure is shown in Fig. S7(d), indicating that the characteristic patterns of atomic stress and volume distribution are universal in the equilibrium configurations of the studied HEAs. Compared with the bulk structure, it is observed that the stress of Cr in Σ5 GB presents much lower values. Owing to the larger Voronoi volume of Cu occupation and Ni at GBs, the Cr atoms were pushed away from the GB region, reflected by the clustering phenomenon of the Cu–Ni pairs shown in Fig. S3.

Fig. 4 Atomic Voronoi volumes and atomic stress diagrams in Σ5 grain boundary structures of CoCuFeNi (left) and CoCrFeNi HEA (right).

Fig. 5 Relationship between atomic volume and atomic stress in CoCrCuFeNi. (a) SQS bulk structure, (b) SRO bulk structure, and (c) Σ5 SRO structure.

For the bulk structure of CoCrFeNi, most of the Cr atoms occupy sub-lattice positions in the SRO structure, thus forming an L1₂ structure. This phenomenon arises from a stress-induced mechanism in alloys. In the SRO structure, the average Bader atomic volumes for Cr, Fe, Co, and Ni are 10.01, 10.58, 10.83 and 11.15 ų, respectively, which increase with the electronegativity of the element. It can be observed from Fig. 3(b) that owing to the overall large tensile stresses of Cr among all constituent elements of CoCrFeNi, Cr preferentially occupies the second-nearest neighbors of other Cr atoms in the SRO structure rather than the first-nearest neighbors (in contrast to the Cu clustering phenomenon). Unlike the CoCuFeNi system, a clear linear relationship is not apparent between atomic stress and atomic Voronoi volume, as illustrated in Fig. 3(b). It is noteworthy that the components of CoCuFeNi are all 3d late transition metals, while Cr in CoCrFeNi is positioned in the middle of the 3d transition metal series, with outer electrons in a half-filled state. The difference in electronegativity between Cr and the other three elements is very large, which differs from the case of Cu, as shown in Fig. 2(a) and (b). The non-linear correlation in CoCrFeNi should be attributed to the relatively flexible nature of Cr with the smallest electronegativity compared with the other elements.

Oxygen solution behavior in quaternary HEAs

The relationship between the solution energies for the bulk structures of CoCuFeNi and CoCrFeNi and the volume variation of the interstitial polyhedron (ΔV/V) is illustrated in Fig. S8. It is apparent that the volume variation in interstitial polyhedrons caused by oxygen at the OCT position is relatively small. To better visualize the solution energy distribution in the matrix, we employed the box and normal statistical method, with which CoFeNi–X (X = Cr or Cu) bulk systems are conveniently compared, as shown in Fig. 6. Overall, the associated solution energies at the OCT positions are lower than those at the TET positions. The variance of the solution energies is small for OCT compared to when oxygen is at TET, implying that OCT is favorable for oxygen occupation. Because there are four nearest neighboring host atoms around a TET site, while an OCT site has only two, moving metal atoms out of the former site causes higher strain energy than out of the latter. It is evident that the variation in TET volumes (ΔV/V) is obviously larger than that of OCT for all systems, as shown in Fig. S8. Therefore, the oxygen atom is more favorable in terms of energy to squeeze into the OCT positions. It is noteworthy that a considerable number of TET positions in both alloys also have negative energies for oxygen solution, indicating the capture capability of TET positions for oxygen. Experimentally, Liu et al.⁶⁵ directly observed oxygen at the OCT positions and provided evidence of oxygen occupancy at the TET positions using multi-slice electron phase diffraction techniques in the TiNbZr alloy.

Fig. 6 Solution energy distribution of oxygen atoms at TET and OCT positions in the SQS and SRO structures of CoCuFeNi and CoCrFeNi alloy bulk systems, presented using a statistical style. (a) CoCuFeNi and (b) CoCrFeNi. The left panel corresponds to the SQS structures, and the right panel corresponds to the SRO structures.

Moreover, the upper quartiles of the solution energy at TET for both the SRO and SQS structures are consistently higher than those at OCT in CoCuFeNi, as shown in Fig. 6(a). However, in CoCrFeNi depicted in Fig. 6(b), the oxygen occupancy behavior at OCT positions is chemical ordering dependent, as shown in the box and normal plots of solution energy for the SRO structure. The average, upper/lower bounds, and upper/lower quartiles are significantly lower than those for the SQS structure. This representation directly demonstrates the trend of oxygen occupying the OCT sites in the SRO structure. Comparing CoCuFeNi shown in Fig. 6(a) and CoCrFeNi shown in Fig. 6(b), it is evident that the variation in the distribution of solution energies from the SQS structure to the SRO exhibits distinct behaviors in these two systems. The solution energies of oxygen at the TET and OCT of CoCuFeNi significantly overlap in both the SQS and SRO structures, which are rooted in the inherent lattice distortion of HEAs. In stark contrast, for CoCrFeNi, the overlap phenomenon of solution energies disappears and the solution energies of OCT and TET are discernible when the SQS structure is transformed into the SRO of the equilibrium state. As shown in Fig. 6, the average, upper/lower bounds, and upper/lower quartiles of OCT are significantly lower than those of TET. The statistics directly demonstrate the trend of oxygen occupying the OCT sites in the ordered structure of CoCrFeNi. This again reflects that Cr occupying sub-lattice positions in the SRO structure efficiently alleviates lattice distortion with respect to CoCuFeNi. In CoCrFeNi, the solution behavior of oxygen is, to some extent, similar to that of traditional alloys without any lattice distortion.

We quantitatively computed the distribution of host atoms within the first coordination shell around oxygen atoms at the TET and OCT sites to further analyze the chemical environment surrounding oxygen atoms at different interstitial sites. The relationship between the solution energy in the quaternary CoCuFeNi/CoCrFeNi alloy vs. coordinated element concentration around interstitial sites is presented in Fig. S9–S16. In the main text, we selected typical diagrams of the local atomic environment, as illustrated in Fig. 7. We can observe from Fig. S9–S12 that all Cu-containing systems present an increasing trend in the solution energy of oxygen as the number of coordinated Cu atoms increases, indicating a disfavor for an oxygen solution in the local environment with high Cu concentrations. Notably, as illustrated in Fig. 7(b), the solution energy of OCT in the SRO structures of CoCuFeNi HEA exhibits a pronounced increasing trend with the number of coordinated Cu, where a brief linear fit was applied as a guide to the eye. In contrast to Cu-rich regions, regions enriched with Cr are conducive to the solution of oxygen, particularly in interstitial sites with an ample free volume of OCT sites. Fig. S13–S16 demonstrates that as the local concentration of Cr atoms increases, the solution energy decreases. This suggests that regions with a high concentration of Cr atoms enhance the propensity to form oxides, thereby improving corrosion resistance. Experimental research has shown that Cr₂O₃ played the role of a protective barrier against further oxidation.^33,66 Intriguingly, Ni presents a different response in affinity for oxygen when the constituent Cu element is replaced with Cr. Either in SRO structures or SQS structures, an increase in the Ni concentration is found to be detrimental to the solution of oxygen in the Cr-containing systems, as illustrated in Fig. S13, while an inverse trend occurs in the Cu-containing systems.

Fig. 7 Typical example of element concentration in the first-nearest neighbour around oxygen atoms at different interstitial sites. OCT site in the SRO bulk structure of (a) CoCrFeNi and (b) CoCuFeNi.

Oxygen solution in quinary alloys

The tendencies of solution energy variation in the quinary CoCrCuFeNi vs. the surrounding Cr/Cu concentrations for all the CoCrCuFeNi systems are similar to the previous quaternary alloys. As depicted in Fig. S17–S20, it can be observed that solution energy increases with Cu concentration in the first-nearest neighbors of the interstitial site, while an increase in Cr concentration reduces solution energy. However, for the quinary systems, the solution energy of oxygen at either OCT or TET interstitial sites rarely exhibits negative values even in the environment of high Cr coordination. This implies that Cu participation even weakens the affinity of oxygen to Cr and further disadvantages the formation of passivating oxide concerning Cr, such as Cr₂O₃. The weakening of the Cr affinity is attributed to the d band variation in Cr in the presence of Cu, as shown in Fig. S21. It can be observed that the number of unfilled d-orbitals of Cr near the Fermi level obviously decreases after introducing Cu. Instead, Cr in CoCrFeNi has significant conduction band d-states around the Fermi level to be filled, enhancing the affinity of oxygen. Experimentally, Song et al.⁶⁷ studied a series of nanocrystalline Al_xCu_y(FeCrNiCo)_100−x−y HEAs, which suggests that reducing the Cu/FeCrNiCo ratio enhances the alloy's corrosion resistance, supporting our results to some extent. Therefore, the precise adjustment of the constituent elements significantly affects the oxidation resistance and corrosion resistance of the alloys.

Oxygen solution in GB structures

The investigations of the oxygen solution above are primarily focused on the bulk structure of HEAs. We know that oxidation and corrosion always occur at GBs or interfaces. A recent experiment⁶⁸ demonstrated that Cu segregation at the alloy interface and near GBs can reduce the oxygen partial pressure to a level where the Cr element cannot be directly oxidized to form Cr₂O₃. GB structures are regarded as large lattice distortions that possess a large excess free volume. With respect to the bulk structure, it is necessary to conduct a systematic study on the effect of GBs on the solution of oxidation and the oxygen affinity for different elements.

The solution energy of oxygen in the bulk structures and the Σ5 GB structures of CoCrFeNi, CoCuFeNi, and CoCrCuFeNi HEAs is summarized in Fig. 8, where the solution energies at OCT and TET in the bulk structure are denoted in red, and those at diverse sites in the Σ5 GB structure are denoted in blue. Again, we can observe from Fig. 8 that the occupancy probability of oxygen at TET is lower than that at OCT in the bulk structure for all the HEAs. The average solution energies of the OCT sites in CoCuFeNi and CoCrCuFeNi are positive, contrasting with the solution energy in CoCrFeNi, which is due to the difference in the affinity of oxygen for Cr and Cu. For the Σ5 GB structures, even in Cu-containing HEAs, all interstitial sites, such as BTE, CTP and PBP, favor the accommodation of the oxygen atom with negative solution energy. The solution energies of oxygen at different interstitial positions at GB are still influenced by the local arrangement of the different elements, as evidenced by the differences in the solution energy at the BTE interstitial sites in these HEAs.

Fig. 8 Oxygen solution energy at different interstitial sites in different structures of CoCrFeNi, CoCuFeNi and CoCrCuFeNi HEAs. Note: red represents OCT and TET sites in the bulk structure, while blue represents BTE, CTP, and PBP sites in the grain boundary structures.

Effect of inherent electronic structure on oxygen solution

The comparison of the solution energies depicted in Fig. S8 indicates that the Cu-HEA is generally unfavorable in energy for oxygen occupancy both at TET and OCT sites, while Cr-HEA can accommodate oxygen, specifically in OCT sites of the SRO structure. This may be because Cu occupies a large atomic volume and a full d shell with fairly inactive valence electrons and oxygen addition causes larger strain energy. However, the small volume of Cr with flexible d electrons facilitates the insertion of oxygen into adjacent OCT sites. Note that the volume of OCT is larger than that of TET in general. This local strain mechanism affects oxygen atoms and is also applicable to the solubility behaviors of the other elements in HEAs. Our earlier analysis of the solution energy for hydrogen atoms at TET and OCT interstitial sites in HEAs indicates that the mechanical contributions to the solution energy resulting from lattice expansion are comparable with the chemical contributions arising from compositional differences.⁶⁹

To unveil the inherent mechanism from a chemical perspective, we calculated the d band centre, which is defined as the centre of the projected density of states for the d electrons of a metal element relative to the Fermi level. It can be employed to assess the impact of the atomic interactions within the d band, such as the adsorption and solution of atoms. The expression for the d band centre is given by

where x and ρ(x) represent the energy and density of the d orbital states, respectively. Here, ε_d is characterized as a scaling relationship for solution energy. It can be clearly observed from Fig. 9 that ε_d of Ni in the Cu-containing HEA shifts right and the density of d states near the Fermi level is obviously high compared with the Cr-containing HEA. This behavior remarkably enhances the affinity of oxygen to Ni in CoCuFeNi. The in-depth mechanism should be traced back to the variation in the relative electronegativity of Ni with respect to Cu and Cr, which causes different charge transfers, as shown in Fig. 2.

Fig. 9 d DOS of Ni in CrCoFeNi and CuCoFeNi HEA.

Effect of GBs on oxidation

The aforementioned analysis reflects that serving as rapid pathways for oxygen ingress into the matrix, the oxygen solution at GBs is relatively unaffected by the specific elemental composition. Utilizing a grain boundary fracture model, Chan⁷⁰ predicted embrittlement and oxidation behavior in Ni-based alloys, which indicates that oxygen in air at 600 °C rapidly diffuses along GBs, reducing the cohesive strength of the alloy and leading to intergranular cracking. However, the formation of oxides at the GBs, to some extent, regulates subsequent oxygen diffusion into the matrix, mitigating the corrosive effect of oxygen on the alloy. Hu et al.⁷¹ observed larger grain-sized Al₂O₃ in Hf-doped AlCoCrFeNiYTa HEAs, effectively impeding oxygen diffusion inward along GBs. In Cr-containing HEAs, apart from the formation of a passive Cr₂O₃ film at GBs,³³ carbides, such as Cr₂₃C₆ phases, are also found to hinder the outward diffusion of metal ions and the inward diffusion of oxygen as carbon atoms infiltrate the coating during the oxidation process.⁷²

However, unlike the stable passive films of Cr₂O₃, elements with larger atomic sizes, such as Ni and Cu, appearing near GBs tended to form less stable oxides. For example, the oxides of Ni appearing on the surface or the outer layer have been confirmed to be less stable compared with Cr₂O₃.^37,73 Similarly, the oxides, including CuO and Cu₂O, in the outer oxide layer have been found to fall out as the oxidation temperature and reaction time increased.³⁶

Preferred oxidation

Now, we can depict a scenario of the oxidation process in HEAs by combining the analysis of the GB effect and the solution energy of oxygen without considering other factors. Given the element segregation due to alleviating the high stress, Cu and Ni tend to segregate at the alloy surface, interface and GBs and form oxides in the initial stages of oxidation (as confirmed by ref. 37 and 66); when oxygen penetrates only the matrix, owing to the superior oxidative activity of Cr compared to Fe, Co, Ni and Cu, oxygen reacts with Cr atoms at the forefront of the oxidation zone, producing a passive Cr₂O₃ film.

Building upon the elucidation of the influence of crystal structure on oxidation behavior discussed earlier, we need to explore the preferred oxidation behavior of multicomponent alloys based on the intrinsic oxygen affinity of the elements to better understand the thermodynamic processes of oxidation in HEAs. Here, we employ a change in Gibbs free energy to comprehend the preferred oxidation. According to the research by Backman et al.,^74–76 the change in the free energy of oxidation reactions determines the oxidation products, and slight differences in the free energy change between different oxidation reactions can lead to significant variations in selective oxidation. Based on the experimentally measured oxidation products,^37,77Fig. 10 illustrates the changes in Gibbs free energy related to oxidation reactions of the five elements Co, Cr, Cu, Fe and Ni calculated using the FactSage program.⁷⁸ It is noteworthy that the Gibbs free energy change here represents the energy change resulting from the oxidation reaction of 1 mole of oxygen with the metal. For example, at 800 °C, the oxidation reaction of 1 mole of oxygen with 4/3 mole of Cr yields a Gibbs free energy change of −565.995 kJ, expressed as follows:

Fig. 10 Change in Gibbs free energy for different oxidation reactions.

As depicted in Fig. 10, the sequence of Gibbs free energy changes for oxidation products, from smallest to largest, is Cr₂O₃, FeO, Fe₃O₄, Fe₂O₃, CoO, NiO, Co₃O₄, Cu₂O, and CuO. This order aligns with the affinity sequence of the elements for oxygen dissolution in this study. From this ranking, it is inferred that a lesser number of d orbital electrons in 3d transition metals is advantageous for their affinity towards oxygen. In alloy systems, differences in the affinity of metals for oxygen induce competitive oxidation reactions. For instance, as discussed earlier, the variation in the coordinated element species around Ni affects the shifting of its d-band center, consequently influencing the dissolution of oxygen around Ni atoms within different alloys. In the bulk structure, we investigate the change in Gibbs free energy resulting from oxidation reactions, which originate from the distinct oxygen affinity characteristics of different elements, a conclusion supported by recent studies. Bai et al.³⁷ demonstrated through surface segregation energy calculations for FeCoNiTiCu constituents that Cu and Ni tend to reside within the matrix during oxidation, contrary to their surface presence under oxygen-free conditions.

The oxidation process also exhibits both kinetic attributes. The outward diffusion of metals is a kinetic process that follows Wagner theory, covering the parabolic law⁷⁹ within the oxidation mechanism, closely related to the diffusion coefficients of the alloying element.⁶⁶ Thus, the different diffusion rates of the principal elements make a difference in the order of the oxides. For instance, the rapid diffusion of Cr in the FeCrNi alloy facilitated the formation of Cr₂O₃ layers through the BCC and FCC phases, thereby promoting the oxidation resistance of HEAs.³³ However, it was reported that the diffusion coefficient of Cu was faster in the CoCrCuFeNi HEA, resulting in an outer layer mainly composed of CuO, which inhibits the probability of the formation of Cr₂O₃. Besides, the kinetic process is somewhat influenced by the temperature. Taking the Al-containing HEA as an example, the temperature was too low for the Al atoms to diffuse fast enough to form protective oxides.⁸⁰

Inspiringly, our results, based on the elemental chemical character involving the electronic interaction and atomic-level interaction, which are essentially associated with thermodynamic Gibbs free energy, show that the preferred oxidation, which determines the passive quality of HEAs, is influenced by the chemical affinity of the element. Moreover, apart from the inherent character of the previous two factors (kinetics and thermodynamics), not only does the microstructure, such as GB, provide channels for oxygen atoms, but structure-induced segregation also facilitates oxygen to preferentially react with the adjacent atoms segregated at GBs, which is the primary focus of our current research efforts. Therefore, diverse oxidation products can be observed. In summary, the driving forces from chemical affinity, microstructure and diffusion rate play essential roles in the passive quality of HEAs, as demonstrated in Fig. 11. In-depth investigations are needed because the synergistic and competitive effects of these processes at different oxidation stages remain vague in the field of HEAs. We employ multiscale calculations to investigate this crucial issue further, thereby offering theoretical insights into experimental endeavors, such as optimizing the fabrication of high entropy nanoparticles and high entropy oxides by incorporating the Kirkendall effect.

Fig. 11 Schematics of the essential factors determining the passive quality of HEAs.

Discussion on experimental observation

Identifying the SRO structure in experiments has been fairly challenging for multi-principal element alloys. In practice, chemical short-range order is the most difficult to decipher in experiments. Mössbauer spectroscopy can distinguish the SRO phenomenon. Unfortunately, it is limited to binary alloys, particularly in the Fe–Cr system. The reason is that Mössbauer spectroscopy can probe SRO in Fe-containing binary alloys on the premise that the neighboring atom (the species different from Fe) around the ⁵⁷Fe nuclei probe is assumed to be measurable by analyzing the variation in spectral hyperfine parameters. However, for multi-principal elements, the spectral hyperfine parameters become too complex to distinguish individual element effects around the ⁵⁷Fe nuclei probe. Recently, Chen et al.⁸¹ stated that SRO was observed using transmission electron microscopy (TEM), where V or Cr forms a sublattice in a ternary alloy (VCoNi or CrCoNi). However, in a rapid sequence, Walsh et al.⁸² refuted it and claimed that superlattice intensities occurring in pure elements cannot represent chemical ordering in TEM observations. This is because the measurement process may be affected by local fluctuations, leading to uncertainties regarding the plausibility of the conclusions drawn. As the number of constituent elements increases, the difficulty in characterizing and identifying the SRO structure severely increases using this kind of crystallographic method. To date, direct evidence for observing the SRO in CoCrFeNi and CoCuFeNi, or other HEAs composed of more constituent elements, is still missing. We also expect more novel experimental methods to feasibly and explicitly demonstrate the SRO phenomenon in multi-principal element alloys in the future.

Here, we focused on the oxygen solution and its behavior in the different microstructures (the grain boundary and intragranular region) to compare the microstructural effects. Generally, the formation of oxide passivation in multi-principal element alloys is associated with selective oxidation. In principle, the oxidation process is driven by the combined effects of thermodynamic (e.g. oxidation Gibbs free energy) and kinetic factors (e.g. elemental diffusion rates), that is, the two factors determine the preferential oxidation of elements in HEAs. In limited experiments on CoCrFeNi and CoCuFeNi systems, CoCrFeNi has been reported to possess high corrosion resistance, which was attributed to the presence of Cr.⁸³ Our calculation showed that Cr has the lowest Gibbs free energy among oxidation reactions with different elements, as shown in Fig. 10, which supports experimental observations from the thermodynamic viewpoint. Regarding chemical composition variation at grain boundaries, it is experimentally challenging to elucidate whether it originates primarily from elemental segregation (thermodynamic factor) or differences in elemental diffusion rates (a kinetic factor) since accurately identifying the SRO structure is fairly challenging, as mentioned above. Obviously, the high diffusivity of element atoms preferentially promotes the formation of the corresponding oxide for oxidation at the grain boundaries. In HEAs, when elemental atoms preferentially segregate at grain boundaries (this also results in chemical composition variation), they analogously promote the preferential oxidation reaction when oxygen is introduced along the grain boundaries. Unfortunately, few experimental studies have not provided enough evidence to decipher these factors, resulting in the preferential oxidation of elements at grain boundaries. In any case, the experimental work observed a grain-boundary phase that was enriched in Cr in a CoCrFeNi HEA,^84–86 which supported our calculations. Moreover, an experimental study on tracer diffusion in CoCrFeNi indicated that Cr diffusion is the fastest among all elements.⁸⁷ Cr enrichment at grain boundaries, along with its affinity to oxygen and fast diffusion, facilitates the formation of a passivation film.

Conclusions

This study systematically investigates the relationship between the elemental composition, chemical short-range order, microstructure and oxidation behavior of high entropy alloys, comparing and analyzing the differences in oxygen affinity among different elements and their effects on other elements. The introduction of a large free volume in crystal structures promotes oxygen dissolution in complex multicomponent alloys. Considering element segregation at GBs and the interface, and the affinity of oxygen, we differentiate selective oxidation behaviors in the HEAs studied herein. Specifically, the findings are as follows:

(1) The Cr element in CoCrFeNi and the Cu element in CoCuFeNi exhibit contrasting chemical short-range orders. Atomic stress regulates the segregation behavior of atoms of different sizes.

(2) Overall, the octahedral interstitial sites are more suitable for oxygen atom occupancy in the bulk structures compared to the tetrahedral interstitial sites due to the former's larger spatial volume.

(3) Regions with a high concentration of coordinating Cr atoms facilitate oxygen atom occupancy, while Cu-rich regions repel oxygen atom occupation. The dissolution of oxygen related to the affinity of coordinated Ni atoms is subject to modulation, either inhibition due to the effect from surrounding Cr or promotion from surrounding Cu.

(4) In CoCrCuFeNi, severe enrichment of Cu elements affects the tendency for Cr atom occupancy in sub-lattice lattice sites. Additionally, Cu insensitivity to oxygen challenges oxygen dissolution to form a high-quality passive film, such as Cr₂O₃, in this alloy system.

(5) High free volume, such as in grain boundary regions, provides ample space for large-sized principal atoms while facilitating oxygen atom dissolution. This indicates that grain boundaries serve as rapid pathways for oxygen ingress into the matrix and diffusion and are not significantly influenced by the types of elements present.

(6) We summarize that a synergetic effect of chemical affinity, microstructure and diffusion rate exists in oxidation associated with the passive quality of HEAs.

Author contributions

Jianqiao Yu and Panhua Shi contributed equally to this work: conceptualization, methodology, validation, writing – original draft, data curation, investigation, writing – review and editing. YiYing Yang: data curation and investigation. Rongchuan Li: investigation. Jiaxuan Si: validation and investigation. Fan Yang: investigation. Yuexia Wang: conceptualization, writing – review and editing, and supervision.

Conflicts of interest

There are no conflicts to declare.

Data availability

The data supporting this article have been included in the SI.

Potential energies with the evolution of MC steps, SRO parameters, interstitial sites, atomic Voronoi volumes and atomic stress diagrams, solution energies, d DOS, and DFT caculations parameters. See DOI: https://doi.org/10.1039/d5cp02698a.

Acknowledgements

This work was supported by the National Natural Science Foundation of China under Grant No. 12275055 and U2067218.

Notes and references

1. S. P. Shu, Y. B. Miao, B. Ye, K. Mo, L. Jamison, S. Bhattacharya, A. Oaks, A. M. Yacout, J. Harp, L. A. Nimmagadda and S. Sinha, J. Nucl. Mater., 2023, 587, 154734.
2. J. W. Yeh, S. K. Chen, S. J. Lin, J. Y. Gan, T. S. Chin, T. T. Shun, C. H. Tsau and S. Y. Chang, Adv. Eng. Mater., 2004, 6, 299–303.
3. B. Cantor, I. T. H. Chang, P. Knight and A. J. B. Vincent, J. Mater. Sci. Eng. A, 2004, 375, 213–218.
4. C. H. Lee, F. Maresca, R. Feng, Y. Chou, T. Ungar, M. Widom, K. An, J. D. Poplawsky, Y. C. Chou, P. K. Liaw and W. A. Curtin, Nat. Commun., 2021, 12, 5474.
5. Z. M. Li, K. G. Pradeep, Y. Deng, D. Raabe and C. C. Tasan, Nature, 2016, 534, 227.
6. M. Y. Jiao, Z. F. Lei, Y. Wu, J. L. Du, X. Y. Zhou, W. Y. Li, X. Y. Yuan, X. C. Liu, X. Y. Zhu, S. D. Wang, H. H. Zhu, P. P. Cao, X. J. Liu, X. B. Zhang, H. Wang, S. H. Jiang and Z. P. Lu, Nat. Commun., 2023, 14, 806.
7. T. M. Butler, K. J. Chaput, J. R. Dietrich and O. N. Senkov, J. Alloys Compd., 2017, 729, 1004–1019.
8. N. Birbilis, S. Choudhary, J. R. Scully and M. L. Taheri, npj Mater. Degrad., 2021, 5, 14.
9. M. Li, H. Henein, C. Zhou and J. Liu, J. Mater. Sci. Technol., 2024, 174, 133–144.
10. J. S. Luo, W. T. Sun, R. X. Duan, W. Q. Yang, K. C. Chan, F. Z. Ren and X. S. Yang, J. Mater. Sci. Technol., 2022, 110, 43–56.
11. G. R. Holcomb, J. Tylczak and C. Carney, JOM, 2015, 67, 2326–2339.
12. X. Q. Wu, S. Xu, J. B. Huang, E. H. Han, W. Ke, K. Yang and Z. H. Jiang, Mater. Corros., 2008, 59, 676–684.
13. Y. X. Xu, W. Y. Li, L. Z. Qu, X. W. Yang, B. Song, R. Lupoi and S. Yin, J. Mater. Sci. Technol., 2021, 68, 172–183.
14. C. D. Marioara, S. J. Andersen, J. Royset, O. Reiso, S. Gulbrandsen-Dahl, T. E. Nicolaisen, I. E. Opheim, J. F. Helgaker and R. Holmestad, Metall. Mater. Trans. A, 2014, 45a, 2938–2949.
15. X. Ye, M. Ma, Y. Cao, W. Liu, X. Ye and Y. Gu, Phys. Proc., 2011, 12, 303–312.
16. D. C. Zhang, C. Kenel and D. C. Dunand, Acta Mater., 2022, 238, 118187.
17. W. Kai, W. L. Jang, R. T. Huang, C. C. Lee, H. H. Hsieh and C. F. Du, Oxid. Met., 2005, 63, 169–192.
18. R. T. Mittireddi, R. K. Singh and E. Panda, Corros. Sci., 2023, 224, 111518.
19. J. Yan, P. Zhang, Y. Gu, Y. Yuan, B. Da, L. Ma and Y. Li, Corros. Sci., 2020, 169, 108601.
20. H. T. Zheng, R. R. Chen, G. Qin, X. Z. Li, Y. Q. Su, H. S. Ding, J. J. Guo and H. Z. Fu, J. Mater. Sci. Technol., 2020, 38, 19–27.
21. Z. H. Zhong, Q. Xu, K. Mori and M. Tokitani, Philos. Mag., 2019, 99, 1515–1527.
22. T. He, Y. Qi, Y. Ji and M. Feng, Int. J. Mech. Sci., 2023, 238, 107828.
23. X. Wang, H. Xie, Y. Deng, D. Zhu, X. Yang, D. Dong and T. Ma, J. Alloys Compd., 2023, 938, 168511.
24. K. N. Campo, C. C. de Freitas, L. Fanton and R. Caram, J. Mater. Sci. Technol., 2020, 52, 207–217.
25. M. C. Troparevsky, J. R. Morris, P. R. C. Kent, A. R. Lupini and G. M. Stocks, Phys. Rev. X, 2015, 5, 011041.
26. Y. Wang, J. Yin, H. Y. Hou, S. Y. Hu, X. B. Liu, F. Xue, C. H. Henager and J. T. Wang, Nucl. Instrum. Methods Phys. Res., Sect. B, 2019, 456, 32–36.
27. M. N. Hasan, J. Gu, S. Jiang, H. J. Wang, M. Cabral, S. Ni, X. H. An, M. Song, L. M. Shen and X. Z. Liao, Scr. Mater., 2021, 190, 80–85.
28. D. Litvinov, K. J. Lu, M. Walter and J. Aktaa, Mater. Today Commun., 2023, 35, 106030.
29. C. M. Barr, J. E. Nathaniel, K. A. Unocic, J. Liu, Y. Zhang, Y. Wang and M. L. Taheri, Scr. Mater., 2018, 156, 80–84.
30. T. Chookajorn, H. A. Murdoch and C. A. Schuh, Science, 2012, 337, 951–954.
31. M. Xu, L. Han, C. Shen, Q. Jiang, G. Yang, C. Tang, A. Gharbi, W. Xu and J. Yi, Vacuum, 2024, 219, 112723.
32. J. Cao, F. Li, Q. Zhang, M. Li, C. Wang, W. Wang and Y. Liu, Scr. Mater., 2023, 234, 115545.
33. H. Sun, M. Li, H. Zhang and J. Liu, Corros. Sci., 2024, 227, 111685.
34. X. Duan, T. Han, X. Guan, Y. Wang, H. Su, K. Ming, J. Wang and S. Zheng, J. Mater. Sci. Technol., 2023, 136, 97–108.
35. T. A. Listyawan, M. P. Agustianingrum, Y. S. Na, K. R. Lim and N. Park, J. Mater. Sci. Technol., 2022, 129, 115–126.
36. M. J. Huang, J. F. Jiang, Y. Wang, Y. Z. Liu and Y. Zhang, Corros. Sci., 2021, 193, 109897.
37. H. Bai, R. Su, R. Z. Zhao, C. L. Hu, L. Z. Ji, Y. J. Liao, Y. N. Zhang, Y. X. Li and X. F. Zhang, J. Mater. Sci. Technol., 2024, 177, 133–141.
38. Y. Liu, G. P. Zheng and M. Li, J. Alloys Compd., 2020, 843, 156060.
39. J. Hafner, J. Comput. Chem., 2008, 29, 2044–2078.
40. J. P. Perdew, K. Burke and M. Ernzerhof, Phys. Rev. Lett., 1996, 77, 3865–3868.
41. H. Grimmer, Acta Crystallogr., 1976, 32, 783–785.
42. A. Zunger, S. H. Wei, L. G. Ferreira and J. E. Bernard, Phys. Rev. Lett., 1990, 65, 353–356.
43. J. W. Yeh, S. K. Chen, J. Y. Gan, S. J. Lin, T. S. Chin, T. T. Shun, C. H. Tsau and S. Y. Chang, Metall. Mater. Trans. A, 2004, 35a, 2533–2536.
44. X. L. Ren, P. H. Shi, W. W. Zhang, X. Y. Wu, Q. Xu and Y. X. Wang, Acta Mater., 2019, 180, 189–198.
45. M. Rieth, S. L. Dudarev, S. M. G. de Vicente, J. Aktaa, T. Ahlgren, S. Autusch, D. E. J. Armstrong, M. Balden, N. Baluc, M. F. Barthe, W. W. Basuki, M. Battabyal, C. S. Becquart, D. Blagoeva, H. Boldyryeva, J. Brinkmann, M. Celino, L. Ciupinski, J. B. Correia, A. De Backer, C. Domain, E. Gaganidze, C. Garcia-Rosales, J. Gibson, M. R. Gilbert, S. Giusepponi, B. Gludovatzj, H. Greuner, K. Heinola, T. Höschen, A. Hoffmann, N. Holstein, F. Koch, W. Krauss, H. Li, S. Lindig, J. Linke, C. Linsmeier, P. López-Ruiz, H. Maier, J. Matejicek, T. P. Mishra, M. Muhammed, A. Muñoz, M. Muzyk, K. Nordlund, D. Nguyen-Manh, J. Opschoor, N. Ordás, T. Palacios, G. Pintsuk, R. Pippan, J. Reiser, J. Riesch, S. G. Roberts, L. Romaner, M. Rosinski, M. Sanchez, W. Schulmeyer, H. Traxler, A. Ureña, J. G. van der Laan, L. Veleva, S. Wahlberg, M. Walter, T. Weber, T. Weitkamp, S. Wurster, M. A. Yar, J. H. You and A. Zivelonghi, J. Nucl. Mater., 2013, 432, 482–500.
46. G. Bonny, A. Bakaev, R. Olsson, C. Domain, E. E. Zhurkin and M. Posselt, J. Nucl. Mater., 2017, 484, 42–50.
47. N. Metropolis, A. W. Rosenbluth, M. N. Rosenbluth, A. H. Teller and E. Teller, J. Chem. Phys., 1953, 21, 1087–1092.
48. D. M. C. Nicholson, G. M. Stocks, Y. Wang, W. A. Shelton, Z. Szotek and W. M. Temmerman, Phys. Rev. B: Condens. Matter Mater. Phys., 1994, 50, 14686–14689.
49. J. M. Cowley, Phys. Rev., 1950, 77, 669–675.
50. P. A. Flinn, Phys. Rev., 1956, 104, 350–356.
51. A. J. Samin, J. Appl. Phys., 2020, 128, 215101.
52. M. F. Ashby, F. Spaepen and S. Williams, Acta Metall., 1978, 26, 1647–1663.
53. M. Mizuno, K. Sugita and H. Araki, Results Phys., 2022, 34, 105285.
54. M. Mizuno and H. Araki, Phys. Rev. Mater., 2024, 8, 013601.
55. Suprianto and C. L. Chen, Met. Mater. Int., 2023, 29, 420–428.
56. Y. Yu, N. N. Xu, S. Y. Zhu, Z. H. Qiao, J. B. Zhang, J. Yang and W. M. Liu, J. Mater. Sci. Technol., 2021, 69, 48–59.
57. R. F. W. Bader and R. F. Bader, Atoms in Molecules: A Quantum Theory, Clarendon Press, 1990.
58. J. B. Mann, T. L. Meek, E. T. Knight, J. F. Capitani and L. C. Allen, J. Am. Chem. Soc., 2000, 122, 5132–5137.
59. T. Egami and Y. Waseda, J. Non-Cryst. Solids, 1984, 64, 113–134.
60. Q. F. He, J. G. Wang, H. A. Chen, Z. Y. Ding, Z. Q. Zhou, L. H. Xiong, J. H. Luan, J. M. Pelletier, J. C. Qiao, Q. Wang, L. L. Fan, Y. Ren, Q. S. Zeng, C. T. Liu, C. W. Pao, D. J. Srolovitz and Y. Yang, Nature, 2022, 602, 251–257.
61. G. Yao, W. Y. Wang, C. Zou, K. Ren, P. Li, X. Gao, D. Lin, J. Wang, S. Yang, Y. Wang, H. Song and J. Li, J. Am. Ceram. Soc., 2022, 105, 4260–4276.
62. Q. Xu, H. Guan, S. Huang, Z. Zhong, H. Watanabe and M. Tokitani, J. Alloys Compd., 2022, 925, 166697.
63. L. Koch, F. Granberg, T. Brink, D. Utt, K. Albe, F. Djurabekova and K. Nordlund, J. Appl. Phys., 2017, 122, 105106.
64. S. Chae, L. Williams, J. Lee, J. T. Heron and E. Kioupakis, Npj Comput. Mater., 2022, 8, 95.
65. C. Liu, J. Z. Cui, Z. Y. Cheng, B. Z. Zhang, S. Y. Zhang, J. Ding, R. Yu and E. Ma, Adv. Mater., 2023, 35, e2209941.
66. M. F. Li, H. Zhang, Y. M. Zeng and J. Liu, Acta Mater., 2022, 240, 118313.
67. Y. P. Song, L. C. Yan, X. L. Pang, Y. J. Su, L. J. Qiao and K. W. Gao, Corros. Sci., 2023, 213, 110983.
68. G. Huang, T. Chou, S. Liu, J. Ju, J. Gan, T. Yang, X. Feng, Y. Shen, J. Gao, M. Li and F. Meng, Appl. Surf. Sci., 2024, 652, 159285.
69. X. L. Ren, P. H. Shi, B. D. Yao, L. Wu, X. Y. Wu and Y. X. Wang, Phys. Chem. Chem. Phys., 2021, 23, 27185–27194.
70. K. S. Chan, Metall. Mater. Trans. A, 2015, 46a, 2491–2505.
71. J. Hu, C. Gu, J. Li, C. Li, J. Feng and Y. Jiang, Corros. Sci., 2023, 212, 110930.
72. C. Pan, X. Li, H. Luo, S. Qu, C. Yang and D. Zhu, Corros. Sci., 2022, 194, 109956.
73. Q. Zhao, X. Huang, Z. X. Zhan, S. L. Zhou, J. X. Liu, P. C. Zhu, L. H. Wei, X. Y. Li, C. Y. Li and Y. J. Xie, Corros. Sci., 2023, 220, 111291.
74. L. Backman, J. Gild, J. Luo and E. J. Opila, Acta Mater., 2020, 197, 20–27.
75. L. Backman, J. Gild, J. Luo and E. J. Opila, Acta Mater., 2020, 197, 81–90.
76. L. Backman and E. J. Opila, J. Eur. Ceram. Soc., 2019, 39, 1796–1802.
77. B. Song, Y. Yang, M. Rabbani, T. T. Yang, K. He, X. B. Hu, Y. F. Yuan, P. Ghildiyal, V. P. Dravid, M. R. Zachariah, W. A. Saidi, Y. Z. Liu and R. Shahbazian-Yassar, ACS Nano, 2020, 14, 15131–15143.
78. C. Bale, P. Chartrand, S. A. Degterov, G. Eriksson, K. Hack, R. Ben Mahfoud, J. Melançon, A. D. Pelton and S. Petersen, Calphad, 2002, 26, 189–228.
79. F. Y. Cho, S. W. Tung and F. Y. Ouyang, Mater. Chem. Phys., 2022, 278, 125678.
80. Z. Peng, J. Sun, H. W. Luan, N. Chen and K. F. Yao, Intermetallics, 2023, 155, 107845.
81. X. F. Chen, Q. Wang, Z. Y. Cheng, M. L. Zhu, H. Zhou, P. Jiang, L. L. Zhou, Q. Q. Xue, F. P. Yuan, J. Zhu, X. L. Wu and E. Ma, Nature, 2021, 592, 712.
82. F. Walsh, M. Zhang, R. O. Ritchie, A. M. Minor and M. Asta, Nat. Mater., 2023, 22, 926–929.
83. X. Shang, Z. Wang, F. He, J. Wang, J. Li and J. Yu, Sci. China: Technol. Sci., 2017, 61, 189–196.
84. S. Guo, C. Ng, Z. J. Wang and C. T. Liu, J. Alloys Compd., 2014, 583, 410–413.
85. K. C. Hoi, W. Lei, Y. N. Liu, C. H. Shek, J. T. G. Ferreira, N. F. T. Cortez, C. T. Kwok, Y. Sun, V. A. M. Cristino and K. Lo, Wear, 2023, 530, 204990.
86. R. Zhou, W. Y. Chen, W. P. Li, T. H. Chou, Y. H. Chen, X. P. Liang, J. H. Luan, Y. T. Zhu, J. C. Huang and Y. Liu, npj Mater. Degrad., 2023, 7, 8.
87. M. Vaidya, K. G. Pradeep, B. S. Murty, G. Wilde and S. V. Divinski, Acta Mater., 2018, 146, 211–224.

---