Accelerated kinetic Monte Carlo: A case study; vacancy and dumbbell interstitial diffusion traps in concentrated solid solution alloys

Because of their unique mechanical and physical properties, concentrated solid solution alloys (CSAs)—including high-entropy alloys (HEAs)—have generated a substantial amount of recent scientific interest.^1–6 Notably, Ni-containing fcc CSAs have demonstrated stability at high temperature and in ion-irradiated environments,^7–10 which is promising with respect to the design of fourth generation nuclear reactor materials.^11–13 This tolerance to ion irradiation and high-temperature is thought to be associated with so-called sluggish diffusion—a concept that is addressed in the next paragraph—in these chemically complex CSAs.^14–16 One should note that tolerance to ion-irradiation does not necessarily translate to tolerance to neutron irradiation. For example, Co and—to a lesser extent—Ni have large neutron absorption cross sections,¹⁷ which leads to significant buildup of He under irradiation.¹⁸ In other words, ion-irradiation of CSAs demonstrated that radiation damage can be managed by tuning chemical composition of the alloys—in addition to tuning its micro-structure—essentially by introducing diffusion traps.¹⁴ However, since we cannot directly use these Co- or Ni-containing HEAs in reactors, a better understanding of the underlying atomistic diffusion mechanisms is necessary in order to design materials that possess similar diffusion traps, without introducing Co to the alloy.

A number of models for diffusion in alloys exist mostly in the context of dilute random^19–21 and ordered alloys.²² However, these models cannot be easily generalized to CSAs, given the great number of local environments that appear during atomic diffusion. This motivates the development of computational strategies that are able to predict atomic diffusion in CSAs, which need to be validated by a comprehensive set of experiments. Of particular interest is the concept of “sluggish diffusion” in CSAs. There are conflicting evidence and models with regards to sluggish diffusion in CSAs, as some researchers reported sluggish diffusion^23–26 and some contradict this phenomenon.^5,27,28 For instance, the conclusion of diffusion-couples studies depends on the theoretical framework employed to extract diffusion coefficients. Furthermore, there is no accepted quantitative measure of the decrease in diffusion coefficients that defines “sluggishness.”

In this context, atomistic simulations of point defects diffusing in CSAs can bring clarity, as they directly probe the particular micromechanisms of interest. In the past few years, μs-scale molecular dynamics (MD)²⁹ and kinetic Activation Relaxation Technique (k-ART)^30–34 simulations of the diffusion of vacancies and self-interstitial atoms (SIAs) were performed in Ni_xFe_1−x systems, where x ranges from 0 to 1. These studies showed that the tracer diffusion coefficients of both vacancies and SIAs vary non-monotonically as a function of composition.^35,36 Furthermore, a remarkably simple lattice-kinetic Monte Carlo (kMC) model was able to accurately capture the main trends in atomic tracer diffusion coefficients. In this simple kMC model, the transition state energies used by the kMC are obtained by averaging Ni and Fe atomic jump barriers calculated for each composition. The main results of these MD and kMC simulations are presented in Sec. IV. In other words, these studies showed that exploring this simplified static energy landscape captures the essential features of the complex, dynamic landscape explored by MD, at a much lower computational cost. Another advantage of kMC is that the transition state energies can be changed intentionally in order to study diffusion phenomena and understand the composition dependence of vacancy and SIA diffusion coefficients. The work presented in this article is based on the following idea: we artificially modify the simplified energy landscape and verify if these modifications enhance or inhibit the non-monotonicity of diffusion coefficients.

The articles mentioned above suggest that configurational traps enhance non-monotonicity—which will typically be associated with an increase in the proportion of anti-correlated jumps. Escaping these traps using thermodynamically defined probabilities can require thousands, or millions, of kMC steps. This is highly inefficient and more sophisticated techniques are needed. Accelerated kMC methods aim at characterizing such multi-state processes, sometimes known as escaping “flickers.” Accelerated algorithms can reach longer times than the conventional kMC method in a computationally efficient manner. These include parallel kMC,^37,38 the Poisson and binomial τ-leap methods,^39,40 and statistical mechanics based methods.^41–44 Each has its own drawbacks. Parallel kMC does not perform well if large time separation exists between events, while statistical mechanics based methods need a minimum time separation. Accelerated algorithms also include absorbing Markov chain kMC (AMC-kMC). This technique was brought to light by Novotny et al. almost 25 years ago. There are two main AMC-kMC implementations: the first passage time analysis (FPTA) and mean rate method (MRM). AMC-kMC works well with or without time separation and requires matrix diagonalization or inversion. These matrix operations can be a computationally prohibitive when dealing with large transition matrices (Markov matrices).^45–47 Finally, let us mention the accelerated superbasin-kMC (AS-kMC),⁴⁸ a method that modifies the rates when the system is trapped in a superbasin or—in other words—if a site or sites are visited a certain number of times. It may interrupt the dynamics of frequent events but lets the kMC execute rare events. AS-kMC is thought to be a computationally efficient method. Currently, two commonly employed adaptive kMC methods—EON⁴⁹ and k-ART³⁰—employ the MRM-KMC to escape configurational traps.

In this article, the origin of non-monotonic composition-dependence of diffusion coefficients of mono-vacancy and ⟨100⟩ dumbbell interstitial is studied. Three accelerated kMC methods are adopted and benchmarked: FPTA-kMC, MRM-kMC, and AS-kMC. Some counterfactual situations were designed and explored using conventional and accelerated kMC simulations of vacancies and SIAs in Ni_xFe_1−x to reveal the role of geometrical percolation, local defect energies, and migration barrier composition-dependent crossover.

