© G. Adjanor et al., Published by EDP Sciences, 2025

1. Introduction

Over the course of their operational lifespan, reactor components materials experience significant alterations in their physical properties because of the combined effects of neutron irradiation and operating temperatures. The reactor pressure vessel (RPV) is often considered as the most critical component as it is the only component that cannot be changed. Other critical components, like those forming the internal structures of the reactor, are potentially replaceable but they have other issues since they are exposed too much larger doses. In all cases, the evolution of their properties originates from neutron irradiation which results in a production of large amounts of crystalline point defects such as vacancies and self-interstitial atoms after the displacement cascades processes. Throughout the various timescales of cascades’ evolution (thermal spike, intra-cascade and then inter-cascade reactions), these defects migrate, recombine or agglomerate to form vacancy clusters and loops of interstitials. This can cause a variety of radiation induced segregation type phenomena which in turn modify mechanical properties by radiation hardening or embrittlement. Usually, the evaluation of RPV mechanical integrity relies on the extrapolation of semi-empirical models calibrated against data obtained from surveillance programs. In the context of extending plant operational lifetimes, computational simulations offer a promising complement to empirical extrapolations, providing a physically based understanding of the coupled effects of materials’ chemistry and incident particle fluxes [1]. The two main simulations methods which are atomistic kinetic Monte Carlo (kMC) methods and cluster dynamics (CD) have almost opposite strengths and weaknesses: the former rely on a model of the crystalline lattice and thus readily model parameterized atom-defect reactions, while the latter consist in efficiently solving with large time increments systems of ordinary differential equations associated to the mean field interactions of defect clusters classes.

Intrinsically, kMC being eventdrawing based, its efficiency is limited by the ratio between the highest and the lowest event frequencies in the system. Indeed, predicting the formation and evolution kinetics of nanometric clusters of crystalline defects which form in irradiated metallic components requires covering considerable timescales: from the picosecond timescale of the Brownian motion of these defects on the crystal lattice, up to the timescale of decades for the evolution of the mechanical properties that they influence. The usual so-called “kinetic” Monte Carlo methods, which directly simulate these complex random walks (which can be one-dimensional or three-dimensional randoms walk in the 3D crystal lattice, or more complex combinations of the latter [2]), offer a great capacity for modelling elementary atomic mechanisms, but are limited to short cumulated times (being single event-based methods). Nevertheless, they enable us to straightforwardly implement the energy models underlying the interactions between clusters of defects and solutes (with either density functional theory methods or classical molecular dynamics calculations [3]) and potentially take into account the vast and complex combinatorics associated with them. One typical application of this method is simulating the formation of solute atoms atmospheres in RPV steels, which leads to radiation embrittlement. Being a low-alloy steel (the major solutes involved in aging being Mn, Ni, Si and Cu), this kind of dilute system considerably reduces the computational burden for building up the atomic jumps parameterization. This leaves us with the challenges of the very long-term simulations. Typically, the very long timescales to be simulated can exceed 1.3 × 10⁹ s (corresponding to about 40 years of reactor operation), while elementary events such as interstitial jumps only increment the physical clock by steps of 10⁻¹² s. In fact, the computational efficiency is severely impacted by the overwhelming number of non-reactive diffusion events (i.e., which leave the microstructure intact). These transitions correspond to very low energy barriers (also called “flickers”) which monopolize the computational power. This is in contrast with rare-events where solute clusters actually evolve. In accelerated Monte Carlo methods, a randomization procedure based on the factorization of absorbing Markov matrices enables the direct drawing of an exit location from the trap, stochastic numbers of visits to each trapping site and physical times spent waiting at these sites. The randomization procedure guarantees that the total waiting time within the trap is distributed according to the correct first-passage time distributions. Each possible exit location corresponds to a dissociation site between the atomic cluster and the migrating defect. The large sequences of non-reactive events are directly drawn. The latter methods are exact and allow for accelerations of many orders of magnitude, as presented at Section 2, but they are currently limited to one type of solute atoms, and modest densities of defects.

Beyond that, very long-term kinetics is within the sole reach of mean field methods such as “cluster dynamics”, but at the cost of drastic simplifications on the treatment of complex clusters: homogenization, neglecting of spatial correlations, coarse-graining and use of steady state solutions of diffusion equations for absorption rates. The systems of ordinary differential equations associated with these methods are integrated using deterministic time-steppers [4] or using stochastic methods (Langevin dynamics or random walks in the space of cluster sizes). More recently, a hybrid approach combining deterministic and stochastic integration has also been proposed [5]. An alternative strategy [6] consists in developing new algorithms for a full stochastic resolution of rate-equations. Indeed, we show at Section 3 thanks to the concept of partial propensities [7–9], of internal times and the use of binary heap to limit reaction lists updates [10, 11]. As shown at Section 3, many orders of magnitudes speed-up are obtained compared to the standard stochastic algorithm (SSA) while maintaining the capacity to mitigate the computational burden linked to the multidimensional clusters composition space.

Finally, for some prospective applications such as high entropy alloys and other concentrated alloys, variants such as kinetic Monte Carlo on a relaxed lattice with on-the-fly prefactor calculations account for the vibrational degrees of freedom and the stabilizing distortion of the lattice. This allows to simulate the physics of slow diffusing defects in complex concentrated alloys with greater details than kMC on a rigid lattice.

2. First passage kinetic Monte-Carlo methods

To briefly sketch first passage analysis to kinetic Monte Carlo, we must first recall that the dynamics of any Markovian system can be described by a master equation:

