npj | computational materials Article

Published in partnership with the Shanghai Institute of Ceramics of the Chinese Academy of Sciences

https://doi.org/10.1038/s41524-024-01322-6

Tunable interstitial and vacancy diffusivity

by chemical ordering control in CrCoNi
medium-entropy alloy
Check for updates

Yangen Li , Jun-Ping Du , Shuhei Shinzato & Shigenobu Ogata

In this study, we utilized a quantitative atomistic analysis approach to investigate the impact of
chemical ordering structures on the diffusion behavior of interstitials and vacancies within the CrCoNi
medium entropy alloy (MEA), employing an advanced neural network interatomic potential (NNP). We
1234567890():,;
1234567890():,;

discovered that the degree of chemical ordering, which can be precisely controlled through annealing
at elevated temperatures, significantly influences both interstitial and vacancy diffusion. This
phenomenon contributes to the notable sluggish diffusion characteristic of CrCoNi, largely
attributable to the restriction of diffusion pathways in regions with lower degree of chemical ordering.
We also emphasized the crucial role of operating temperature on diffusion, which should be remained
well below the annealing temperature to preserve the sluggish diffusion effect. Our research sheds
light on the interplay between chemical ordering and defect diffusion in MEAs, and it proposes
effective strategies for tailoring the diffusivity of MEAs by altering their chemical ordering. These
insights are instrumental in the development of next-generation materials, which are optimized for use
in challenging environments, such as high-temperature and irradiation conditions.

Multi-principal element alloys, which encompass both medium-entropy necessarily sluggish and can even be enhanced at a given temperature19.
alloys (MEAs) and high-entropy alloys (HEAs), have garnered significant Additionally, recent theoretical work20 has suggested that there is no direct
attention due to their outstanding mechanical properties1–3. Recently, their relationship between high configurational entropy and sluggish diffusion in
exceptional radiation resistance4–7 has also come into focus, as they exhibit MEAs and HEAs. These diverse findings underscore the need for further
the ability to withstand harsh environments, such as high temperatures, research to fully comprehend the phenomenon of sluggish diffusion in
high pressures, and intense irradiation, for extended periods. Consequently, MEAs and HEAs.
they are being considered as potential next-generation materials for nuclear Currently, some research has highlighted the formation of chemical
structures8–10. During irradiation, the bombardment by high-energy parti- ordering structures in MEAs and HEAs after annealing at elevated
cles creates numerous interstitials and vacancies11. The evolution of these temperatures21–23, which could potentially influence the behavior of inter-
defects through atomic diffusion leads to the formation of vacancy clusters, stitial and vacancy diffusion as well as the radiation resistance24–27. Notably,
dislocations, and voids, ultimately resulting in materials failure12,13. The an experimental study indicated that the local chemical order in CrCoNi can
exceptional radiation resistance of MEAs and HEAs is believed to stem from delay the formation and evolution of interstitial and vacancy, leading to
their unique defect diffusion kinetics, characterized by sluggish diffusion improved radiation resistance at relatively low irradiation doses25. Atomistic
dynamics14–17, which can impede the rapid formation of large defects, simulations supported this by showing the local chemical order can reduce
thereby preventing materials failure. interstitial and vacancy diffusivity through increased diffusion energy bar-
However, the precise mechanism behind this sluggish diffusion in riers. Nevertheless, it is vital to recognize that these simulations are quali-
MEAs and HEAs remains unclear and is a subject of debate. Consequently, tative, primarily due to the difficulties in capturing the genuine chemical
the ability to manipulate and control sluggish diffusion has not yet been ordering structures and diffusion energy barriers in chemically ordered
achieved. For instance, Tsai et al. reported that experimental results CrCoNi using the embedded atom method (EAM) potential. Some other
demonstrate a slower diffusion coefficient and higher activation energy in vacancy diffusion simulation works have also pointed out that the formation
CoCrFeMnNi HEAs compared to reference metals18. In contrast, Vaidya of chemical short-range order can reduce vacancy diffusion by decreasing
et al. concluded that diffusion in CoCrFeNi and CoCrFeMnNi HEAs is not jump frequency and enhancing diffusion correlation28,29. However, the

Department of Mechanical Science and Bioengineering, Osaka University, Osaka 560-8531, Japan.
e-mail: ogata@me.es.osaka-u.ac.jp

npj Computational Materials | (2024)10:134 1

https://doi.org/10.1038/s41524-024-01322-6 Article