Separate kMC simulations are performed for mono-vacancy and ⟨100⟩ dumbbell interstitial diffusion in Ni_xFe_1−x fcc CSAs. When possible, conventional-kMC is used. Accelerated kMC methods are used when conventional-kMC cannot escape configurational traps. FPTA-kMC, MRM-kMC, and AS-kMC can access longer simulation times. Note, however, that these accelerated methods share most of the main physical limitation of conventional-kMC: atomic vibrations and comprehensive chemical effects are not considered. Therefore, all kMC methods are benchmarked by MD and k-ART. Our FPTA-kMC implementation follows most of Puchala’s prescriptions, with an important exception: the final absorbing state is picked differently. Therefore, both FPTA-kMC and MRM-kMC are fully described in  Appendixes A and  B in order to clarify the necessity of this modification. In addition, detailed and step-by-step conventional-kMC algorithms for mono-vacancy and dumbbell interstitial are explained at the end of this section in order to give better insight into point defects, predefined moves, and simulation steps.

Implementing AS-kMC is fairly straightforward, starting from a conventional-kMC code. The main modification is adding a database that contains the number of times each site has been visited. Equation (A17) can be used to find the minimum number of times to visit the same state continuously before modifying the rate by Eq. (A18). The parameters of Eq. (A17) are chosen conservatively, which favors accuracy over cost-effectiveness. This accuracy is related to user-defined uncertainty value (0 < δ < 1). AS-kMC becomes more efficient as sites are visited for a couple of N_f times. If sites are visited for n time of N_f, it is shown that the error is <nδ.⁴⁸ 

In this study, the AMC-kMC is implemented in a way that mirrors our AS-kMC implementation: if any state is visited ten times during a conventional-kMC run, the AMC-kMC is “turned on.” The number of states in the Markov chains handled by the algorithm is indicated in Table V. If this number is too small, the algorithm cannot escape the traps. If it is too large, the computational cost can become problematic. Note that in MRM-kMC, there is no correlation between mean exit time and the state in which the superbasin is entered. On the contrary, FPTA-kMC is an exact method, which solves the master equation directly and does correctly capture such correlations.^45–47 Therefore, FPTA-kMC will be used to benchmark MRM-kMC and AS-kMC.

In Subsections II A 1 and II A 2, conventional, lattice-based, kMC simulations of mono-vacancy and dumbbell interstitial diffusion are explained separately. The main difference between them is in the calculation of migration energies.

• 1.

  A 10 × 10 × 10a³ (a is the lattice parameter) simulation box with atoms forming the fcc lattice structure with periodic boundary conditions is defined. Then, one atom is removed.

• 2.

  At first, all the atoms are Ni (one type). Then, Ni atoms are randomly replaced by Fe atoms until the required chemical composition of the Ni_xFe_1−x alloy is obtained. Eleven different compositions are considered.

• 3.

  A monovacancy has 12 possible jumps to the nearest neighbor sites in an fcc crystal. As a consequence, 12 paths are considered.

  • 4.1.

    The nudged elastic band (NEB)^50,51 method is used to calculate migration barriers. For each chemical composition, on average, more than 3000 configurations are considered to calculate the migration energies (E^m) of Ni and Fe jumps (Table I). A constant pre-exponential factor of 10 THz is assumed based on matching the kMC result with molecular dynamic simulation for pure Ni at 1100 K. More details of the method can be found in Ref. 52. Rates are calculated using the Arrhenius equation, as indicated by the following equation:

where ν is the jump frequency, R_i→j is the jump rate from state i to j, and E_i→j is the migration barrier between state i and j. k_b and T are the Boltzmann constant and absolute temperature, respectively.

• 4.2.

  The Bonny potential is used in this study.⁵³ This potential is chosen because it captures important energetic features of vacancies and dumbbell interstitials in equiatomic NiFe, as confirmed by electronic density functional theory (DFT) calculations.⁵⁴ For instance, vacancy tends to jump to Fe sites rather than Ni sites. Furthermore, on average, Ni–Ni dumbbells are more energetically stable than Ni–Fe dumbbells, and Ni–Fe dumbbells are more energetically stable than Fe–Fe dumbbells.

• 5.

  Conventional-kMC steps are as follows:

• 5.1.

  Calculate rates for the 12 available paths using Eq. (1).

• 5.2.

  Make a catalog by taking summation over all the 12 rates as follows:

where i is the current state of vacancy and j goes over the 12 nearest neighbor sites.

• 5.3.

  A random number r in the range of (0,1] is defined and multiplied by R_total, and pick the path as follows:

Based on Eq. (3), path j is the path that is going to be picked among the 12 available paths.

• 5.4.

  Define another random number r in the range of (0,1] and increase the time by using the following equation:

• 6.

  Steps 5.1–5.4 are repeated for at least 50 000 jumps, and the tracer diffusion coefficient, D^*, is calculated through the atomic square displacement. Effective migration energies are calculated using an Arrhenius treatment.

• 1.

  A 10 × 10 × 10a³ (a is the lattice parameter) simulation box with atoms forming the fcc lattice structure with periodic boundary conditions is defined. Additionally, an atom is added to the box to form a ⟨100⟩ dumbbell interstitial. In total, there are 4001 atoms in the box.

• 2.

  At first, all the atoms are Ni (one type). Then, Ni atoms are randomly replaced by Fe atoms until the required chemical composition of the Ni_xFe_1−x alloy is obtained. Fourteen different compositions are considered.

• 3.

  For each of the two atoms forming the ⟨100⟩ dumbbell, four possible jumps to nearest neighbor sites are considered. Therefore, there are a total of eight possible paths (2 × 4 = 8).