$$\begin{array}{c}\hfill \dot{p_t^T}=p_t^{T}K,\end{array}$$(1)

where the probability vector $p_t^T$ represents the probability distribution at time t over the discrete set of states and K is the transition rate matrix. Its off-diagonal elements K_ij encode the transition rates from state i to state j and its diagonal elements satisfy the relation K_ii = −∑_j ≠ iK_ij so that the total probability is conserved, i.e., the elements of the probability vector always add up to 1. The evolution operator P from time t to time t+τ can be formalized as:

$$\begin{array}{c}\hfill P(t,t+\tau )=exp\left(\tau K\right).\end{array}$$(2)

The element P_ij(t,t+τ) corresponds to the probability of observing the system in state j at time τ + t knowing that it was in state i at time t. The operator exponential is defined by its Taylor expansion. The n-th order term of the expansion yields the contribution arising from n elemental transitions (moves). In the standard kinetic Monte Carlo (kMC), the residence time algorithm performs a single move at each step based on the first term of the Taylor expansion (see [12] for instance). The associated transition probability matrix then writes

$$\begin{array}{c}\hfill \tilde{P}=I+K{\tau }_{\max },\end{array}$$(3)

with ${\tau }_{\max }=-K_{\mathit{ii}}^{-1}$, I being the identity matrix. The physical time elapsed before the elemental transition occurs is drawn from the exponential distribution Δt ∝ exp(−Δt/τ_max). Simulating the physical time using an auxiliary random variable guarantees the correctness of the algorithm.

In practice, the efficiency of the kinetic Monte Carlo method is drastically reduced whenever the transition rate matrix describing the evolution of the system exhibits a wide spectrum. In these situations, the system transitions a huge number of times between configurations separated by small energy barriers, as illustrated in Figure 1.

Fig. 1. (a) A cluster of B solute atoms within a bulk of A atoms is sketched. In insets (b,c,d) a typical “flicker” is illustrated. It results in having the vacancy (V) returning at its initial position many times. This is due to B–V exchanges being far more likely than A-V ones (more complex flickers are possible (e)) as reflected by their respective transition rates Γ.

These connected configurations form trapping basins from which the average escape time is much larger than the characteristic time for crossing the small barriers inside the basins. This issue is recurrent in kMC simulations [13]. Cavities may form under irradiation and remain stable over a long period of time due to the low vacancy emission rate resulting from the strong attraction between cavities and neighboring vacancies. Kinetic trapping may also be caused by the formation of dynamically stable clusters of manganese or copper substitutional atoms and vacancies in α-Fe. These solute clusters migrate slowly without dissociating owing to ineffective numerous atomic rearrangements, or “flickers”. A second situation causing a drastic loss of efficiency of the kMC method is when the distances a mobile defect must travel before interacting with another mobile defect (e.g. a dislocation, a solute cluster or precipitate) is large. It takes a huge number of kMC hops to observe a small change of the microstructure.

The so-called “first-passage” algorithms allow accelerating the kMC simulations in both aforementioned situations. To achieve this goal, nonlocal events involving mobile defects may be randomly selected using exact probability rules [12, 14, 15]. Avoiding conflicts between defects evolving in parallel requires spatial protection of defects and exact time synchronization. Spatial protection serves to prevent distant walkers from colliding or conversely to enable two neighbouring defects to recombine. To satisfy the time synchronization requirement, the theory of absorbing Markov chains [16] is used to draw first-passage times, and paths to distant states located on the periphery of the protection, which acts as an absorbing sink.

The transition rate matrix K ^a associated with the absorbing Markov process is defined as:

$$\begin{array}{c}\hfill K^a=\left(\begin{array}{cc}-A\hfill & A\overrightarrow{1}\hfill \\ 0^T\hfill & 0\hfill \end{array}\right).\end{array}$$(4)

Here, matrix A is a symmetric, positive-definite matrix of size N (the number of transient states), and its elements are minus the transition rates between transient states, while $\overrightarrow{1}$ is a vector of ones (see Fig. 2).

Fig. 2. The trapping sites within the solute cluster (shown in blue) determine the N transient states, while the perimeter states determine the absorbing states. The former are associated with the A blocks of the matrix and the latter with the null blocks. At this stage, all perimeter sites are treated as a single absorbing state with an index of N + 1.

The evolution operator of the absorbing Markov chain is obtained through exponentiating the associated transition rate matrix:

$$\begin{array}{c}\hfill exp\left(K^{a}t\right)=\left(\begin{array}{ccc}exp[-At]\hfill & \hfill & (I-exp[-At])\overrightarrow{1}\hfill \\ 0^T\hfill & \hfill & 1\hfill \end{array}\right).\end{array}$$(5)

The probability vector restricted to the transient states then evolves according to $p_t^T=p_0^T\exp (-At)$, while the survival probability at time t, defined as the probability that the system remains in one of the transient trapping states is $p_t^T\overrightarrow{1}$. Then, the mean survival time or mean first passage time (MFPT) is defined as

$$\begin{array}{c}\hfill {⟨t⟩}_{\mathit{FP}}={\int }_0^{\infty }{p}_t^T\overrightarrow{1}\mathrm{d}t=p_0^T{\int }_0^{\infty }\exp [-At]\overrightarrow{1}\mathrm{d}t=p_0^T{A}^{-1}\overrightarrow{1}.\end{array}$$(6)