specific correlation between genuine chemical ordering structures and to an overestimation of the initial vacancy concentration, especially when
defect diffusion pathways has not been carefully explored. The precise compared to the experimental Time–Temperature–Electric resistance
nature of these chemical ordering structures is still debated, and a unan- relation32. In our current study, we took a different approach. We determine
imous conclusion remains elusive. Not to mention the study how the the initial vacancy concentration that best aligns with the experimental
modification of defect diffusion can be achieved by controlling the forma- Time–Temperature–Electric resistance relation, which relies on electric
tion of these chemical ordering structures. The chemical ordering effect on resistance measurements of the CrCoNi MEA at various annealing times
defect diffusion deserves an in-depth understanding. and temperatures (see Supplementary Fig. 1).
Recently, we successfully developed a precise neural-network potential In Fig. 1, the degree of chemical ordering is quantified using the
ij
(NNP) for CrCoMi MEA grounded in the density functional theory (DFT) Warren–Cowley order parameters31 of the 1st neighbor shell, α1 , with both i
produced training dataset. The NNP enabled us not only to specify realistic and j representing Cr (refer to the “Method” section for details). Each
local chemical ordering structures but also to quantitatively predict their C-shaped curve within the diagram stands for a specific chemical ordering
formation kinetics in CrCoNi MEA30,31. Moreover, the NNP accurately degree. The plateau of each C-shape, in relation to time, indicates the
replicates the diffusion kinetics of interstitials and vacancies, largely because equilibrium degree of chemical ordering reached at a given temperature.
it well-reproduce the DFT activation energy for the interstitial and vacancy Significant changes in chemical ordering can be observed within a reason-
diffusion. able annealing time at T = 673 K as shown in Fig. 1. At higher temperatures
In the present study, utilizing the NNP-driven atomistic simulations, than 673 K, chemical ordering is expected to be limited and not strongly
such as molecular dynamics (MD) simulations for interstitial diffusion and affects vacancy and interstitial diffusion, despite rapid chemical ordering
artificial-neural-network-accelerated diffusion-kinetic-Monte Carlo formation. In contrast, at temperatures lower than 673 K, the development
(ANN-kMC) simulations31 for vacancy jumps, we first precisely depict the of chemical ordering is considerably slower, though strong chemical
genuine chemical ordering structures relative to annealing temperature and ordering will eventually form. Therefore, our subsequent analysis of inter-
duration in CrCoNi MEA. Following this, we analyze its impact on inter- stitial and vacancy diffusion focuses on models annealed at T = 673 K for
stitial and vacancy diffusion, revealing a clear link between the degree of various durations. Hereafter we referred to the RSS model and the annealed
chemical ordering and the diffusion kinetics. Our findings indicate that the models as Model_RSS and Model_AN 1 ~ Model_AN 5, respectively.
genuine chemical ordering structures play a crucial role in defect sluggish Next we discerned the chemical ordering structures via a chemical
I
diffusion by restricting interstitial and vacancy diffusion regions. In addi- ordering matching analysis31, yielding a chemical ordering degree, M i d for
tion, the changes in diffusivity depending on the operational temperature of specific ordering type Id (Id = types I, II, and III) around each atomic site i
alloys was discussed, which is pivotal for understanding how alloys behave (see the “Method” section for the matching analysis details). The three
under practical operating conditions. Consequently, our results offer ave- chemical ordering types (Id = types I, II, and III, as depicted in Fig. 2a) were
nues for adjusting the diffusion kinetics of interstitials and vacancies by detailed in ref. 31. In these structures, Cr-rich and Co/Ni-rich layers are
tuning chemical ordering through annealing conditions. oriented along {110}, {100}, and {113} for types I, II, and III respectively.
Specifically, type I displays an Immm symmetry, with a propensity for Cr
Results atoms separation. Type II is reminiscent of the L10 structure, where Cr
Control of chemical ordering structures by thermal annealing atoms are inclined to cluster in one layer and Co/Ni atoms in another. Type
In our previous work31, we introduced a Time–Temperature–Chemical- III is akin to the L11 structure, where Cr atoms locate at the eight vertices (see
ordering-degree (TTC) diagram for the equimolar CrCoNi MEA using Supplementary Fig. 2 for the distribution of each chemical ordering struc-
ANN-kMC simulations combined with NNP. This diagram enables control ture in all models). The remaining structures all together such as other than
over the degree of chemical ordering by adjusting the annealing time from a the types I, II, and III are labeled “other” structures here. It should be noted
random solid solution (RSS) state at a fixed temperature. For the sake of that “other” structures also contain some weak chemical ordering, so they
clarity in the subsequent discussions, we have revisited this diagram in Fig. 1 should be distinguished from the RSS structures (see Supplementary Fig. 3
utilizing the computational strategy previously published in ref. 31 with a for details). Figure 2b displays distribution of the chemical ordering sites i
I
modification of initial vacancy concentration by using experimental electric having maxId ¼I;II;III M i d >0:85 in an equimolar RSS and models annealed for
resistance change data32. It is important to note that the previous TTC 4.4, 13.2, 209, 293, and 407 h from the RSS state at T = 673 K. This figure
diagram in reference31 was constructed under the assumption that the illustrates the spatial and temporal evolution of chemical ordering. The
CrCoNi MEA had cooled down from 1500 K, and the equilibrium vacancy chemical ordering increases with annealing time, with a noteworthy
concentration at 1500 K represented the initial vacancy concentration just observation being that after 407 h of annealing, the volume fraction of
before annealing commenced. However, this assumption seems to have led chemical ordering approaches nearly 50% (see also Fig. 3b).
Figure 3a, b quantifies the evolution over time of the Warren–Cowley
order parameters and the volume fraction of regions exhibiting type Id
chemical ordering, represented as f Ivd ¼ mId =N. Here mId denotes the
I
number of sites with M i d > 0.85, and N represents the total atom count in the
model. As depicted in Fig. 3a, the 407-h anneal resulted in the Cr-Cr and Co-
Ni pairs exhibiting positive Warren–Cowley order parameters of αCrCr 1 = 0.38
and αCoNi
1 = 0.25, respectively, indicating repulsion tendencies between Cr-Cr
and Co-Ni in the first nearest neighbor. Conversely, the Cr-Co and Cr-Ni
pairs manifested negative parameters of αCrCo 1 = −0.22 and αCrNi
1 = −0.14,
respectively, indicating attraction tendencies between Cr-Co and Cr-Ni in the
first nearest neighbor. These observations align with the DFT results33. Figure
3b further reveals that the volume fractions of types I and II chemical ordering
structures expand with annealing duration. Notably, after 407 h of annealing,
the type I structure’s volume fraction approached 39.6%, and substantial
chemical domain structures appeared in CrCoNi. In contrast, type III
diminished, likely due to its inherent instability and gradual decomposition
Fig. 1 | Time–Temperature–Chemical-ordering-degree diagram. Star shows state over extended annealing, as suggested by Zhou et al.34 and attributed to the
of models RSS and Model_AN 1 ~ Model_AN 5 annealed at T = 673 K for various energy penalty from increased Cr-Cr interactions31.
durations.

npj Computational Materials | (2024)10:134 2

https://doi.org/10.1038/s41524-024-01322-6 Article

Fig. 2 | Chemical ordering structures and their

distribution. a Structure information of three
representative chemical ordering (the red spheres
represent Cr atoms, and the gray spheres represent
Co/Ni atoms). b Distribution of chemical ordering
structure after annealing at T = 673 K for different
time (please refer to the text for details).

Fig. 3 | Evolution of Warren-Cowley order para-