• 4.

  Migration energies (E^m) and jump frequencies are required to calculate rates (1). A constant pre-exponential factor of 2.5 THz that reproduced the MD result for the dumbbell interstitial diffusion coefficient in pure Ni at 500 K is used.⁵² 

  • 4.1.

    Migration energies are calculated based on the Bell–Evans–Polanyi principle. This is possible, in part, because—on average—the reaction path of a dumbbell interstitial from one site to another has a fixed distance.⁵⁵ Such an assumption does not necessarily hold for simulating more complex arrangements of defects. As indicated in Fig. 1, the change in the migration energy follows the change in the formation energy (E^f). Migration energies for Ni (or Fe) jumps in pure Ni (or Fe) systems can be calculated using static or dynamic methods. Static methods include the activation relaxation technique nouveau⁵⁶ and the nudged elastic band method. Dynamics methods include Arrhenius treatment of MD runs. In this study, values from static methods are used.⁵² The average Ni–Ni, Ni–Fe, and Fe–Fe dumbbell formation energies E^f for each composition are considered.

  • 4.2.

    If the type of the dumbbell is not changed by a jump, then the jump barrier corresponds to that of pure Ni or pure Fe (for instance, if the dumbbell is Ni–Fe before and after the jump, and it has a Ni jump, then E^m for Ni jump in pure Ni should be used). However, if the bond type is changed, such as Ni–Ni to Ni–Fe or Ni–Fe to Ni–Ni, then the migration energy is $ENim±ENi–Nif$ − $ENi–Fef$⁠. Consequently, all the eight possible migration energies [Ni–Ni → Ni–Ni, Ni–Ni → Ni–Fe, Fe–Fe → Fe–Fe, Fe–Fe → Ni–Fe, Ni–Fe → Ni–Fe (Ni jump), Ni–Fe → Ni–Fe (Fe jump), Ni–Fe → Ni–Ni, and Ni–Fe → Fe–Fe] can be calculated. All possible migration energies are collected in Table II.

• 5.

  Conventional-kMC steps are similar to those used to simulate the monovacancy except that instead of 12 paths, 8 paths are used by the dumbbell interstitial.

• 6.

  For each chemical composition, at least 100 000 jumps are simulated, and the tracer diffusion coefficient is calculated from atomic square displacement. Effective migration energies can be calculated using an Arrhenius treatment.

In kMC simulations, the defect migration energies play an important role. The computationally simple approach to calculate barriers described in conventional-kMC section reproduced the main features of the energy landscape of defect diffusing in Ni_xFe_1−x systems, as reported in Refs. 35, 36, and 52. Here, we want to separate the effect of geometrical site percolation from the effect of a migration barrier crossover, which includes the composition dependency of migration and formation energies of monovacancies and dumbbell interstitials. In order to do this, we introduce counterfactual scenarios in which the formation and migration energies follow different rules. These scenarios lead us toward not only the origin of non-monotonic behavior of diffusion coefficients but also an evaluation of the efficiency and accuracy of accelerated kMC techniques. To our knowledge, there is little direct comparison of such accelerated methods—applied to such complex energy landscapes—in the literature.

Two different vacancy scenarios are considered:

• Scenario 1:

  This is the unmodified scenario described in monovacancy subsection in conventional-kMC. The migration energies are calculated by the NEB^50,51 method and are indicated in Table I.

• Scenario 2:

  Constant migration energies for Ni and Fe jumps (migration energy of the equiatomic Ni–Fe alloy, $ENim$ = 1.305 eV, $EFem$ = 1.039 eV) are considered for all chemical compositions.

Three different dumbbell interstitial scenarios are considered:

• Scenario 1:

  This is the unmodified scenario described in ⟨100⟩ dumbbell interstitial subsection in conventional-kMC. The transition state energies are calculated by using the Bell–Evans–Polanyi principle,⁵² taking into account the dumbbell local energy and migration barrier of the particular component atom.

• Scenario 2:

  The migration energies of pure Ni and pure Fe are used at all compositions (⁠$ENim$ = 0.35 eV and $EFem$ = 0.27 eV). In other words, the formation energies of different dumbbells have no effect on this setup.

• Scenario 3:

  Similar to the first scenario, however, constant dumbbell formation energies are applied to all Ni_xFe_1−x alloys. The dumbbell formation energies of the equiatomic Ni–Fe alloy is used in this setup. These values correspond to $ENi–Nif−ENi–Fef$ = −0.44 eV and $ENi–Fef−EFe–Fef$ = −0.22 eV (see Fig. 6). The migration energies for all chemical composition are as indicated in Table III (which are obtained using the formulas presented in Table II).

If the migration energy is negative, it is replaced by a small value (0.001 eV). As an example, in Table III, −0.09 eV is replaced by 0.001 eV.

In this section, benchmarks of the kMC models are presented, followed by simulation results for each scenario separately. In this study, the definition of tracer diffusion coefficient is based on all Ni and Fe atoms as tracing elements. This can also be seen as tracing diffusion of the point defect.

Figure 2 shows tracer diffusion coefficients, D^*, as a function of chemical composition for mono-vacancy and dumbbell interstitial. The mono-vacancy data were previously presented in Ref. 35. The MD and conventional-kMC dumbbell interstitial data were previously presented in Ref. 36. We reproduce these values here for completeness and to illustrate an important point: the simplified kMC scheme is able to accurately capture the most important features of the energy landscape governing diffusion in these concentrated alloys. Indeed, the simple conventional-kMC model is a (surprisingly) reliable tool to study diffusion in CSAs, as it gives the same trend as MD and k-ART. Figure 2 also indicates that accelerated and conventional kMC simulations for dumbbell interstitial yield a comparable diffusion coefficient when traps are not present in the system.

Note that k-ART is an on-the-fly (adaptive) kMC technique. At each time step, it autonomously explores all transition states and builds a catalog. It can model the defect trajectory that approaches MD in accuracy, as it considers a large number of realistic transition states, as long as dynamic effects are negligible. In other words, while MD explores the dynamic free-energy landscape, k-ART explores the static potential energy landscape. The difference between diffusion barriers predicted by MD and k-ART can be attributed to dynamic, finite-temperate effects. In this version, k-ART uses a constant attempt frequency prefactor (or a pre-exponential factor), but quantitative dynamic information can be obtained by using prefactors extracted from harmonic transition state theory or fitting to MD data. Note that these k-ART simulations are accelerated by the basin-autoconstructing mean rate method that is based on the mean rate method (MRM).^30,35,52,57 The k-ART dumbbell interstitial diffusion simulations are presented here for the first time, for completeness. The idea is to show that the potential energy landscape drives the main features of dumbbell interstitial diffusion. More details about these k-ART simulations, and the method, will be given in a future report.