The above linear system can be solved by using either a direct solver, such as Gaussian elimination, or an iterative solver, such as the conjugate gradient method. First-passage time and escape from the exact probability distributions can be sampled using randomization based on direct factorization of the absorbing transition rate matrix [12], or on its eigenvalue decomposition [13–15]. The former randomization technique is based on the probabilistic interpretation of the factorization in terms of paths [17] and is referred to as kinetic path sampling (kPS). The application to the long-term evolution of the FeCu system was proposed in [12]. A parametrization based on electronic structure calculations by Soisson et al. [18] was used to model a system with a single vacancy and 1.34 at% Cu.

In this work [17], the first passage time formalism is used to perform nonlocal transitions of the vacancies and their escape from the trapping basins of the copper rich clusters. The physical time is increased by a first passage time associated with this escaping transition. Figure 3 illustrates the considerable speed-ups that these non-local techniques allow for. From high to low temperatures, speed-ups (from standard kMC) from 10² up to 10⁸ can be respectively achieved. Higher speed-ups at low temperature simply correspond to the situation in which the regular kMC algorithm is most penalized by trapping.

Fig. 3. Integrated speed-up plotted as a function of the physical time simulated by kPS at three different temperatures T0 = 273 K, T1 = 373 K, and T2 = 473 K [12].

Other applications of kPS further exploit the ability of the direct factorization to perform matrix inversion, as well as the fact that Green’s functions, which are pseudoinverses of transition rate matrices, can also be interpreted as geometric sums of path probabilities and used in atomic transport theory. This connection enables atomic diffusion coefficients to be computed very efficiently using kPS in conjunction with mass transport theory [19]. This opens the way to the efficient development of parameterization that accounts for the mobility of clusters, as demonstrated for vacancy-trapping Mn clusters in iron [13]. Figure 4 shows the extent of mobilities and the various curves for a wide range of temperatures, demonstrating that the efficiency of kPS is insensitive to temperature.

Fig. 4. Diffusion coefficients of Mn solute atoms as a function of cluster size for various temperatures T. For clarity, scaling is done relative to monomer diffusivity at T0 and its high-temperature activation energy Ea = 0.646 eV. Reference temperature is T0 = 600 K and reference diffusivity is D(1,T0)=6.716 × 10−15 m2 s−1 at vacancy concentration of 0.05% [13].

Its limitations stem from the computational overhead associated with the transition rate matrix manipulations, which involves O(n³) matrix operations and substantial RAM usage. To efficiently apply kPS to the largest clusters in this study, several strategies were employed, including storing the matrices and labelling clusters based on a connectivity graph. Indeed, analysis shows that whenever trapping becomes severe, the simulation often repeatedly visits a limited number of basins and performs the same factorizations many times. Physically, at lower temperatures, the system gets trapped in super basins containing several shapes that are repeatedly adopted by the cluster. Consequently, we proposed a first acceleration method involving storing and recycling the transition matrix factorizations, using a hash key and a dictionary to index path factorization matrices. An additional method consists in constructing a graph on-the-fly to connect cluster configurations, detect super basins and re-use the factorized transition matrices.

As Figure 5 shows, these strategies enabled acceleration by approximately two orders of magnitude compared to regular kPS and by around five orders of magnitude compared to standard kMC. This demonstrates that the computational overhead associated with kPS on large clusters can be mitigated.

Fig. 5. Comparison of algorithm efficiencies as a function of cluster size. Efficiency is measured by the ratio of simulated physical time to that of F-KMC (kPS), based on a fixed wall-clock time. Shaded regions around the curves represent the 68% confidence interval.

3. Stochastic resolution of rate-equations for cluster dynamics

In the mean-field formulation of cluster dynamics (CD) [4], one writes a system of coupled balance equations relating the concentrations C_(n,p) of clusters made of |n| single defects (where negative n denotes a vacancy cluster and positive n describes an interstitial cluster) and p solute atoms:

$$\begin{array}{c}\hfill \frac{\mathrm{d}{C}_{(n,p)}}{\mathrm{d}t}={\varphi }_{(n,p)}(t,\overrightarrow{C}),\end{array}$$(7)

where $\overrightarrow{C}$ is the concentration vector of all C_(n,p) concentrations. The right-hand side term ${\varphi }_{(n,p)}(t,\overrightarrow{C})$ is theenumeration of all reactions for affecting class (n, p), which reads (for immobile clusters):

$$\begin{array}{cc}& {\varphi }_{(n,p)}(t,\overrightarrow{C})=G_{(n,p)}+\sum _{j=-m_v}^{m_i}\sum _{q=0}^{m_s}{A}_{(n+j,p+q),(j,q)}{C}_{(n+j,p+q)}\hfill \\ \hfill & -A_{(n,p),(j,q)}{C}_{(n,p)}+\sum {}_{j=-m_v}^{m_i}\sum {}_{q=0}^{m_s}{B}_{(n-j,p-q),(j,q)}\hfill \\ \hfill & C_{(n-j,p-q)}{C}_{(j,q)}-B_{(n,p),(j,q)}{C}_{(n,p)}{C}_{(j,q)},\hfill \end{array}$$(8)