meters and volume fraction of chemical ordering.
a Warren–Cowley order parameter changes during
annealing simulation, dash lines represent Mod-
el_AN 1 ~ Model_AN 5 with different annealing
time. b Volume fraction of types I, II, and III che-
mical ordering structures at different annealing
time. The label “other” represents volume fraction of
structure does not contain the types I, II, and III.

Impact of chemical ordering on interstitial and vacancy diffusion see Supplementary Fig. 4a, for details on chemical ordering variations
Interstitials are usually formed in radiation damage due to the high for- throughout the diffusion simulation).
mation energy. Materials will form a large number of point defects after One the other hand, vacancy diffusion is activated at high temperature,
radiation damage, including interstitials and vacancies. Since the high and atoms will continue to exchange positions with vacancy through vacancy
migration barrier of vacancy, interstitial diffusion plays an important role in diffusion mechanism. These behaviors usually occur during the annealing
the evolution of irradiation defects. Interstitials prefer to form dumbbells in process. Due to the different interaction between elements, a certain chemical
CrCoNi and migrate along with the dumbbell in [100] direction and form ordering structures will be formed after long-term annealing26. Therefore, for
new dumbbell to complete interstitial migration. Most of interstitials will the vacancy diffusion simulation, due to the higher migration energy barrier,
combine with vacancies and annihilate, and the remaining interstitials will we adopted the same simulation methodology as that used for the annealing
aggregate to form clusters after long-term diffusion, resulting in elements simulation, namely ANN-kMC combined with NNP (refer to “Method”
segregation within material17. For the interstitial diffusion simulation in this section for details). Here, a vacancy was generated by removing an atom from
work, an atom was arbitrarily incorporated into models with a count of the model, and the simulation temperature was set to T = 1200 K. A total of
N = 2916 atoms. The diffusion temperature was established at T = 1200 K to 2000 accepted vacancy jumps were executed in ANN-kMC simulation, and
ensure an adequate diffusion process during our 1 ns MD simulation time. similar to the interstitial diffusion simulation, no substantial change in the
Notably, no remarkable change in the existing chemical ordering were existing chemical ordering were detected within this MD simulation time
observed during this simulation time window, avoiding the potential window, thus avoiding the impact of the change in chemical ordering
influence of the segregation/ordering change on the diffusion results (please structures formed by vacancy diffusion (please see Supplementary Fig. 4b).

npj Computational Materials | (2024)10:134 3

https://doi.org/10.1038/s41524-024-01322-6 Article

Fig. 4 | Interstitial and vacancy diffusivities.

a Interstitial and b vacancy MSD - time curves in
models with different degrees of chemical ordering.
c Interstitial and d vacancy diffusivity in models with
different degree of chemical ordering.

It is worth mentioning that both prolonged interstitial and vacancy Supplementary Table 1). It is discernible that diffusivity exhibits a mono-
diffusion lead to changes in element distribution, such as segregation and tonic decline with amplifying chemical ordering degree, leading to sluggish
chemical ordering, but with different kinetic mechanisms. Since interstitial diffusion.
diffusion is much faster than vacancy, timescale required to change element
distribution varies. Here, both interstitial diffusion MD simulations and Correlation between chemical ordering region and diffuse region
vacancy jump kMC simulations were executed on the Model_RSS and To more precisely probe the relationship between the diffusive region of
Model_AN 1 ~ Model_AN 5, and chemical ordering impact on interstitial interstitial and vacancy and areas exhibiting chemical ordering, we under-
and vacancy diffusion was studied. took 50 independent diffusion simulations (both interstitial MD and
To quantify the diffusion behavior of both interstitials and vacancies, vacancy ANN-kMC diffusion) for the Model_RSS and Model_AN 5
we calculated the average mean squared displacement (MSD) of all atoms models. Each simulation had varied initial positions for the interstitials and
for these models. The gradients of the MSD-time plots were employed to vacancies.
calculate diffusivity (for more details on MSD calculation, please see the For a finer understanding, we sectioned the model volume into a
“Method” section). It should be noted that the calculated diffusivity depends detailed mesh (see the “Method” section for the details). We then calculated
on the defect concentration in the simulation, and the concentration of the spatial distribution for both the chemical ordering degree and the
interstitial and vacancy in this work are set to c ¼ N1 , where N = 2916 is the residence densities of interstitials and vacancies during the simulation. Note
number of atoms in the defect free model. For self-interstitial defects, as that since interstitials evolve into a dumbbell pair structure as they diffuse,
mentioned before, they are usually formed under high energy irradiation we utilized the coordinate centers of these dumbbell pairs to pinpoint the
environment due to the high formation energy of self-interstitial. Therefore, interstitial’s position.
the relationship between the actual self-interstitial concentration and the As shown in Fig. 5a, a (010) plane was selected to display the chemical
calculated diffusivity will not be further discussed. For vacancy defects, the ordering degree and interstitial and vacancy residence density results for
actual concentration depends on conditions such as the vacancy formation better understanding in the following discussions. Figure 5b, c renders the
energy and the local temperature environment. If the actual vacancy con- distribution
P ofItotal chemical ordering degree of type I + type II + type III,
centration (cact) can be obtained, e.g., by fitting the experimental data as Pk ¼ III P
Id ¼I k
d
on the (010) plane and “other”, Pok ¼ (1  Pk Þ on the (010)
I
described in Supplementary Fig. 1, one can correct the actual vacancy dif- plane, respectively, derived from Model_AN 5. Pkd is an average chemical
act
fusivity as Dact ¼ cccal Dcal according to the calculated vacancy diffusivity ordering degree of type Id in a mesh k (please see also “Method” section).
(D ) and concentration (ccal).
cal
Concurrently, Fig. 5d, e renders interstitial and vacancy residence densities
Figure 4a, b displays the MSD results for interstitials and vacancies, of mesh k; Qk on the (010) plane (please see also “Method” section). The
respectively. The calculated diffusivities for both interstitials and vacancies outcomes underscore a heterogeneous distribution pattern for interstitial
across all models are shown in Fig. 4c, d. It is evident that all models after and vacancy residence densities, with pronounced densities discerned in
annealing (Model_AN 1~ Model_AN 5) exhibit lower MSD gradients specific zones. This suggests a propensity for interstitials and vacancies to
compared to the Model_RSS. Additionally, these MSD slopes diminish with selectively diffuse within these zones. Contrarily, the Model_RSS exhibits a
escalating chemical ordering or prolonged annealing. For instance, the more homogeneous residence density, signaling a non-selective, random
Model_RSS yielded interstitial and vacancy diffusivities of 3.6 × 10−12 m2 s−1 diffusion pattern across regions (see Supplementary Fig. 5, for details on
and 3.6 × 10−13 m2 s−1, respectively. In contrast, the thoroughly annealed diffusion results in the Model_RSS).
Model_AN 5 recorded diffusivities of just 2.0 × 10−12 m2 s−1 and Upon juxtaposing the residence density distribution of interstitials and
1.2 × 10−13 m2 s−1 for interstitials and vacancies, respectively (please see also vacancies with the chemical ordering degree distribution (type I + type

npj Computational Materials | (2024)10:134 4

https://doi.org/10.1038/s41524-024-01322-6 Article

Fig. 5 | Distribution of chemical ordering degree and defects residence density. chemical ordering structures (type I + type II + type III) and “other” structures in
a A (010) plane was selected to display the chemical ordering degree and interstitial the Model_AN 5 on the (010) plane, respectively. d, e Interstitial and vacancy
and vacancy residence density results. b, c Chemical ordering degree distribution of residence density in the Model_AN 5 on the (010) plane, respectively.

II + type III chemical ordering structures and “other” structures), a pro- exhibit elevated formation energies for interstitials and vacancies. This
minent concurrence emerges between mesh regions boasting higher resi- results in a tendency for these defects to migrate toward regions with less
dence densities for interstitials and vacancies and those showing chemical ordering. Furthermore, the MSDs in both chemical ordering and
pronounced chemical ordering degrees. This infers that defects gravitate “other” structures were evaluated, as illustrated in Fig. 7. This analysis
toward “other” regions devoid of higher degree of chemical ordering reveals a slower diffusivity in the chemical ordering structures than in the
structures. Put differently, the inception of chemical ordering structures can “other” structures.
delimit the diffusion pathways for defects. For additional confirmation, refer These results point to a significant conclusion: the main reason for the
to Supplementary Figs. 6 and 7, which depict outcomes from two different decreased diffusion of both interstitials and vacancies within chemical
models bearing the same chemical ordering degree Model_AN 5 but ordering structures can be attributed to constraints in regions characterized
annealed from different RSS configurations, and Supplementary Figs. 8–10 by low defect formation energy, such as the “other” structure regions with
for those from models with chemical ordering degree that are less pro- limited chemical ordering. Additionally, even when interstitials and vacan-
nounced than in Model_AN 5. These results also support above conclusion. cies transition from the “other” region to the chemical ordering region, their
To quantify the interplay between diffusion regions for interstitials and diffusivity significantly decreases, preventing quick diffusion. As a result, the
vacancies and the spread of chemical ordering structures, we computed overall diffusivity slows down in relation to the degree of chemical ordering.
correlation factors between them. This was done using the relation between The findings robustly indicate that by manipulating the chemical
average chemical ordering degree Pk , and residence density Qk , expressed as ordering degree, achieved by adjusting the annealing duration at a given
P P Id
G ¼ a n1 nk¼1 III Id ¼I P k Qk , where n (=36 × 36 × 36) is the total number of
annealing temperature, we can deftly regulate the diffusivity of both inter-
mesh regions and a is a normalization factor (see “Method” section for the stitials and vacancies. It is worth noting that as demonstrated in Fig. 1, to
details). achieve a specific degree of chemical ordering and the corresponding target
diffusivity, one can select the annealing temperature to match the nose-
Figure 6a–d presents distribution of the correlation for both interstitials temperature of the iso-degree curve for the target degree in the TTC dia-
and vacancies, differentiated by chemical ordering structures (type I + type gram. This approach is aimed at minimizing the annealing time.
II + type III) and “other” structures on the (010) plane. The total correlation Note that when we compare the MSD results of the “other” structures
factors for these groups were found to be: G = 0.41 for interstitials with (Fig. 7) with those of the RSS structures (Fig. 4), we can find that the
chemical ordering structures, G = 0.38 for vacancies with chemical order- diffusivity of both vacancies and interstitials is certainly slower in the “other”
ings structures, G = 0.59 for interstitials with “other” structures, G = 0.62 for structures. This suggests that even the “other” structures contain some
vacancies with “other” structures. These results underline a stronger cor- degree of weak chemical ordering and its diffusivity is slower.
relation of residence densities of both interstitials and vacancies with “other”
structures. It reaffirms that the preferred diffusion regions for these defects Impact of operating temperature on diffusion and chemical
are those with less chemical ordering degree. ordering
To unravel why these diffusion-favorable regions emerge, we com- In our simulations, even though diffusion processes were modeled at high
puted the formation energy distributions for interstitials and vacancies in temperatures, they operated in short bursts. This was to ensure we analyzed
the Model_AN 5 on the same (010) plane. As showcased in Fig. 6e, f (with diffusivity within a fixed state of chemical ordering. But when we use MEAs
more detailed histograms in Supplementary Fig. 11 for both Model_RSS and and HEAs over prolonged durations at elevated temperatures, diffusion can
Model_AN 5), regions with pronounced chemical ordering structures alter the existing chemical ordering structures.

npj Computational Materials | (2024)10:134 5

https://doi.org/10.1038/s41524-024-01322-6 Article

Fig. 6 | Correlation between chemical ordering degree and defects residence correlation factors between chemical ordering degree of chemical ordering struc-
density. a, b Correlation factor distribution of interstitial and vacancy for chemical tures and “other” structures and interstitial (vacancy) residence density were cal-
ordering structures (type I + type II + type III) in the Model_AN 5 on the (010) culated as 0.41 (0.38) and 0.59 (0.62), respectively. e, f Distribution of interstitial and
plane, respectively. c, d Correlation factor distribution of interstitial and vacancy for vacancy formation energy in the Model_AN 5 on the (010) plane, respectively.
“other” structures in the Model_AN 5 on the (010) plane, respectively. The total

The interplay of the annealing temperature and the operating tem- notable disparity in chemical ordering degree between the initial and the
perature can be clearly depicted in the Time–Temperature–Chemical- equilibrium states.
ordering-degree diagram, illustrated in Fig. 8. The starting state of the In essence, if the operating temperature exceeds the annealing tem-
MEA, post-annealing at a particular temperature and time, is denoted by a perature, the sluggish diffusion effect progressively diminishes. On the other
black dot. The blue arrows depict scenarios where T o > T a . Here, when the hand, when the operating temperature is below the annealing temperature,
MEA is suddenly subjected to a higher operating temperature, it rapidly the sluggish diffusion effect becomes incrementally more pronounced.
achieves the thermal equilibrium for that temperature, driven by the
diffusion of interstitials and vacancies. Conversely, the red arrows Discussion
represent situations where T o < T a . In such cases, upon initiation of With the progress of MEAs and HEAs research, the impact of chemical
operation, the MEA gravitates toward increasing its chemical ordering to ordering on their performances has gradually attracted attention. It is widely
match the equilibrium structure of the operating temperature. This believed that the formation of chemical ordering can suppress defect dif-
transition is also facilitated by interstitial and vacancy diffusion but tends fusion, and its impact mainly lies in the diffusion barrier and correlation
to be prolonged due to reduced diffusivity at lower temperatures and the coefficient. For example, a MD simulation result concludes that chemical

npj Computational Materials | (2024)10:134 6

https://doi.org/10.1038/s41524-024-01322-6 Article

Fig. 7 | Mean-Squared-Displacement (MSD) -

time curve of defect diffusion. a Interstitial diffu-
sion and b vacancy diffusion in chemical ordering
(type I + type II + type III) and “other” structures in
the Model_AN 5, respectively.

Here, we explored the influence of chemical ordering on interstitial and

vacancy diffusion in CrCoNi MEA, illustrating the potential to modulate
this diffusivity by adjusting the chemical ordering through controlling
annealing conditions, such as temperature and duration. Our simulations
show that distinct chemical ordering structures emerge in CrCoNi MEA
when annealed at certain temperatures. These structures effectively repel
interstitials and vacancies, resulting in a restricted diffusion region and,
consequently, slower effective diffusion rates of entire system. Notably,
diffusivity correlates directly with the degree of chemical ordering, which in
turn is influenced by the annealing duration at a given temperature. This
underscores the potential to regulate the diffusivity of interstitials and
vacancies by modulating the degree of chemical ordering, achieved by fine-
tuning the annealing temperature and duration. We conclude by high-
lighting the pivotal role of operating temperatures on diffusion, emphasizing
Fig. 8 | Schematic diagram of the relationship between chemical ordering degree that to preserve the slow diffusion effect, the operating temperature should
αCrCr
1 and annealing temperature T a and operating temperature T o . One remain well below the annealing temperature.
C-shaped curve stands for the same chemical ordering degree. The longer the annealing
time or the lower the annealing temperature, the higher the chemical ordering degree.
Methods
Kinetic Monte Carlo simulation
The annealing process in CrCoNi was modeled using the vacancy jump
kinetic Monte Carlo (kMC) simulation. One vacancy was created in the
model with FCC structure containing 2916 atoms. Vacancy will jump to one
of its 12 nearest neighbor lattice sites, which is called
 as one
 kMC step. The
jump frequency was calculated as vi ¼ v0 exp  kΔETi , where the pre-
B

exponential factor v0 = 4.16 × 1013 s−1 was adopted by fitting the MD

vacancy diffusion simulation (details of the fitting are discussed below), kB is
Boltzmann’s constant, T is the absolute temperature, and ΔEi is the energy
barrier for a vacancy jump to the nearest neighbor site i, which is given by
artificial neural network (ANN) predictions. The ANN predictions can
significantly accelerate vacancy jump kMC simulations, especially when
performing long kMC simulations31.
In order to determine the pre-exponential factor v0 for the kMC
vacancy diffusion simulation, long-term vacancy diffusion MD simulation
was executed. The MD vacancy diffusion temperature was set to T = 1200 K,
and total simulation time is 10 ns. For the MD simulations, vacancy diffu-
sivity in Model_RSS and Model_AN 1 ~ 5 are calculated as
3.2 × 10−13 m2 s−1, 2.2 × 10−13 m2 s−1, 2.1 × 10−13 m2 s−1, 1.7 × 10−13 m2 s−1,
1.6 × 10−13 m2 s−1 and 1.1 × 10−13 m2 s−1, respectively. The modified pre-
Fig. 9 | Mean-Squared-Displacement (MSD)-time curve of vacancy diffusion exponential factor for the vacancy kMC diffusion simulation was fitted by
computed by MD simulations and kMC simulations. Modified pre-exponential Pk¼6  MD 2
 g, where k is the index of models with
factor was applied for kMC simulation in Model_RSS and Model_AN 1 ~ 5, respectively. v0 ¼ minf k¼1 Dk  DkMC k
different chemical ordering degree (k = 1 and k = 2 ~ 6 represent Mod-
el_RSS and Model_AN 1 ~ 5, respectively). This modified pre-exponential
factor includes factors such as temperature depended vibration and
ordering caused diffusion barrier changes will lead to the vacancy trapping enthalpy effects. Based on the MD simulation results, the modified pre-
effect28. Meanwhile, DFT calculation shows that chemical ordered structure exponential factor is calculated as v0 = 4.16 × 1013 s−1.
will reduce defects jump frequencies and correlation factors16. However,
how chemical ordering structures directly affect defect diffusion behavior Figure 9 shows the MSD-time curve with modified pre-exponential
and how to control defect diffusion has not been carefully explored. factor and vacancy diffusivity in Model_RSS and Model_AN 1 ~ 5 are

npj Computational Materials | (2024)10:134 7

https://doi.org/10.1038/s41524-024-01322-6 Article

calculated as 3.6 × 10−13 m2 s−1, 2.1 × 10−13 m2 s−1, 1.8 × 10−13 m2 s−1, Chemical domain structure identifying
1.4 × 10−13 m2 s−1, 1.3 × 10−13 m2 s−1 and 1.2 × 10−13 m2 s-1, respectively. In All chemical ordering structures were identified by using the chemical
this way, the vacancy diffusivity calculated by kMC simulation shows good domain structure matching analysis. PThe matching degree of atom site i was
defined as M i d ¼ maxI ;T  N1
I NI
agreement with MD simulation results within at most 18% error in MSD. K¼1
d
δ αðKÞβðKÞ , where Id is the chemical
d Id I
The ANN predicts the vacancy jump barriers by the inputs of on-lattice ordering type (type I, II, and III),d TId is the space group operation of che-
atomic configuration of whole model, and then ANN can output 12 energy mical ordering type Id (e.g., rotating 90° around y axis for type I in Fig. 2a, so
barriers of the vacancy jump to the 12 nearest neighbor sites of FCC. The that the Cr-rich plane normal to [101] direction become to the normal [101]
energy barriers in the training dataset are computed using Nudged Elastic direction), N Id is the total number of atom sites in type Id chemical ordering
Band (NEB) method with the fully relaxed atomic configurations before and structure, αðKÞ represents element type at K site in the chemical ordering,
after the vacancy jump. The average error between the ANN predictions and and βðKÞ represents the actual element type at K site in the model. δ is
NEB calculation in different vacancy configurations is only 0.05 eV and Kronecker delta, which means that if βðK Þ ¼ αðKÞ, then δ αðKÞβðKÞ ¼ 1,
I
shows good accuracy of ANN predictions31. The P average incubation time for a otherwise δ αðKÞβðKÞ ¼ 0. When M i d >0:85, we define site i to belong to the
vacancy jump was calculated Pas Δt kmc ¼ 1= 12
i¼1 vi . The annealing time was chemical ordering type Id . This matching analysis method has been suc-
calculated as t a ¼ cExp:cvðT Þ Δt kmc , where cv ¼ N1 is the actual vacancy cessfully used to find the chemical ordering structures in our previous
concentration during annealing simulation, N ¼ 2916 is the total number of
v
work31.
atoms, cv ðT Þ is the initial vacancy concentration at temperature T deter-
Exp:

mined by experimental Time–Temperature–Electric resistance relation32. Residence density distribution and correlation with chemical
It is crucial to understand that using a unified pre-exponential factor is ordering
an approximation. The local chemical structure (local chemical ordering) Defect residence density and chemical ordering degree of each small mesh
surrounding a vacancy significantly influences the vibrational entropy region was calculated using mesh region division. The model was divided
contributing to the free energy barrier of a vacancy jump. Specifically, this into small mesh regions using a 36 × 36 × 36 mesh. The total sampling size
local chemical structure shapes the potential energy landscape and dictates of both MD and ANN-kMC simulations is 100,000. The average chemical
I
the vibrational modes at the initial and saddle point configurations of the I nd
ordering degree of each mesh region Pkd ¼ Nk was assumed to be constant
vacancy jump. Consequently, changes in the local chemical structure can k

alter the pre-exponential factor due to shifts in activation vibrational during diffusion, where k is the mesh region index, Id is the type of structure
I
entropy. (type I, II, III, and “other”), nkd is the total number of atoms of type Id around
As discussed in this paper, vacancy diffusion pathways vary according (1st nearest neighbor) mesh region k, and N k is the total number of atoms
I P I
to the chemical ordering degree of the alloy. Therefore, during the vacancy around mesh region k (0 ≤ Pkd ≤ 1, Id Pkd ¼ 1). We next defined a con-
diffusion process, the encountered local chemical structure statistics
1; rðtÞ 2 Ωk
expected to vary based on the chemical ordering degree of the alloy. To ditional function Ck ðr Þ ¼ to calculate defect residence
0; rðtÞ 2
= Ωk
address above discrepancy between KMC and MD simulations, we propose
numbers in mesh region k, where Ωk is the mesh region k, and rðtÞ is
to incorporate activation vibrational entropy into the training dataset of
interstitial or vacancy defect position at time t. When defects located in the
energy barrier and to train the ANN using this dataset. This methodological
mesh region Ωk at sampling time t, mesh region k will be counted once. In
enhancement will be further explored in our future research.
this case, the normalized total defect residence density in mesh region k was
R tF
Interstitial diffusion molecular dynamics simulation calculated as Qk ¼ a 0 C k ðrðtÞÞΔt dt, where a ¼ tnF is the normalization
P
One interstitial atom was created in the model containing 2916 atoms to factor, t F ¼ nk¼1 Qk is the total diffusion simulation time, n
perform the interstitial diffusion MD simulation. The interstitial atom was (=36 × 36 × 36) is the total number of mesh regions in the model, and Δt is
randomly inserted in the model and migrated by constantly forming new the residence time for each sampling. For interstitial diffusion simulation,
dumbbell pairs with surrounding lattice atoms. The periodic boundary the residence time for each sampling is equal to the sampling interval, which
conditions were applied along X, Y and Z directions. Interstitial diffusion is constant at 0.5 ps, while for vacancy diffusion simulation, the actual
with a total time of 1 ns was performed, and the timestep was set to 0.001 ps. residence time for each sampling should be calculated based on the kMC
The simulation temperature was set to T = 1200 K to ensure an adequate jump time Δt kmc . Finally, we defined the average correlation factor between
diffusion process within the limited simulation time window. The NPT the defect residence density and the chemical ordering degree of type Id as
P I
ensemble was adopted and temperature and pressure were controlled by GId ¼ an1 nk¼1 Pkd Qk . A higher correlation factor GId indicates that defects
35
Nosé-Hoover Pthermostat . The2 MSD was calculated as prefer to diffuse in type Id region. In this work, adequate diffusion simu-
< R ðt Þ > ¼ N i¼1 Ri ðtÞ  Ri ð0Þ , where Ri ðtÞ and Ri ð0Þ are the posi-
2 1 N 
lations were performed to capture the characteristic of interstitial and
tion vectors of atom i at time t = t and t = 0, respectively, and N is the total vacancy diffusion in different regions. See Supplementary Fig. 12 for con-
number of atoms in the model with interstitial. Interstitial diffusivity can be vergence information on residence density and number of diffusion
assessed by calculating the gradients of the MSD-time plots. Diffusion simulations.
trajectories were visualized by using the software OVITO, and interstitial
positions were identified by using Wigner-Seitz method36–38.
Chemical ordering parameter ij
ij
Vacancy diffusion ANN-kMC simulation The Warren–Cowley order parameter αn ¼ 1  PCn was used to descript the
j
Similar to the annealing simulation, ANN-kMC vacancy jump simulation chemical ordering parameters for the element pairs i and j39, where Pijn is the
was performed to study vacancy diffusion in CrCoNi, and ANN was used to probability of occurrence of element j in the nth neighbor coordinate shell of
predict the vacancy jump barrier in each kMC step. One vacancy was element i, and Cj is the concentration of element j in the system. A positive
introduced in the model and periodic boundary conditions were applied in parameter indicates a tendency of repulsion pairs of each other (e.g., Cr-Cr
X, Y and Z directions. In each vacancy diffusion simulation, 2000 accepted pair has a positive αCrCr
1 in this work), while negative parameters indicate a
vacancy jumps were performed. The vacancy diffusion temperature was set tendency of attraction with each other (e.g., the Cr-Co and Cr-Ni pairs have
as T = 1200 K. Base on the kMC vacancy jump, the MSD of vacancy dif- a negative αCrCo
1 and αCrNi
1 in this work). If the value equal to 0, it means that
fusion can be calculated, and the vacancy diffusivity can be evaluated by pairs of i and j elements are randomly distributed.
calculating the gradients of MSD-time plots. The software OVITO and
Winger-Seitz method were used to visualize simulation results and identify Data availability
vacancy position, respectively. The data used in this work are available from the authors upon request.

npj Computational Materials | (2024)10:134 8

https://doi.org/10.1038/s41524-024-01322-6 Article

Code availability 22. Inoue, K., Yoshida, S. & Tsuji, N. Direct observation of local chemical
The computer codes used in this work are available from the authors upon ordering in a few nanometer range in CoCrNi medium-entropy alloy by
request. atom probe tomography and its impact on mechanical properties.
Phys. Rev. Mater. 5, 85007 (2021).
Received: 27 December 2023; Accepted: 14 June 2024; 23. Chen, X., Yuan, F., Zhou, H. & Wu, X. Structure motif of chemical short-
range order in a medium-entropy alloy. Mater. Res. Lett. 10,
149–155 (2022).
References 24. Liu, L. et al. Local chemical ordering and its impact on radiation
1. Yeh, J.-W. et al. Nanostructured high‐entropy alloys with multiple damage behavior of multi-principal element alloys. J. Mater. Sci.
principal elements: novel alloy design concepts and outcomes. Adv. Technol. 135, 13–25 (2023).
Eng. Mater. 6, 299–303 (2004). 25. Zhang, Z. et al. Effect of local chemical order on the irradiation-
2. Cantor, B., Chang, I. T. H., Knight, P. & Vincent, A. J. B. Microstructural induced defect evolution in CrCoNi medium-entropy alloy. Proc. Natl
development in equiatomic multicomponent alloys. Mater. Sci. Eng. A Acad. Sci. 120, e2218673120 (2023).
375–377, 213–218 (2004). 26. Xu, B. et al. Influence of short-range order on diffusion in multiprincipal
3. Enserink, M. Malaria treatment: ACT two (Science (560)). Science 318, element alloys from long-time atomistic simulations. Phys. Rev.
1866 (2007). Mater. 7, 1–11 (2023).
4. Lu, C. et al. Enhancing radiation tolerance by controlling defect 27. Li, Y. et al. Chemical ordering effect on the radiation resistance of a
mobility and migration pathways in multicomponent single-phase CoNiCrFeMn high-entropy alloy. Comput. Mater. Sci. 214,
alloys. Nat. Commun. 7, 1–8 (2016). 111764 (2022).
5. Xu, Q., Guan, H. Q., Zhong, Z. H., Huang, S. S. & Zhao, J. J. Irradiation 28. Xing, B., Wang, X., Bowman, W. J. & Cao, P. Short-range order
resistance mechanism of the CoCrFeMnNi equiatomic high-entropy localizing diffusion in multi-principal element alloys. Scr. Mater. 210,
alloy. Sci. Rep. 11, 1–8 (2021). 114450 (2022).
6. Luo, H., Li, Z., Mingers, A. M. & Raabe, D. Corrosion behavior of an 29. Xing, B., Zou, W., Rupert, T. J. & Cao, P. Vacancy diffusion barrier
equiatomic CoCrFeMnNi high-entropy alloy compared with spectrum and diffusion correlation in multicomponent alloys. Acta
304 stainless steel in sulfuric acid solution. Corros. Sci. 134, Mater. 266, 119653 (2024).
131–139 (2018). 30. Yu, P., Du, J. P., Shinzato, S., Meng, F. S. & Ogata, S. Theory of
7. Schuh, B. et al. Mechanical properties, microstructure and thermal history-dependent multi-layer generalized stacking fault energy—a
stability of a nanocrystalline CoCrFeMnNi high-entropy alloy after modeling of the micro-substructure evolution kinetics in chemically
severe plastic deformation. Acta Mater. 96, 258–268 (2015). ordered medium-entropy alloys. Acta Mater. 224, 117504 (2022).
8. El-Atwani, O. et al. Outstanding radiation resistance of tungsten- 31. Du, J. P. et al. Chemical domain structure and its formation kinetics in
based high-entropy alloys. Sci. Adv. 5, eaav2002 (2019). CrCoNi medium-entropy alloy. Acta Mater. 240, 118314 (2022).
9. Pickering, E. J. et al. High-entropy alloys for advanced nuclear 32. Li, L. et al. Evolution of short-range order and its effects on the plastic
applications. Entropy 23, 98 (2021). deformation behavior of single crystals of the equiatomic Cr-Co-Ni
10. Li, Y., Li, R., Peng, Q. & Ogata, S. Reduction of dislocation, mean free medium-entropy alloy. Acta Mater. 243, 118537 (2023).
path, and migration barriers using high entropy alloy: Insights from the 33. Ding, J., Yu, Q., Asta, M. & Ritchie, R. O. Tunable stacking fault
atomistic study of irradiation damage of CoNiCrFeMn. energies by tailoring local chemical order in CrCoNi medium-entropy
Nanotechnology 31, 8 (2020). alloys. Proc. Natl Acad. Sci. USA 115, 8919–8924 (2018).
11. Derlet, P. M. & Dudarev, S. L. Microscopic structure of a heavily 34. Zhou, L. et al. Atomic-scale evidence of chemical short-range order in
irradiated material. Phys. Rev. Mater. 4, 1–20 (2020). CrCoNi medium-entropy alloy. Acta Mater. 224, 117490 (2022).
12. Wirth, B. D. How does radiation damage materials? Science 318, 35. Nosé, S.A molecular dynamics method for simulations in the
923–924 (2007). canonical ensemble A molecular dynamics method for simulations in
13. Yang, T. et al. Irradiation responses and defect behavior of single- the canonical ensemblet. Mol. Phys. 52, 255–268 (1984).
phase concentrated solid solution alloys. J. Mater. Res. 33, 36. Stukowski, A. Visualization and analysis of atomistic simulation data
3077–3091 (2018). with OVITO-the Open Visualization Tool. Model. Simul. Mater. Sci.
14. Huang, W. & Bai, X. M. Machine learning based on-the-fly kinetic Eng. 18, 015012 (2010).
Monte Carlo simulations of sluggish diffusion in Ni-Fe concentrated 37. Stukowski, A. Structure identification methods for atomistic
alloys. J. Alloy. Compd. 937, 168457 (2023). simulations of crystalline materials. Model. Simul. Mater. Sci. Eng. 20,
15. Xu, B. et al. Mechanism of sluggish diffusion under rough energy 045021 (2012).
landscape. Cell Rep. Phys. Sci. 4, 101337 (2023). 38. Stukowski, A., Bulatov, V. V. & Arsenlis, A. Automated identification
16. Zhao, S., Osetsky, Y. & Zhang, Y. Preferential diffusion in and indexing of dislocations in crystal interfaces. Model. Simul. Mater.
concentrated solid solution alloys: NiFe, NiCo and NiCoCr. Acta Sci. Eng. 20, 085007 (2012).
Mater. 128, 391–399 (2017). 39. de Fontaine, D. The number of independent pair-correlation functions
17. Lu, C. et al. Radiation-induced segregation on defect clusters in in multicomponent systems. J. Appl. Crystallogr. 4, 15–19 (1971).
single-phase concentrated solid-solution alloys. Acta Mater. 127,
98–107 (2017). Acknowledgements
18. Tsai, K. Y., Tsai, M. H. & Yeh, J. W. Sluggish diffusion in Co-Cr-Fe-Mn- Y.L., S.S., and S.O. were used computational resources of supercomputer
Ni high-entropy alloys. Acta Mater. 61, 4887–4897 (2013). Fugaku provided by the RIKEN Center for Computational Science (Project
19. Vaidya, M., Trubel, S., Murty, B. S., Wilde, G. & Divinski, S. V. Ni tracer IDs: hp230205 and hp230212), the large-scale computer systems at the
diffusion in CoCrFeNi and CoCrFeMnNi high entropy alloys. J. Alloy. Cybermedia Center, Osaka University, the Large-scale parallel computing
Compd. 688, 994–1001 (2016). server at the Center for Computational Materials Science, Institute for
20. Koyama, T., Tsukada, Y. & Abe, T. Simple approach for evaluating the Materials Research, Tohoku University, and Research Center for Compu-
possibility of sluggish diffusion in high-entropy alloys. J. Phase tational Science, Okazaki, Japan (Project: 24-IMS-C503). S.S. and S.O.
Equilib. Diffus. 43, 68–77 (2022). acknowledge the support by projects of MEXT of Japan (grant numbers
21. Zhang, R. et al. Short-range order and its impact on the CrCoNi JPMXP1122684766, JPMXP1020230325 and JPMXP1020230327). S.S.
medium-entropy alloy. Nature 581, 283–287 (2020). acknowledges the support by projects of JSPS KAKENHI (grant numbers

npj Computational Materials | (2024)10:134 9

https://doi.org/10.1038/s41524-024-01322-6 Article

JP21K14042, JP24K17170, and JP24H00994). S.O. acknowledges the Reprints and permissions information is available at
support by projects of JSPS KAKENHI (grant numbers JP17H01238, http://www.nature.com/reprints
JP21K18675, JP23H00161, and JP23K20037).
Publisher’s note Springer Nature remains neutral with regard to
Author contributions jurisdictional claims in published maps and institutional affiliations.
Yangen Li: methodology, software, formal analysis, and original draft writing;
Jun-Ping Du: software and methodology; Shuhei Shinzato: software and Open Access This article is licensed under a Creative Commons
methodology; Shigenobu Ogata: conceptualization, writing—review & Attribution 4.0 International License, which permits use, sharing,
editing, supervision, and project management. adaptation, distribution and reproduction in any medium or format, as long
as you give appropriate credit to the original author(s) and the source,
Competing interests provide a link to the Creative Commons licence, and indicate if changes
The authors declare no competing interests. were made. The images or other third party material in this article are
included in the article’s Creative Commons licence, unless indicated
Additional information otherwise in a credit line to the material. If material is not included in the
Supplementary information The online version contains article’s Creative Commons licence and your intended use is not permitted
supplementary material available at by statutory regulation or exceeds the permitted use, you will need to
https://doi.org/10.1038/s41524-024-01322-6. obtain permission directly from the copyright holder. To view a copy of this
licence, visit http://creativecommons.org/licenses/by/4.0/.
Correspondence and requests for materials should be addressed to
Shigenobu Ogata. © The Author(s) 2024

npj Computational Materials | (2024)10:134 10