In order to calculate tracer diffusion coefficients, D^*, atomic square displacement as a function of time plots is required. The effective migration barriers can be calculated treating temperature dependence of D^* and applying the Arrhenius law. All steps to find D^* and effective migration barrier E^eff are presented for 50% Fe alloy in a mono-vacancy case in the supplementary material. For other alloy compositions and scenarios, only diffusion coefficients are reported, for sake of clarity. Effective migration barriers for all scenarios are indicated in the supplementary material.

Tracer diffusion coefficients vs Fe concentration for four different temperatures (500 K, 700 K, 900 K, and 1100 K) are shown in Fig. 3. Non-monotonic behavior is observed at all temperatures. The 500 K partial tracer diffusion coefficients are $DNi*$ and $DFe*$⁠. Only the 500 K partial tracer diffusion coefficients, $DNi*$ and $DFe*$⁠, are presented, and the total tracer diffusion coefficient, D^*, at other temperatures is merged in one plot, for brevity. The total tracer diffusion coefficient, D^*, at other temperatures is merged in one plot, for brevity. Separated and detailed plots are available in the supplementary material. The non-monotonic behavior is less pronounced at higher temperature. The minimum in diffusion coefficient is around 20%–25% Fe, which is similar to that observed in MD and k-ART simulations.³⁵ Here, the tracer diffusion coefficient is the slope of a line fitted to the square displacement as a function of time plot.

As we previously explained, constant migration energies are applied for Ni (⁠$ENim$ = 1.31 eV) and Fe (⁠$EFem$ = 1.04 eV) jumps for all chemical compositions in scenario 2. In this scenario, traps are observed in conventional-KMC for 10% and 20% Fe alloys at 500 K (check Fig. 4). Therefore, accelerated methods (AS-KMC, FPTA-KMC, and MRM-KMC) need to be applied in order to generate long enough vacancy diffusion trajectories. Additionally, accelerated methods applied to all chemical compositions at 500 K and 700 K in order to check the accuracy of these methods when there are no traps in the system. Tracer diffusion coefficients at different Fe concentrations are shown in Fig. 5. In order to simplify comparisons, only the total tracer diffusion coefficients, D^*, are shown (see the supplementary material for $DNi*$ and $DFe*$⁠). Diffusion coefficients demonstrate a monotonic behavior as the Fe concentration increases for all temperatures. In the 10% Fe, 500 K case, FPTA-KMC leads to a diffusion coefficient value bounded by the AS-KMC- and MRM-KMC-generated values. However, at 700 K, no trap observed in conventional-KMC; the diffusion coefficient from FPTA-KMC and AS-KMC is similar to those obtained using conventional-KMC, while MRM-KMC gives a little smaller diffusion coefficient in comparison to the other methods.

In order to calculate dumbbell interstitial migration energies, we first characterized distribution of their formation energies. Molecular statics are used to find ⟨100⟩ dumbbell formation energies. For each alloy composition, the averaged cohesive energy for each type of dumbbell is calculated based on more than 700 different configurations. A separate set of simulations is carried out to find the chemical potential of Ni and Fe for each specific chemical composition. The approach used to calculate chemical potential is based on the work of Piochaud et al.⁵⁸ It is explained in  Appendix B. More than 1000 different configurations are used to find the minimum potential energies. Dumbbells formation energies are reported in Fig. 6.

Tracer diffusion coefficients, D^*, $DNi*$ and $DFe*$ as a function of Fe concentration at four temperatures (500 K, 700 K, 900 K, and 1100 K) are shown in Fig. 7. Non-monotonic behavior for D^* with a minimum around 56%–60% is observed at all temperatures. This is in agreement with MD and k-ART values, as we indicated above (Fig. 2). When temperature increases, the minimum shifts to lower Fe concentration, as was observed in Ref. 52. $DFe*$ exhibit a plateau from 0% Fe up to 50%–60% Fe and a sudden increase above 80% Fe. $DNi*$ exhibit a minimum around 65% Fe and an increase at 70% Fe, followed by a decrease at higher Fe content. Partial tracer diffusion coefficients in low concentration alloys (10% Fe or 10% Ni) are not reported due to insufficient statistics (not enough solute jumps). $DNi*$ and $DFe*$ for 700 K, 900 K, and 1100 K are presented in the supplementary material.

In this scenario, constant values of migration energies (the same as in pure Ni and Fe) are applied for all the chemical compositions. D^* changes monotonically with Fe concentration, as shown in Fig. 8. $DNi*$ and $DFe*$ are presented in the supplementary material.

In this scenario, superbasins are observed at 60%, 65%, 70%, 80%, and 90% Fe concentration. Therefore, accelerated KMC methods are employed. All the plots are presented in Fig. 9 in order to facilitate comparisons. For clarity, partial diffusion coefficients (⁠$DNi*$ and $DFe*$⁠) are omitted. Diffusion coefficients are a non-monotonic function of the Fe content.

This work has two main aims. First, different accelerated kMC methods applied to modeling highly correlated defect migration in concentrated alloys are compared. Second, atomic transport by mono-vacancy and ⟨100⟩ dumbbell interstitial is investigated, in order to understand non-monotonic dependence of diffusion coefficients on alloy composition.

Differences between accelerated kMC schemes are more apparent in the case of modeling interstitial dumbbell migration. FPTA-kMC is the most accurate simulation since the master equation [Eq. (A3)] is solved directly at each AMC-kMC step. It is also the most expensive and time consuming of the methods since several matrix exponentials must be calculated. MRM-kMC is the most cost-efficient method. Table IV indicates the actual performance of accelerated methods for a deep superbasin (90% Fe at 500 K in scenario 3 of dumbbell interstitial). The results are presented in Table IV. All accelerated methods provide a significant time-boost relative to the conventional-kMC. MRM-kMC is the most efficient method in simulated physical time per cpu time ratio. However, AS-kMC is arguably the easiest of these accelerated kMC methods to implement.