where m_i, m_v and m_s are respectively the maximum number of interstitials, vacancies and solutes in mobile clusters. The absorption rates and emission rates of a mobile cluster μ = (j, q) by a cluster (n,p) are respectively noted B_{(n,p), (j,q)} and A_{(n, p)(j, q)}.The source term G_(n, p) corresponds to instantaneous production of clusters due to irradiation and the resulting cascade process. Absorption terms (sink strengths) for mobile species due to fixed sinks like grain boundaries or dislocation networks may be included, as well as many additional physical refinements such as geometric and elastic factors related to clusters shapes [20, 21], the effect of neighbouring clusters [22] or the 1D/3D character of mobile’s clusters diffusion [23, 24] which are physically relevant but are out of this papers’ scope.

Formulating cluster dynamics in this mean-field manner enables very long times to be simulated, even when the non-linear ordinary differential equation system is stiff [4]. Traditionally, this initial value problem is solved using a multistep implicit formulation and Newton’s method to solve the associated nonlinear system. In this approach, repeated solving of a huge linear system is the most computationally intensive part. This drastically limits the extent of the cluster size space for the deterministic resolution of the CD system. Several improvements are possible, such as reformulating n as a continuous variable (yielding a Fokker-Planck equation) [4, 25, 26], grouping of equations [27, 28], or applying various linear algebra strategies for iterative solvers [6]. These can significantly increase the resolution domain, but the horizon of deterministic resolution of CD is currently limited to two-dimensional size spaces. This means that when treating mixed clusters in irradiated alloys, cluster space is limited to (n, p), for instance a cluster consisting of |n| vacancies and p copper atoms surrounding the vacancy cluster. Describing a more complex cluster in which vacancies are also associated with nickel atoms (i.e. considering a three-dimensional cluster space for (n,p,r) clusters, r being the number of nickel atoms associated with the copper-vacancy cluster, for example) is beyond the scope of deterministic simulations in conditions in which these clusters could be experimentally observed. This issue is often referred to as the dimensionality curse.

Building on the previous work of Bulatov et al. [29] and Gherardi et al. [30], who attempted to adapt Lanore and Gillespie’s original stochastic simulation algorithm (SSA) [31–33] of chemical rate equations to CD, we will now consider the stochastic resolution of the CD system. This method involves propagating one or more stochastic trajectories in cluster size space [32, 33]. The great advantage of this stochastic approach is that it is almost immune to the dimensionality curse as it is a Monte Carlo method. This means that one can simulate dynamical systems with very large numbers of degrees of freedom, which is impossible to do using purely deterministic methods.

The work of Bulatov et al. [29] and Marian et al. [34–38], among others, has developed a dedicated framework and pioneered the adaptation of the stochastic simulation algorithm (SSA) of Lanore and Gillespie [31–33] to cluster dynamics for irradiated materials. Notably, this approach has been successfully applied to the irradiation of W, accounting for the multispecies (He and H) for the first time. These applications demonstrated the usefulness of SSA-type methods in addressing the dimensionality issue of multispecies systems, thanks to sparse stochastic sampling from the underlying rate equation. However, SSA has been shown to be quite inefficient in systems exhibiting major stiffness issues, like typical RPV applications, as the time increment between two simulated reactive events becomes very small. Implementing hybrid deterministic/stochastic approaches only partially solves this issue [30]. The approach proposed by Vasav et al. [6] further exploits the concepts of partial propensities [7–9], of internal times and the use of binary heap to limit updates to reaction lists [10, 11]. Order-of-magnitude speed-ups are expected compared to the standard stochastic algorithm (SSA) while maintaining the capacity to mitigate the computational burden linked to the multidimensional clusters composition space.

To that end, CD equations must be reformulated as chemical master equations. Let us consider a homogenized system within a volume V with N clusters types (chemical species) S₁, …, S_N involved in M reactions R₁, …, R_M. The state of the system is modeled by the vector $X\left(t\right)=(X_1\left(t\right),...X_n\left(t\right))$, where X_α represents the number of clusters of type α ≡ ψ(n, p) and ψ is the function that gives the cluster index based on its chemical composition.

The reaction probability for each reaction R_j is specified in terms of a propensity function a_j = a_j(X(t), t), which equals the product of the number of possible combinations of reactive clusters involved in reaction R_j and a rate constant K_j, representing the emission or absorption rate.

Once a reaction R_j is carried out, the number of molecules for each species is updated according to the state change vector ν_j, that is, X(t)←X(t)+ν_j. These notations allow us to formulate the chemical master equation describing the evolution of the probability occupancy vector P over time:

$$\begin{array}{cc}\hfill \frac{\mathrm{d}}{\mathrm{d}t}P[X\left(t\right)=x]& =\sum _{j=1}^M{a}_j(x-{\nu }_j)P[X\left(t\right)-{\nu }_j=x]\hfill \\ \hfill & -a_j\left(x\right)P[X\left(t\right)=x].\hfill \end{array}$$(9)

For an emission or absorption reaction of a mobile cluster μ, S_α → S_α′ + S_μ, S_α + S_μ → S_α′, we obtain a_j = A_α, μX_α(t) and a_j = B_α, μX_α(t)X_μ(t)/V, respectively [32, 33]. This connection between cluster reaction rates highlights the dependence of propensities on the number of clusters of each reacting species in the associated reaction. The propensities associated with bimolecular reactions (controlling the agglomeration of two clusters for example), therefore depend on the product of the number of clusters of the two species reacting by a process controlled by the diffusion of a mobile cluster. Consequently, after a bimolecular reaction, the propensities associated with each of the two defect clusters are recalculated and associated reaction times removed from the priority list. These tasks represent a substantial part of the overall computational cost of standard SSA. We showed how to avoid these tasks by combining an internal time scale and a binary stack to each mobile cluster, inspired by the “next-reaction algorithm” by Gibson and Bruck [10] and Slepoy et al. [11]. Its algorithmic complexity is linear with respect to the number of new reactions to evaluate after each reaction and logarithmic with respect to the total number of potential reactions [10] if a binary stack is used to prioritize the times of all possible reactions. Depending on the cluster type and type of reactions, we can divide all reactions into separate classes:

• The class {0} contains all emission reactions from clusters but also the source terms due to irradiation (see [6] for details on the corresponding queue)

• A class is considered for each type of mobile cluster S_μ. The reactions in this class correspond to the potential absorptions of the mobile cluster S_μ by all the immobile clusters and by the mobile clusters S_μ′ such that μ′≥μ.

We then construct a priority queue P_μ containing all reactions of this class. This classification of reactions into separate classes and the use of internal times allow us to avoid re-calculating the propensities (rates) for all the reactions, just the ones that were involved in the previous reaction. For reaction between clusters α and μ, given the internal time τ_μ associated with P_μ, the potential internal time increments $\delta {\tilde{\tau }}_{\alpha ,\mu }$ are drawn from the exponential distribution of parameter${\tilde{a}}_{\alpha ,\mu }=B_{\alpha ,\mu }{X}_{\alpha }$, the partial propensity that was made independent of X_μ on purpose:

$$\begin{array}{c}\hfill \delta {\tilde{\tau }}_{\alpha ,\mu }\propto Exp({\tilde{a}}_{\alpha ,\mu }).\end{array}$$(10)

The next event candidate from the priority queue P_μ is the reaction (α, μ) such that:

$$\begin{array}{c}\hfill \delta {\tilde{\tau }}_{\mu }=\delta {\tilde{\tau }}_{\alpha ,\mu }.\end{array}$$(11)

At a given physical time t, the physical time increment until the next reaction corresponds to the minimum of internal times rescaled with X_μ(t), the current number of clusters in each class:

$$\begin{array}{c}\hfill \delta t=\min (\delta {\tilde{\tau }}_{\mu }/X_{\mu }(t)).\end{array}$$(12)

Where M is the set of mobile clusters classes.

Owing to its logarithmic complexity with respect to immobile clusters, and linear complexity with respect to the number of mobile species, this algorithm happens to be very efficient when the number of mobile species is small, which corresponds to the typical cases of cluster dynamics. One typical application of the method is the microstructure evolution in the FeHe system. The characteristics of the parameterization we use can be found in [30]. Notably, this parameterization does not allow the formation of interstitial loops containing helium atoms. The helium atoms can only be found in the matrix or in vacancy clusters, thus forming helium bubbles.

The first simulation was carried out with a volume of V = 5 × 10⁻⁸ cm³, which allows for the simulation of a physical time of 1353 seconds after 24 hours calculation on a single core. Figure 6 shows the distribution profiles of clusters as a function of their helium and vacancy composition (top panel) for the microstructure obtained at the end of the simulation. The color of the symbols indicates the number of helium-vacancy clusters. A logarithmic color scale is used to highlight the contrast in cluster densities according to their size. The distribution of interstitial dislocation loops is plotted as a function of the number of interstitial atoms in the loop in the bottom panel. Figures 7 and 8 also show the distribution profiles obtained for times of 2.9 × 10⁴ s and 2.8 × 10⁷ s, with simulation volumes of 2 × 10⁻¹⁰ cm³ and 2 × 10⁻¹² cm³,respectively.

Fig. 6. Distribution of He-vacancy and He-interstitial clusters at time t = 1353 s during the stochastic cluster dynamics simulation of the FeHe system with a volume of V = 5 × 10−8 cm3.

Fig. 7. Distribution of He-vacancy and He-interstitial clusters at time t=1353 s during the stochastic cluster dynamics simulation of the FeHe system with a volume of V=2 × 10−10cm3.

Fig. 8. Distribution of He-vacancy and He-interstitial clusters at time t = 2.859 × 107 s during the stochastic cluster dynamics simulation of the FeHe system with a volume of V = 2 × 10−12 cm3.

We observe that using smaller volumes allows for the simulation of longer times, which is easily explained by the fact that the smaller the system volume, the fewer the number of clusters. We also observe that with a small simulated volume, it is only possible to infer the regions of the phase space with high cluster densities, but post-processing is required to estimate the isodensity lines.

We now show that the number of clusters simulated by stochastic dynamics follows a linear scaling law with the simulated system volume beyond a certain threshold. To demonstrate this, we plotted on the left graph of Figure 6 the densities of large clusters as a function of physical time during simulations performed at constant volume. The density of large clusters is obtained by dividing the number of recorded clusters by the system volume. We observe that the density curves obtained with different volumes overlap well in Figure 9, highlighting the presence of a linear scaling law for the number of clusters as a function of the simulated volume beyond the threshold of 5 ×10⁻¹⁰ cm³. Below this threshold, significant side effects occur. This scaling law validates the stochastic cluster dynamics approach, whose goal is to reproduce deterministic cluster dynamics simulations that consider the limit of an infinite-volume system.

Fig. 9. Left: Master curve of the density of large clusters as a function of physical time during the simulation. Right: Estimated acceleration factor during the simulation. A cluster is considered large if it contains more than 250 vacancies and more than 250 helium atoms.