For both vacancy and interstitial dumbbell migration simulations (Figs. 5 and 9), AS-kMC and MRM-kMC are overestimating and underestimating the value of tracer diffusion coefficients, respectively. We interpret this result as follows: MRM-kMC effectively considers that all the states within the superbasin have been visited according to their equilibrium occupation probability. It can overestimate the length of trajectories within the superbasin. To the opposite, AS-kMC, by raising barriers around the states near the “entrance” of the superbasin, favors an overly rapid escape. At lower temperatures, the MRM-kMC is more accurate than AS-kMC, as the number of jumps in the superbasin will tend to increase and favor an equilibrated state occupation probability. At a higher temperature, AS-kMC is more accurate than MRM-kMC, as the Boltzmann term “washes out” the errors introduced by increasing the migration barriers.

If accuracy is the priority, the FPTA-kMC method is the best choice. In the case of on-the-fly kMC, where the transition states’ catalog generation is much more computationally challenging than in lattice-based kMC, the additional computational cost of FPTA-kMC with respect to MRM-kMC will likely not add much to the overall computational cost of the simulation. In that case, using FPTA-kMC would be preferable.

In our system, only one defect is present. However, in many systems of interest where kMC and adaptive kMC are used (e.g., displacements cascade annealing), several defects will be present. In these cases, it is often necessary to simulate several local superbasins simultaneously.^59–62 Conceptually, AS-kMC barrier-modification can easily be implemented and interpreted, in comparison to managing several independent local AMCs.

Puchala et al. stated that MRM and FPTA lead to identical diffusion coefficients after a great number of steps are calculated.⁴⁷ Our benchmarking contradicts this statement. This is because our implementation is somewhat different from Puchala’s prescriptions. Our FPTA follows Puchala’s prescriptions up to the point where exit times are calculated [Eq. (A8)]. Puchala et al.⁴⁷ determined weighted rates using Eq. (5). The transient states probabilities at the exit time are multiplied by the rates toward absorbing states. In contrast, our approach involves one more matrix exponential, which considers each absorbing state separately [Eq. (A9)],

where R_i→j is the rate from transient state i to absorbing state j. Then, a random number r in the range of (0,1] is picked and r·R_total can be used to pick one path, as indicated in Eq. (3).

This difference leads to a different diffusion behavior. Consider a case with no traps. While computationally inefficient in comparison to conventional kMC, both FPTA and conventional-kMC should lead to identical results. Take a monovacancy in pure Ni at 500 K (scenario 2 of monovacancy). Let us use AMC containing 1 to 30 states. Also, let us reduce the number of visiting a state to start AMC from 10 to 0. In other words, all the time evolution is handled by AMC. In this scenario, FPTA and conventional-kMC should lead to identical results, while the MRM should diverge. The diffusion coefficient as a function of the AMC size is shown in Fig. 10. Our FPTA-kMC implementation shows a good agreement with conventional-kMC. As expected, MRM-kMC provides increasingly inaccurate results as the size of the AMC is increased. Puchala’s implementation of FPTA leads to results similar to those produced by the MRM. Applying Eq. (5) instead of calculating one more matrix exponential significantly reduces the accuracy of FPTA.

The vacancy migration barrier due to the jumps of Ni or Fe atoms was obtained by NEB in different alloys. Their composition dependence leads to a specific composition behavior of the diffusion coefficient. The minimum tracer diffusion coefficient is in ∼20%–25% Fe alloy, as can be seen in Fig. 3. This value is close to the site percolation threshold in the fcc crystal⁶³ and consistent with the previous MD, k-ART, and conventional-kMC study.³⁵ However, the non-monotonic behavior cannot be explained solely by the site percolation threshold, as has already been shown.³⁵ A simple explanation for this is presented in scenario 2 for vacancy. It should be noted that accelerated kMC methods are required to find reliable diffusion coefficients under this scenario. In scenario 2, constant migration energies (⁠$ENim$ = 1.31 eV and $EFem$ = 1.04 eV) are applied to all Ni_xFe_1-x alloys. Scenario 2 shows that diffusion coefficients increase as the concentration of the fast diffuser increases (Fe in this case since $EFem<ENim$⁠). Monotonic D^* vs Fe-concentration dependence is exhibited, as illustrated in Fig. 5. This indicates that chemically biased diffusion is a necessary ingredient to provide non-monotonic concentration-dependence in scenario 1. In other words, a migration barrier crossover must be added to geometrical percolation in order to explain non-monotonicity and the sluggish diffusion effect. One more question that should be answered is what will happen if the time separation between two events increases. This is another way to show that site percolation threshold on its own cannot explain non-monotonic concentration dependence of diffusion coefficients. In order to do so, scenario 2 is simulated again but with migration energies of the Ni₂₀Fe₈₀ alloy in which the difference between Ni and Fe migration energies is more than the Ni₅₀Fe₅₀ alloy (see Table I). The results are shown in Fig. 11, and once again, a monotonic behavior can be seen at all temperatures.

Transport properties of ⟨100⟩ dumbbell interstitial in the fcc crystal structure show a more complex behavior than that of monovacancy. Monovacancies’ formation energies depend on the chemical composition of CSAs. The formation energy of interstitial dumbbells relates to the chemical composition of both the alloy and dumbbell (Ni–Ni, Ni–Fe, and Fe–Fe in this study). In other words, vacancies in a solid solution alloy do not have a chemical type, while dumbbell interstitials do. Another difference is that vacancies tend to swap with Fe atoms,⁶⁴ while interstitial atoms tend to create the Ni–Ni configuration.^33,65,66 This has important consequences in component distribution—negative vacancy-Fe and positive interstitial-Ni transport coefficients as demonstrated in Ref. 33.

The dumbbell interstitial diffusion coefficients are a non-monotonic function of Fe concentration. This is seen in scenario 1 (Fig. 7). This non-monotonic behavior is due to both site percolation and chemical complexity, which, in conventional-kMC, is related to migration energies of pure Ni, pure Fe, and formation energies of different dumbbells. Fe is the fast diffusing element, while Ni–Ni is the most stable dumbbell up to 80% Fe concentration. As indicated in Fig. 7, the minimum in the diffusion coefficient is around 65% Fe. Previously, we showed that when artificially just one tenth of the variations in formation of energies are considered, the effect of Fe as a fast diffuser dominates, and a minimum around 20% Fe is observed in the diffusion coefficient vs Fe concentration plot, which is around the Fe site percolation threshold.³⁶ The geometrical site percolation for ⟨100⟩ dumbbell interstitial in the fcc crystal is the same as vacancy in bcc crystal, as geometrically both have eight nearest neighbor sites.

Scenarios 2 and 3 are designed to understand how local chemical complexity has an effect on dumbbell interstitial diffusion through physics that does not apply to vacancy diffusion. Both scenarios 2 and 3 have constant migration energies for all chemical compositions. All three types of dumbbell (Ni–Ni, Ni–Fe, and Fe–Fe) have the same migration energies in scenario 2. The formation energy differences between Ni–Ni, Ni–Fe, and Fe–Fe are completely neglected. In scenario 3, different types of dumbbells have different migration energies (see Table III). As indicated in Fig. 8, scenario 2 leads to monotonic behavior. In scenario 3, non-monotonic behavior is observed (Fig. 9), which is associated with the differences in formation energies between dumbbells of different chemistry. All the accelerated kMC methods indicate a minimum around 80% Fe. This is explained by the constant formation energies that keeps the Ni–Ni dumbbell as the stable dumbbell at all Fe concentrations. Consequently, given that the contribution of Fe to dumbbell diffusion is inhibited, the minimum is observed around 80% Fe, which is close to site percolation threshold (76% Fe); furthermore, the minimum is much more pronounced than in scenario 1.

The aforementioned logic of scenarios 2 and 3 and the non-monotonicity observed in scenarios 1 and 3 (Figs. 7 and 9) sheds light on the important role of formation energies in diffusion. The relative stability of the dumbbells changes as we modify the Ni to Fe ratio, which controls the shape of the non-monotonic plot of diffusion coefficient as a function of chemical composition. This complexity is absent during vacancy transport because vacancies do not have an “internal chemistry” that can change as they jump, whereas dumbbells can be of three types, Ni–Ni, Ni–Fe, and Fe–Fe, which can change as they jump. There are therefore three critical effects to take into account in order to predict tracer diffusion coefficients via dumbbell interstitial mechanism as a function of the Fe content: (1) geometrical site percolation, (2) differences in formation energies between Ni–Ni, Ni–Fe, and Fe–Fe dumbbells, and (3) the abatement of this difference in alloys containing 60% Fe or more. This abatement leads to a composition-dependent crossover of Ni and Fe migration barriers. The first two elements were well explained in Ref. 36. The third element explains why the minimum is not observed at the Ni site percolation threshold (about 20%–25% Ni), rather around 65%. It also explains why the concentration at which dumbbell diffusion is minimal changes as a function of temperature. Exploring the counterfactual scenarios enabled this important insight, which the accelerated kMC made possible.

As it is already mentioned, by estimating formation and migration energies, one can predict diffusion of vacancies and interstitial atoms formed by irradiation. However, radiation effects are more than solely point defect transport and involve defect clusters formation and diffusion, He production/management, radiation induced segregation, stability of micro- and nano-structural elements, and their interplay as the temperature is modified.^67–69 As a consequence, more elements and multi-scale modeling should be considered in order to predict overall radiation damage evolution.

Counterfactual mono-vacancy and ⟨100⟩ dumbbell interstitial atom diffusion scenarios were considered. The diffusion of these defects in Ni_xFe_1−x CSAs was simulated by conventional and accelerated kMC, which provided insights into the non-monotonic concentration-dependence of tracer diffusion coefficients. A detailed and systematic comparison of modeling efficiency and accuracy has been made between different kMC methods. A summary of the most important results is as follows:

• 1.

  Implementation of accelerated kMC methods for modeling point defects diffusion in a fcc crystal structure is explained step-by-step and in details. Following these prescriptions, FPTA-kMC and MRM-kMC lead to different kinetics, even in the long run, a contrario to the work of Puchala et al.

• 2.

  Accelerated kMC methods (AS-kMC, FPTA-kMC, and MRM-kMC) are benchmarked. FPTA-kMC is the most reliable method but can be quite slow when handling large superbasins. MRM-kMC and AS-kMC are faster than FPTA-kMC and can capture all important trends. Generally, AMC-kMC has a larger cost per kMC step than AS-kMC but requires less kMC steps to reach a given simulated time. All of these accelerated methods can be implemented in adaptive or on-the-fly kMC techniques such as k-ART. MRM-kMC is likely the best candidate for systems with large superbasins and a high computational overhead per kMC step, while AS-kMC is more straightforward to implement and interpret, especially in the presence of simultaneous local superbasins.

• 3.

  Non-monotonic concentration dependence of tracer diffusion coefficients cannot be explained solely by geometrical site percolation. Composition- and environment-dependent migration barriers need to be taken into account. This dependence has different origins whether diffusion occurs by mono-vacancy or dumbbell interstitial mechanism.

• 4.

  Vacancies do not have internal chemistry, and their non-monotonic diffusion is defined by composition-dependent barriers for their exchange with different alloying components. If the migration barriers did not vary as a function of Fe concentration, vacancy diffusion would be a monotonic function of Fe concentration.

• 5.

  Dumbbell interstitial formation energies depend on the dumbbell chemical composition. This is an important driver that defines a non-monotonic behavior. Even if Ni–Ni, Ni–Fe, and Fe–Fe dumbbell formation energies were independent of the CSA’s Fe concentration, non-monotonicity would be observed. In fact, in this model system, the abatement of differences in Ni–Ni, Ni–Fe, and Fe–Fe dumbbell formation energies made this non-monotonic behavior less pronounced and shifted the minimum to lower Fe concentrations.

See the supplementary material [available online within a portable document format (pdf) file] for detailed data and analyses that underlie our calculation of tracer diffusion coefficients and effective activation barriers and figures illustrating the partial tracer diffusion coefficients and effective activation barriers observed in all scenarios that were simulated.

The data that support the findings of this study are available from the corresponding author upon reasonable request. The kMC codes that support the findings of this study are available upon request from Keyvan Ferasat (17kf7@queensu.ca).