In the right panel of Figure 9, we show the computational acceleration enabled by the use of binary heaps compared to the direct algorithm, as a function of the simulated physical time. The acceleration factor is occasionally estimated during the simulation by directly summing the propensity functions associated with all possible reactions. We observe a very significant acceleration—nearly two orders of magnitude. This acceleration factor tends to decrease as the simulation progresses. This trend is well explained by the use of a memory-saving procedure that we implemented (“embedded clusters” see [6] for details).

For the application of FeCu alloys [6] up to 4 orders of magnitude increase in the number of reactions were achieved for similar CPU run-times compared to the original Gillespie-Lanore SSA algorithm. This will make it possible to treat multi-solute CD systems efficiently by considering two or more solute atoms types in mixed defect-solute clusters.

4. Kinetic activation relaxation technique with on-the fly prefactor calculations

In the previous sections, we have seen that algorithmic developments on kinetic Monte Carlo methods allow for increased efficiency in physical time simulation. Now, we turn to developments that allow for more advanced physical models. Diffusion rates of simple point defects, such as vacancies and interstitials, are generally assumed to depend on energy barriers only, as evidenced by the large number of kinetic Monte Carlo studies that rely on constant jump rate prefactors (also known as attempt frequencies):

$$\begin{array}{c}\hfill \mathrm{\Gamma }\left(T\right)={\mathrm{\Gamma }}_0{e}^{-{(E}_s-E_i)/kT},\end{array}$$(13)

where Γ is the jump rate, Γ₀ is the attempt frequency, T the temperature, E_i and E_s are the initial and saddle point energies, respectively, and k is the Boltzmann constant.

While this is mostly correct when energy barriers are well separated from each other in systems with little elastic deformation (i.e. typically in simple elemental crystals and ordered alloys), entropic variations between the local minima and associated activated states become increasingly important as energy barriers of diffusion-mediating mechanisms get closer or the elastic strain becomes significant, such as in disordered environments. In many systems, a linear correlation is observed between the prefactor and the activation barriers. This phenomenon corresponds to the long-recognized Meyer–Neldel (MN) compensation rule [39, 40] which, in glasses and amorphous systems, can be explained as a correlation between local geometry and vibrational models [41]. While the disorder in high entropy alloys is only chemical, the unexpected slowdown in defect diffusion observed for these systems, characterized by the presence of many different elements in roughly equal proportions, led us to revisit the role of prefactors.

We use the Activation-Relaxation Technique (ART nouveau) [42, 43] to sample the energy landscape and compute the transition rate for each specific event using the harmonic Transition State Theory (hTST) [44], an approximation that has been shown to capture the essential part of the entropic contributions for these mechanisms [41]. At variance with rigid lattice KMC, this relaxed-lattice kinetic Monte Carlo technique naturally takes into account elastic and lattice distortion effects, which are crucial in disordered systems such as high entropy alloys and disordered systems in general. Monte Carlo Metropolis type trial displacements centered on all first, second and third-neighbour atoms surrounding a vacancy are randomly generated to build an extensive catalogue of events containing the initial, saddle and final configurations. For each event, the prefactor is either computed within the harmonic approximation as described in Equation (3) or given a constant value of 10¹³ s⁻¹. Combining the prefactor with the energy barriers, we obtain the total rate associated with moving out of the current configuration. To facilitate the operations on the event-based catalogue and ensure that specific events are only counted once, we use, more specifically, the kinetic Activation-Relaxation Technique (k-ART) package [45, 46] and launch at least 25 event searches per unique topology for each atoms surrounding the vacancy. This extensive search ensures a rich sampling of the energy landscape surrounding the vacancy. The harmonic Transition State Theory (hTST) [44] defines the transition rate prefactor as

$$\begin{array}{c}\hfill {\nu }_{\mathit{hTST}}=\frac{\prod _{i=1}^N{\nu }_i^m}{\prod _{i=1}^{N-1}{\nu }_i^s},\end{array}$$(14)

where ${\nu }_i^s$ and ${\nu }_i^m$ are real vibrational frequencies at the saddle point and the minimum, respectively, and are obtained by computing the dynamical matrix

$$\begin{array}{c}\hfill D_{i\alpha j\beta }=\frac{1}{\sqrt{m_i{m}_j}}\frac{{\partial }^{2}V}{\partial x_{i,\alpha }\partial x_{j,\beta }},\end{array}$$(15)

where i, j run over all atoms, α, β are the cartesian coordinates (x, y, z), V is the interaction potential, and m_i the atomic mass. The matrix can be obtained through a centred-finite difference formula with small finite displacements of 0.01 Å, for the selected potential. Frequencies correspond to the square root of the dynamical matrix eigenvalues. Combining the Activation-Relaxation Technique nouveau (ART nouveau) and the harmonic approximation for computing diffusion prefactors (referred as hTST), we find that vacancy diffusion prefactors in a 55Fe-28Ni-17Cr concentrated solution alloy modelled with EAM empirical potentials can vary by up to six orders of magnitude, at almost constant energy barrier as illustrated on Figure 10 [47].

Fig. 10. Transition rates calculated using the harmonic Transition State Theory (hTST) prefactor are plotted against rates computed with a constant prefactor (bottom axis) and corresponding energy barriers (top axis). The various dashed lines indicate where the data points would lie if the hTST rates matched those computed with a constant prefactor from 109 to 1015 s−1[47].

This amplitude of variation, mostly associated with changes in local pressure, suggests that prefactor could play a much more important role than previously thought in the defect kinetics of high entropy alloys [48] and of disordered systems such as the concentrated alloys for austenitic stainless steels.