Y.N.O. and Y.Z. were supported as part of the Energy Dissipation to Defect Evolution (EDDE), an Energy Frontier Research Center funded by the U.S. Department of Energy, Office of Science, Basic Energy Sciences, under Contract No. DE-AC05-00OR22725. K.F., L.K.B, and Z.Y. were financially supported by NSERC-UNENE. K.F., L.K.B., and Z.Y. thank Compute Canada for generous allocation of computer resources.

“Flickers” can be observed in the transportation of point defects in CSAs. In other words, there are traps caused by the deep local basin in the energy landscape, which slow diffusion. In some simulation scenarios, conventional-kMC cannot escape these traps. Accelerated kMC methods are required to study these cases. FPTA-kMC, MRM-kMC,^45–47 and AS-kMC are implemented. These methods are explained in the following.

The first passage time analysis (FPTA) and mean rate method (MRM) are subclasses of absorbing Markov chain kMC (AMC-KMC). In the AMC-KMC method, the repeated states are the transient states, and states out of this basin are the absorbing states. The transition probabilities among these states are used to make the absorbing Markov matrix, as written in Eq. (A1).⁴⁵ In the following, there are q absorbing states and s transient states,

where I, 0, R, and T are submatrices of the AMC matrix M. I is the identity matrix and 0 is a zero matrix, which indicates that the probability of leaving an absorbing state is zero. T (transient submatrix) is the transition probabilities among the transient states, and R is the transition probabilities from transient states to absorbing states.

In this method, the time evolution of the system is calculated by solving the master equation for the state probabilities, as indicated in Eq. (A2).⁴⁷ Rates in this equation are calculated using Eq. (1),

where P_i is the probability of being in state i. The master equation can be rewritten in a vectorized format as follows:

where $P¯$ is probability vector and M is the absorbing Markov chain matrix, but instead of transition probabilities, rates are used to form the AMC matrix as follows:

The solution for master Eq. (A3) can be written as

where $P¯(t)$ is the probability vector at time t. $P¯(0)$ is the initial probability vector, and all the terms are zero except the last transient state, which is one. It is beneficial to consider all the absorbing states as one absorbing state, since it reduces the size of matrix M, which lessens the cost of obtaining the exponential of that matrix. $P¯(0)$ is

The last row of Eq. (A6) is for all the absorbing states, which is zero, and the other rows are transient states.

Equations (A4) and (A5) are additional steps that should be used to change the conventional-KMC to FPTA-KMC. The steps in the FPTA-KMC are as follows (here, it is just explained for the vacancy case):

• 1.

  FPTA-KMC runs as a conventional-KMC when the system is not inside a superbasin—i.e., a trap. There are 12 possible jumps, and Eq. (1) can be used to find rates. The catalog is made up of these rates, and one of the paths can be picked.

• 2.

  The FPTA is applied when the system is inside a superbasin. Here, it is assumed that if a state visited more than ten times, FPTA-KMC can be applied. Experience shows that applying ten steps of FPTA is enough to let the vacancy escape from the superbasin. The log of the number of times a state is visited is reset to zero after applying ten steps of FPTA or every 2000 steps, whichever one happens first.

• 3.

  When a state is visited ten times, the FPTA steps apply as follows:

  • 3.1.

    Add the absorbing state to the transient state for the next step.

  • 3.2.

    Build the rate matrix M, as explained in Eq. (A4).

  • 3.3.

    Consider all the absorbing states as one absorbing state.

  • 3.4.

    Define an initial value for the time in Eq. (A5). In this study, the summation of the diagonal of matrix M [Eq. (A4)] was inversed and used as an initial value for time (t₀).

  • 3.5.

    A random number (r) is defined in the range of (0,1]. This random number is the probability of being in the absorbing state after time t.

  • 3.6.

    Initial time (t₀) plus the bisection method is used to find the time (t_exit) that satisfies the probability of being in the absorbing state is equal to defined random number r, as indicated in Eq. (A7). Eigen package⁷⁰ is used to calculate the exponential of the matrix,

where t_max and t_min are the boundaries in the bisection method and t_exit is the average of these two values.

• 3.7.

  Step 3.6 is repeated until the approximation in exit time t satisfies the following condition:⁴⁷ 

• 3.8.

  Up to this point, the exit time t_exit is calculated. Now, one of the absorbing states should be picked with another random number (r′) in the range of (0,1].

• 3.9.

  In this step, the probability vector at t_exit is recalculated again, but instead of considering all absorbing states as one state, they are considered separately [Eq. (A9)]. In other words, the size of M and $P¯(0)$ matrices are increased. One could initially consider absorbing states separately and find t_exit. However, it is more efficient to get the exponential of a smaller matrix,

• 3.10.

  The FPTA-KMC catalog is built using the calculated probabilities of absorbing states at exit time (t_exit) from Eq. (A9). Then, the defined random number r′ can be used to pick one of the absorbing states.

• 3.11.

  Finally, the time in KMC simulation increases by t_exit and one of the absorbing states have been picked.

• 3.12.

  In this study, steps 3.1–3.10 are repeated for 10 times.

• 4.

  The simulation proceeds as the conventional-KMC until another trap is encountered.

It is important to mention that the minimum number of repeated AMC steps is related to the size of the superbasin. Generally, for the mono-vacancy and the dumbbell interstitial, ten steps of AMC are considered. In some cases, simulating the dumbbell interstitial diffusion required the AMC containing 30, 50, 100, and 140 states (see Table V). The more states in the AMC, the slower the simulation will be, because of the large cost of the matrix exponential calculation.

In FPTA-KMC, the exact exit time is calculated. Alternatively, one can calculate the mean exit time. This is what is done in MRM-KMC.⁴⁷ The main merit of MRM is that it is computationally more affordable than FPTA, since just one inverse matrix needs to be calculated, instead of multiple matrix exponentials. Additionally, instead of dealing with the whole AMC matrix M [Eq. (A1)], linear algebra is limited to the transient submatrix T (transition probability matrix) in AMC matrix M. In order to find the mean exit time (⟨t_exit⟩), first, the mean residence time (τ) for each transient state is calculated. The mean exit time is the summation over mean residence time of transient states. In the following, the way to build the transition matrix T and calculate the mean exit time is explained. Also, the MRM steps are explained in the context of vacancy diffusion.

The transition probability matrix T is the transition probabilities among transient states. It can be built as indicated by the following equation:

where i and j are just for transient states, while k goes over both transient and absorbing states. $τi1$ is the mean residence time in transient state i.

As the objects (mono-vacancy and dumbbell interstitial) move among the transient states, the occupation probability for transient states changes. The updated occupation probability after m jumps $P¯(m)$ can be found by applying m-times of transition probability matrix T to the initial occupation probability vector $P¯(0)$⁠, as indicated by the following equations:

Then, all the possible jumps (m) from zero to infinity are considered. This is equivalent to assuming that the system is in equilibrium. Their summation is calculated as follows:

I is the identity matrix. $P¯total$ can be used to find the mean residence time after applying all possible number of jumps (τ_i), as indicated by the following equation:

Finally, τ_i can be used to find the mean rate of exiting to absorbing states and the mean exit time (⟨t_exit⟩) as follows:

where i and k go over all the transient sates and j is for the absorbing states related to transient state i. Calculated mean rates from Eq. (A15) should be used to make the catalog for KMC, and ⟨t_exit⟩ from Eq. (A16) is the time increment for this KMC step.

The MRM-KMC is applied as follows for the monovacancy:

• 1.

  Run conventional-KMC when the system is out of a superbasin.

• 2.

  Activate MRM when the system is inside a superbasin. If a state is visited more than ten times, MRM-KMC is used. The log of the number of times a state is visited is reset to zero after applying ten steps of FPTA or every 2000 steps, whichever one happens first.

• 3.

  When a state is visited ten times, MRM steps are applied as follows:

  • 3.1.

    Add the picked absorbing state to the transient state for the next step.

  • 3.2.

    The transition probability matrix among transient states is made by Eq. (A10). Rates can be calculated by using Eq. (1).

  • 3.3.

    Equations (A13) and (A14) should be used to find the updated residence times for transient states.

  • 3.4.

    Equation (A15) should be used to find the necessary rates to make the KMC catalog. Then, a random number r in the range of (0,1] is defined, and one absorbing state will be picked from the catalog using that random number.

  • 3.5.

    Equation (A16) gives the time increment of the KMC steps.

  • 3.6.

    Then, step 3.1 to 3.5 are repeated ten times.

• 4.

  Conventional-KMC is run until another trap is detected.

Like FPTA-KMC, in some cases, for dumbbell interstitial simulation, instead of 10 AMC steps, 30, 50, 100, and 140 steps are considered (see Table V).

It is interesting to point out that the MRM-kMC is currently used to accelerate different adaptive—i.e., on-the-fly–kMC algorithms and codes.^{30,49,71–75}

Another method to escape superbasin is proposed by Chatterjee and Voter.⁴⁸ The idea is to increase the activation barrier—i.e., reduce the rate—of the “flickering” processes. In order to use this method, a database of processes is required, as well as two input parameters: uncertainty δ and the rate lowering factor α. If a process (vacancy or dumbbell interstitial jump in this study) is observed more than a specific number of times (N_f), then the rate for this process is decreased. N_f is specified by Eq. (A17). In this study, δ and α are chosen to be 0.1 and 2, respectively. Hence, the minimum value of N_f is 24,

The AS-kMC code is based on conventional-kMC, augmented by a database of processes. Two arrays, one temporary and one permanent, are kept. If a process happens more than N_f times in the temporary array, then in the permanent array, the number of repetitions for that process is increased by one, and in the temporary array, the number of repetitions is reset to zero. Additionally, the temporary array is reset every 1000 steps. Temporary arrays of size up to 10 000 are tested and showed that 1000 steps are sufficient. The permanent array is used to build a catalog with modified rates. As an example, in the case of mono-vacancy diffusion, there are 12 paths in the catalog. For each path, the number of repetition in the permanent array is checked, and Eq. (A18) is used to modify the rate,

where $Ri→j′$ is the modified rate and R_i→j is the original rate obtained by applying migration energies in Arrhenius Eq. (1). These rates are calculated once in the code and will be used to make catalog. n is the number of repetitions in the permanent array. A more efficient approach is to make a database for processes modified rates, and in Eq. (A18), instead of dividing rates by αⁿ, a modified rate from the database can be divided by just α. The error that is produced by AS-kMC depends on the uncertainty (δ) value and on how many times the rate is modified [n in Eq. (A17)].⁴⁸ 

The summary of the steps is as follows:

• 1.

  Find N_f by using Eq. (A17).

• 2.

  Find all the paths that are connected to the current state. In the case of monovacancy and ⟨100⟩ dumbbell, there are twelve and eight paths, respectively.

• 3.

  Keep track of processes and states, and make the database in the permanent array.

• 4.

  Modify the rates using Eq. (A18) if the jump leads the system toward a repeated state.

• 5.

  Pick an event and advance the clock using standard kMC rules.

The chemical potential of Ni and Fe can be calculated by solving Eqs. (B1) and (B2). This method is based on the work of Piochaud et al.⁵⁸ In principle, the chemical potential is a function of temperature and can be obtained through an ensemble average of the potential energy change involved in swapping one species for another U^Fe→Ni and U^Ni→Fe (i.e., a Monte Carlo simulation). Here, we use an approximation of this ensemble average at 0 K. Indeed, when the temperature goes toward 0 K, Δ^Fe→Niμ can be estimated by the minimum of the distribution of U^Fe→Ni,⁵⁸ 

where N_Fe and N_Ni are the number of Ni and Fe in the simulation box and Δ^Fe→Niμ is the change in chemical potential when one Fe atom is substituted by a Ni atom, and vice versa for Δ^Ni→Feμ. U₀ is the potential energy before changing atoms.

Another approximation is to use an arithmetic average of this distribution instead of an ensemble average, as explained by Zhao et al.⁵⁴ In our work, the approach of Piochaud et al. is used to calculate dumbbells’ formation energies. For the sake of comparison, the chemical potential as calculated by both approaches is shown in Fig. 12.