We also calculated the vacancy diffusion coefficient by monitoring the mean square displacement (MSD) and how the system’s total energy dissipates through this process [48]. This analysis revealed a consistent trend: vacancy diffusion is not random but rather exhibits a marked preference for migration into specific local atomic environments, which can be understood by considering local strain fields and interatomic interactions.

In particular, the vacancy shows a tendency to diffuse toward Fe-depleted regions, which is linked to the chemical strain induced by Ni–Cr bonds, the longest among the alloy’s binary interactions. These longer bonds introduce local compressive strain, lowering the energy barriers for nearby atomic rearrangements and making these regions more accessible to vacancy migration. This behavior is especially prominent in simulations using a constant attempt frequency, where diffusion pathways are selected solely based on energy considerations.

However, the picture changes substantially when we introduce vibrational entropy effects via harmonic Transition State Theory (hTST). In these simulations, the prefactors for each transition are calculated explicitly based on the local vibrational spectra at the initial and saddle-point configurations. This introduces a broad distribution of prefactors, which effectively acts as an entropic weighting factor on the transition rates. Consequently, while some high-barrier events become more kinetically viable due to large prefactors, others–despite being energetically favorable–are suppressed due to lower vibrational contributions. The net result is a slower system evolution, as reflected in the delayed atomic and vacancy displacement in hTST simulations, and a reduced rate of energy dissipation, as seen in the flatter energy relaxation curves seen in Figure 11.

Fig. 11. Comparison of vacancy diffusion and system energetics from 50 KMC simulations (∼2000 steps each) using the harmonic Transition State Theory (hTST) prefactor (orange) and 50 KMC simulations (∼2000 steps each) using a constant prefactor (blue). (top) Mean square displacement of a single vacancy. (bottom) Evolution of the system’s total energy [48].

This interplay between energy barriers and vibrational prefactors is further illuminated by analyzing the statistical correlation between the two. When considering the full catalog of potential events generated during the simulations, we observe that Fe exhibits a negative correlation, meaning higher energy barriers are not compensated by correspondingly high prefactors. Ni and Cr show little to no correlation. These results, illustrated in Figure 12 (top), deviate from the expectations of the Meyer–Neldel compensation rule, which posits a positive correlation between activation energy and prefactor due to entropic contributions compensating for higher barriers.

Fig. 12. Logarithmic values of hTST prefactors plotted against their corresponding energy barriers for Fe, Ni, and Cr. (top) All the possible events from the simulations are included, with compensation effects determined via linear fits to events with barriers above 0.3 eV, yielding slopes of −0.587(0.006) eV−1 (Fe), −0.070(0.002) eV−1 (Ni), and 0.017(0.004) eV−1 (Cr). (bottom) A curated subset of selected events above 0.3 eV shows compensation slopes of 0.47(0.02) eV−1, 0.49(0.02) eV−1, and 1.29(0.04) eV−1 for Fe, Ni, and Cr, respectively. Values in parentheses represent the standard deviation of the linear fit [48].

However, a critical distinction emerges when we restrict the analysis to kinetically selected events, i.e., those that are actually executed during the kMC evolution. As shown in Figure 12 (bottom), these events do conform to the Meyer–Neldel rule, exhibiting a clear positive correlation between prefactors and barriers across all species. This indicates that the system’s dynamics naturally favor pathways that balance energy and entropy, selecting higher-barrier transitions only when they are entropically supported by large attempt frequencies. Notably, this compensation effect is strongest for Cr, which already has lower average barriers, while Fe and Ni and though less mobile in the full event set, gain prominence in the kinetics due to favorable prefactor-barrier combinations.

This shift underscores an important conceptual point: in some systems, diffusion may be governed not just by the availability of low-energy pathways, but by a complex competition between energy and entropy. Methods such as kART-hTST are needed to capture this balance, enabling the identification of truly relevant diffusion mechanisms rather than those that are merely energetically accessible.

More specifically, these findings emphasize that the explicit treatment of vibrational entropy is essential for accurately modeling diffusion in chemically disordered systems. Simplified, constant-prefactor models can misrepresent not only diffusion rates but also the species-specific roles in defect transport, potentially leading to flawed predictions of microstructural evolution.

In summary, the correlation between energy barriers and prefactors reveals vibrational entropy as a dynamic moderator of atomic diffusion, not a passive thermodynamic quantity. By integrating energetic and entropic effects, tools like kART-hTST bridge atomistic detail and mesoscale behavior, providing powerful capabilities for predicting the long-term evolution of complex alloys, including high-entropy and radiation-tolerant materials.

5. Conclusion

Atomic Monte Carlo methods are central in irradiated microstructures modelling. In the case of rigid lattice kMC, the ability of the algorithm to reach long evolution times is often limited by kinetic trapping of the simulated trajectories within low energy basins presenting small intra-basin barriers. This can be alleviated by implementing non-local jumps relying on the theory of absorbing Markov chains. To reach even larger time scales, many models rely on huge systems of cluster dynamics equations which are usually solved deterministically. Monte Carlo schemes can solve them using dedicated algorithms relying on binary heaps to make them more efficient. To study the origin phenomena like slow diffusion processes in HEA with kMC techniques, one must first address other types of limitations of kMC such as, the rigid lattice approximation, the absence of elastic fields and of vibrational entropy contributions. An improved version of the activation-relaxation technique can be used for such a purpose and reveals substantial effects of the spread in jump rate prefactors.

These kMC methods and developments allow to accurately determine the properties of complex defects and their evolution, which provide inputs for OKMC and CD for coarse graining the microstructure evolution under thermal ageing or irradiation.

Funding

G.A., R.V., T.J., C.D., M.A. thanks I3P (Institut Tripartite with CEA, EDF and Framatome) for co-funding the study. N.M. acknowledges partial support through a Discovery grant from the Natural Science and Engineering Research Council of Canada.

Conflicts of interest

The authors have nothing to disclose.

Data availability statement

Data associated with this article cannot be disclosed due to legal reasons.

Author contribution statement

Conceptualization, M.A. and N.M.; Methodology, M.A., T.J., G.A. and N.M.; Software, R.V., T.J., G.A., M.N., M.A. and M.M.R.; Validation, R.V., T.J., G.A., C.D., M.A., M.N. and M.M.R.; Formal Analysis, M.A., R.V., T.J., G.A., C.D., M.N. and M.M.R.; Investigation, R.V., T.J., G.A., C.D., M.A., M.N. and M.M.R.; Resources, M.A., M.N., T.J., G.A. and C.D.; Data Curation, R.V., T.J., G.A., C.D., M.A., M.N. and M.M.R.; Writing – Original Draft Preparation, G.A., R.V., T.J., C.D., M.A., M.N. and M.M.R.; Writing – Review & Editing, G.A., R.V., T.J., M.A., C.D., M.N. and M.M.R.; Visualization, G.A., R.V., T.J., C.D., M.A., M.N. and M.M.R.; Supervision, M.A., T.J. and M.N.; Project Administration, G.A., T.J., C.D., M.A. and M.N.; Funding Acquisition, M.A., M.N., T.J., G.A. and C.D.

References

All Figures

	Fig. 1. (a) A cluster of B solute atoms within a bulk of A atoms is sketched. In insets (b,c,d) a typical “flicker” is illustrated. It results in having the vacancy (V) returning at its initial position many times. This is due to B–V exchanges being far more likely than A-V ones (more complex flickers are possible (e)) as reflected by their respective transition rates Γ.
In the text

	Fig. 2. The trapping sites within the solute cluster (shown in blue) determine the N transient states, while the perimeter states determine the absorbing states. The former are associated with the A blocks of the matrix and the latter with the null blocks. At this stage, all perimeter sites are treated as a single absorbing state with an index of N + 1.
In the text

	Fig. 3. Integrated speed-up plotted as a function of the physical time simulated by kPS at three different temperatures T0 = 273 K, T1 = 373 K, and T2 = 473 K [12].
In the text

	Fig. 4. Diffusion coefficients of Mn solute atoms as a function of cluster size for various temperatures T. For clarity, scaling is done relative to monomer diffusivity at T0 and its high-temperature activation energy Ea = 0.646 eV. Reference temperature is T0 = 600 K and reference diffusivity is D(1,T0)=6.716 × 10−15 m2 s−1 at vacancy concentration of 0.05% [13].
In the text

	Fig. 5. Comparison of algorithm efficiencies as a function of cluster size. Efficiency is measured by the ratio of simulated physical time to that of F-KMC (kPS), based on a fixed wall-clock time. Shaded regions around the curves represent the 68% confidence interval.
In the text

	Fig. 6. Distribution of He-vacancy and He-interstitial clusters at time t = 1353 s during the stochastic cluster dynamics simulation of the FeHe system with a volume of V = 5 × 10−8 cm3.
In the text

	Fig. 7. Distribution of He-vacancy and He-interstitial clusters at time t=1353 s during the stochastic cluster dynamics simulation of the FeHe system with a volume of V=2 × 10−10cm3.
In the text

	Fig. 8. Distribution of He-vacancy and He-interstitial clusters at time t = 2.859 × 107 s during the stochastic cluster dynamics simulation of the FeHe system with a volume of V = 2 × 10−12 cm3.
In the text

	Fig. 9. Left: Master curve of the density of large clusters as a function of physical time during the simulation. Right: Estimated acceleration factor during the simulation. A cluster is considered large if it contains more than 250 vacancies and more than 250 helium atoms.
In the text

	Fig. 10. Transition rates calculated using the harmonic Transition State Theory (hTST) prefactor are plotted against rates computed with a constant prefactor (bottom axis) and corresponding energy barriers (top axis). The various dashed lines indicate where the data points would lie if the hTST rates matched those computed with a constant prefactor from 109 to 1015 s−1[47].
In the text

	Fig. 11. Comparison of vacancy diffusion and system energetics from 50 KMC simulations (∼2000 steps each) using the harmonic Transition State Theory (hTST) prefactor (orange) and 50 KMC simulations (∼2000 steps each) using a constant prefactor (blue). (top) Mean square displacement of a single vacancy. (bottom) Evolution of the system’s total energy [48].
In the text

Fig. 12. Logarithmic values of hTST prefactors plotted against their corresponding energy barriers for Fe, Ni, and Cr. (top) All the possible events from the simulations are included, with compensation effects determined via linear fits to events with barriers above 0.3 eV, yielding slopes of −0.587(0.006) eV−1 (Fe), −0.070(0.002) eV−1 (Ni), and 0.017(0.004) eV−1 (Cr). (bottom) A curated subset of selected events above 0.3 eV shows compensation slopes of 0.47(0.02) eV−1, 0.49(0.02) eV−1, and 1.29(0.04) eV−1 for Fe, Ni, and Cr, respectively. Values in parentheses represent the standard deviation of the linear fit [48].

In the text