Atomic scale modeling of the thermodynamic and kinetic properties of dilute alloys driven by forced atomic relocations

TL;DR

This paper develops an extended theoretical framework to model the thermodynamic and kinetic properties of dilute alloys under external forces, revealing non-equilibrium effects and deviations from classical relations, with applications to radiation damage.

Contribution

It introduces a modified Self-Consistent Mean Field theory incorporating forced atomic relocations, enabling accurate modeling of driven alloy systems beyond equilibrium assumptions.

Findings

FAR significantly affects vacancy-solute flux coupling.

Onsager reciprocal relations are violated in driven systems.

The extended model agrees with Monte Carlo simulations.

Abstract

Sustained external forces acting on a material provide additional mechanisms to evolve the state of the system, and these mechanisms do not necessarily obey the microscopic detailed balance. Therefore, standard methods to compute the thermodynamic and kinetic properties do not apply in such driven systems. The competition between these mechanisms and the standard thermally activated jumps leads to non-equilibrium steady states. We extend the Self-Consistent Mean Field theory to take into account forced atomic relocations (FARs) as a model of these additional kinetic mechanisms. The theory is applied to the atomic-scale modelling of radiation damage. Using a first-shell approximation of the theory, we highlight the violation of Onsager reciprocal relations in driven systems. An implementation of the extended theory into the KineCluE code yields calculations of the effective Onsager coefficients in agreement with Monte Carlo simulations. A systematic parametric study is performed to emphasize the effect of FAR distances and the solute-defect interaction on the diffusion properties. The effect of FAR on the vacancy-solute flux coupling and the solute tracer diffusivity is non-negligible when: (i) the solute-vacancy thermodynamic attraction is large, (ii) the magnitude of the thermal jump frequencies is lower or comparable to the frequencies of FAR, and (iii) the range of interactions between vacancies and solute atoms is close to FAR distances.

Tables

Table 1. Table 1: List of the parameters needed to estimate the vacancy concentration and their values set in the paper.

Parameter	Value
Lattice parameter $a_0$	0.35 nm
Vacancy formation enthalpy $H_V^{\text{f}}$	1.65 eV
Vacancy formation entropy $S_V^{\text{f}}$	1.82 $k_B$
Number of relocations per dpa $n_{\text{rel}}$	100
Effective sink strength $k^2$	${10}^{15}{\text{m}}^{-2}$

Equations

$${\omega }_{\text{n}\to \tilde{\text{n}}}^{\alpha V}=\nu \text{exp}\left(-\frac{E_{\text{n}\to \tilde{\text{n}}}^{\text{mig}}}{k_{B}T}\right),$$

$$P_{\text{n}}{\omega }_{\text{n}\to \tilde{\text{n}}}^{\alpha V}=P_{\tilde{\text{n}}}{\omega }_{\tilde{\text{n}}\to \text{n}}^{\alpha V},$$

$$p(r)=\frac{1}{r_{\text{m}}}\text{exp}\left(-\frac{r}{r_{\text{m}}}\right),$$

$$\mathcal{P}(i)={\int }_{r_i}^{r_{i+1}}p(r)\text{d}r,$$

$${\mathcal{P}}_L(i)=\frac{\mathcal{P}(i)}{{\sum }_{s=0}^L\mathcal{P}(s)}.$$

$${\Gamma }^{\text{aa}}=n_{\text{rel}}\varphi .$$

$${\Gamma }^{\text{ad}}=\gamma {\Gamma }^{\text{aa}},$$

$$C_V^{\text{ness}}=C_V^{\text{eq}}-\frac{k^2\Omega }{8\pi r_c}+\sqrt{{\left(\frac{k^2\Omega }{8\pi r_c}\right)}^2+\frac{\varphi \Omega }{4\pi r_c{D}_V}},$$

$$C_V^{\text{eq}}=\text{exp}\left(-\frac{H_V^{\text{f}}-T{S}_V^{\text{f}}}{k_{B}T}\right).$$

$$W_{\text{n}\to \tilde{\text{n}}}=W_{\text{n}\to \tilde{\text{n}}}^{AV}+W_{\text{n}\to \tilde{\text{n}}}^{BV}+W_{\text{n}\to \tilde{\text{n}}}^{AB},$$

$$\frac{\text{d}}{\text{d}t}\mathbf{P}=\mathbf{W}\mathbf{P},$$

$$\forall \text{n},\sum _{\tilde{\text{n}}}{W}_{\tilde{\text{n}}\to \text{n}}{P}_{\tilde{\text{n}}}-W_{\text{n}\to \tilde{\text{n}}}{P}_{\text{n}}=0.$$

$${\overrightarrow{J}}_{\alpha }=-\sum _{\beta }{\lambda }_{\alpha \beta }\frac{\overrightarrow{\nabla }{\mu }_{\beta }}{k_{B}T},$$

$$P_{\text{n}}(t)=P_{\text{n}}^{\text{ness}}\times \delta P_{\text{n}}(t).$$

$$\frac{\text{d}{P}_{\text{n}}(t)}{\text{d}t}=\sum _{\tilde{n}}\left[W_{\tilde{\text{n}}\to \text{n}}{P}_{\tilde{\text{n}}}^{\text{ness}}\delta P_{\tilde{\text{n}}}(t)-W_{\text{n}\to \tilde{\text{n}}}{P}_{\text{n}}^{\text{ness}}\delta P_{\text{n}}(t)\right].$$

$$\frac{\text{d}{P}_{\text{n}}(t)}{\text{d}t}=\sum _{\tilde{n}}{W}_{\tilde{\text{n}}\to \text{n}}{P}_{\tilde{\text{n}}}^{\text{ness}}\left[\delta P_{\tilde{\text{n}}}(t)-\delta P_{\text{n}}(t)\right].$$

$$W_i^{AV}$$

$$W_i^{AB}$$

$$W_2^{BV}$$

$$\frac{P_1^{\text{ness}}}{P_0^{\text{ness}}}$$

$${\lambda }_{VB}=$$

$${\lambda }_{BV}=$$

$${\lambda }_{VV}=$$

$${\lambda }_{BB}=$$

$$\left\{\begin{array}[]{lr}C_{BV}^{p}=C_{B}^{0}C_{V}^{0}z^{\text{ness}}\\ C_{V}^{m}=C_{V}-C_{BV}^{p}\\ C_{B}^{m}=C_{B}-C_{BV}^{p},\end{array}\right.$$

$$\left\{\begin{array}[]{lr}C_{B}=C_{B}^{0}+C_{B}^{0}C_{V}^{0}\left(z^{\text{ness}}-z_{o}\right)\\ C_{V}=C_{V}^{0}+C_{B}^{0}C_{V}^{0}\left(z^{\text{ness}}-z_{o}\right),\end{array}\right.$$

$${\lambda }_{BV}={\lambda }_{VB}=-\frac{a_0^2{C}_{BV}^p}{4}\frac{W_2^{BV}\left(13W_3^{AV}-2W_1^{AV}\right)}{7W_3^{AV}+2W_1^{AV}+2W_2^{BV}}.$$

$$\Delta \lambda =\frac{3a_0^2{C}_{BV}^p}{4}\frac{(1-{\omega }_3/{\omega }_4)[{\omega }_a-(1-\gamma )\Gamma ]{\omega }_4\Gamma }{({\omega }_b+11\gamma \Gamma +9\Gamma )\left[{\omega }_4+(1+\gamma )\Gamma \right]},$$

$${\omega }_a$$

$${\omega }_b$$

$$D_B^{\ast }=\frac{{\lambda }_{BB}}{C_B}.$$

$$D_{B,\text{add}}^{\ast }=D_{B,\text{th}}^{\ast }{C}_V^{\text{ness}}/C_V^{\text{eq}}+D_{B,\text{far}}^{\ast },$$

$$\Delta D_B=\frac{D_{B,\text{add}}^{\ast }-D_B^{\ast }}{D_B^{\ast }}.$$

Full text

††preprint: APS/123-QED

I Introduction

Excess point defects (PDs) are massively generated in materials driven away from equilibrium by external forces such as mechanical solicitation and irradiation Martin and Bellon (1996); Was (2007). This phenomenon induces fluxes of PDs toward the microstructural features (e.g., grain boundaries, dislocation lines, etc.) acting as PD sinks. The solute-PD interaction in alloys leads to the coupling between the PD and solute atom fluxes to the sinks. This flux coupling is the main kinetic process controlling the redistribution of solute atoms in alloys driven by an excess of PDs Anthony and Hanneman (1968); Anthony (1969); Okamoto and Wiedersich (1974); Barbu and Ardell (1975); Okamoto (1979); Kato et al. (1992); Bruemmer et al. (1999); Nastar and Soisson (2012); Ardell and Bellon (2016).

Additionally, the external forces may introduce new mechanisms that affect the state of the system at the atomic scale. Examples of these mechanisms include collective motion of atoms or mixing. For instance, the latter mechanisms occur under irradiation producing displacement cascades Haff and Switkowski (1977); Motta and Lemaignan (1992); Averback and De La Rubia (1998) or during severe plastic deformation such as shearing Vo et al. (2013); Ashkenazy et al. (2017), torsion Pouryazdan et al. (2012); Beach et al. (2017) and ball milling P. Pochet, E. Tominez, L. Chaffron and Martin (1995); Klassen et al. (1997); Suryanarayana (2004). Unlike the thermally activated mechanisms leading the system toward equilibrium, these additional mechanisms are mostly athermal. They do not obey the microscopic detailed balance and lead to the disordering of the atomic configurations Martin and Bellon (1996). Note that both the thermally activated and externally forced mechanisms contribute to the mass transport in driven alloys, such that the system is driven to a non-equilibrium steady state (NESS) Martin and Bellon (1996). This situation prevents the use of standard methods to compute the thermodynamic and kinetic properties. In order to understand and predict the effect of external driving forces, it is crucial to be able to model the interplay between the two above-mentioned mechanisms. However, there is no model of flux coupling simultaneously accounting for these two competitive mechanisms at the atomic scale.

Radiation damage is a good school case for the modeling of driven alloy thermodynamics and kinetics because there exists microscopic models that effectively reproduce the atomic mixing caused by the external force Martin (1984); Averback (1986). Under irradiation, atoms are regularly hit by incident particles. The collision transfers kinetic energy from the incident particle to the primary knock-on atom (PKA). If this energy is below the displacement threshold energy (DTE), the PKA will only vibrate around its position unless a PD is located nearby, in which case the atom may exchange its position with the PD Roussel and Bellon (2002). For a recoil energy well above DTE (e.g., typically 1 keV in metals), the PKA will move away from its original site, thereby creating Frenkel pairs and transferring kinetic energy to neighbouring atoms, which will themselves move away from their positions, so on so forth. Locally, this displacement cascade process produces a large excess number of PDs as well as atomic mixing since most of the atoms are displaced Littmark and Hofer (1980). The climax of this process is called the heat spike, where the material is locally liquid-like and atoms are able to rapidly diffuse in this region Vineyard (1976). The excess energy eventually dissipates leading to a quench-like process where the crystalline structure is recovered and only a small excess number of PDs remain Brinkman (1954); Benedek et al. (1987); Nordlund and Averback (1999). The whole phenomenon can be effectively modeled by (1) creating PDs and PD clusters, and (2) shuffling atomic positions. The latter process is often considered random English and Jenkins (1981); Pramanik and Seidman (1986) but there are some evidences that it is in fact affected by thermodynamics Workman et al. (1987). Besides, recent investigations on concentrated alloys Terentyev et al. (2006); Aidhy et al. (2015); Zhang et al. (2017) and high entropy alloys Do and Lee (2018) have shown that during the quenching stage, the spatial distribution of the PDs compared with the solute atoms is partially driven by the thermodynamic short range order (i.e. binding interaction). For more details concerning displacement cascades, please refer to Ref. Nordlund et al. (2018).

The mixing of atomic positions in the displacement cascade was previously modeled effectively by forced atomic relocations (FARs) which consist in forced exchanges of positions between an atom and its nearest neighbour species such as another atom or a PD Martin (1984); Enrique and Bellon (2000); Lear et al. (2017); Soisson et al. (2018). However, details in the cascade such as the spatial correlation between solute atoms and PDs due to the thermodynamic interactions were neglected. Moreover, the thermal mechanism and FAR were considered separately at the macroscopic scale. For instance, the tracer diffusion coefficient of the solute atom was written as the sum of two diffusion coefficients respectively related to the two mechanisms Was (2007). In this case, the interplay between these two mechanisms was neglected. An improved version of the model was proposed in Ref. Roussel and Bellon (2002), where the five-frequency model Lidiard (1955, 1960) in face-centered cubic (fcc) systems was generalized to account for both types of kinetic mechanisms simultaneously. However, solute diffusivity is the only accessible quantity from this model, so it does not provide any information on solute-PD flux coupling. Moreover, long range FAR cannot be accounted for since the five-frequency model is limited to first nearest neighbour (1NN) interactions.

Thermodynamic and kinetic properties such as flux coupling coefficients and tracer diffusion coefficients of an alloy can be calculated from the Onsager matrix of the transport coefficients. Whenever the diffusion mechanism satisfies the microscopic detailed balance, Onsager has demonstrated that this matrix is symmetric Onsager (1931a, b). We may calculate it either from the equilibrium atomic displacement fluctuations using the Allnatt formulae Allnatt (1965, 1982) or from the flow of matter resulting from an applied external force. However, for driven alloys including athermal mechanisms not obeying the microscopic detailed balance, we cannot compute the transport coefficients by means of a Monte Carlo numerical approach based on the Allnatt formulae. Yet, recent statistical theories have shown that it is possible to derive an effective Onsager matrix from the fluctuation theorem Evans et al. (1993); Gallavotti (1996), though the resulting matrix is non symmetric. These theories go beyond the linear response theory. They provide a methodology for the investigation of far from equilibrium kinetics. Such an approach has been applied to the study of a molecular motor driven by forced chemical reactions Lau et al. (2007); Lacoste et al. (2008). However, it is not directly applicable to properly model systems with FAR because there are no notions of alloying effects and kinetic correlations in this model. In the context of research on diffusion in alloys, one knows how to deal with the complexity of calculating a sequence of PD jumps when the frequency of each jump depends on the local environment of the defect as long as the diffusion mechanism satisfies the microscopic detailed balance Allnatt and Lidiard (1993); Nastar (2005); Garnier et al. (2013, 2014); Messina et al. (2014, 2016); Schuler and Nastar (2016); Jain et al. (2019).

To obtain the transport coefficients at NESS, we start from the standard self-consistent mean field (SCMF) theory Nastar et al. (2000); Nastar (2005), which was applied to the calculation of transport coefficients only in systems near equilibrium with jump mechanisms obeying the microscopic detailed balance. In this paper, we generalize the SCMF theory by including the athermal mechanism of FAR. The generalized SCMF is implemented in the KineCluE code Schuler et al. (2020) in order to perform automated calculations of transport coefficients. Different models of FAR are tested and the impact of each model parameter is systematically investigated. This study allows us to understand the conditions in which FAR significantly affect the material thermodynamic and kinetic properties.

The structure of the paper is as follow. In Sec. II we introduce the model of irradiation damage and thermal diffusion. Then, we introduce the mean field kinetic model to estimate the PD concentration under irradiation. Sec. III is devoted to the calculation of the transport coefficients within the SCMF framework. In Sec. IV, we focus on the comparison between the results given by different FAR models in various representative model alloys. The impact of different model parameters are studied in Sec. V. A discussion of the results, including the limitations and possible improvements, is presented in Sec. VI.

II Modeling of diffusion mechanisms under irradiation

II.1 Thermally activated jump frequencies

We use the transition state theory to model thermally activated diffusion Vineyard (1957). We introduce the thermal jump frequency $\omega_{\textbf{n}\rightarrow\widetilde{\textbf{n}}}^{\alpha V}$ associated with the thermally activated exchange of atom $\alpha$ and vacancy V which brings the system from configuration n to $\widetilde{\textbf{n}}$:

[TABLE]

where $\nu$ is the attempt frequency, $k_{B}$ is the Boltzmann constant, $T$ is the temperature and $E^{\text{mig}}_{\textbf{n}\rightarrow\widetilde{\textbf{n}}}$ is the migration barrier from configuration n to $\widetilde{\textbf{n}}$, which can be computed by means of ab initio calculations Tucker et al. (2010); Messina et al. (2014). This mechanism is mediated by PDs and the jump rate depends on the temperature as well as the initial and saddle-point configurations. Note that it satisfies the principle of the microscopic detailed balance:

[TABLE]

where $P_{\textbf{n}}$ is the probability of configuration n.

II.2 Irradiation damage

We follow the ideas of previous studies to model the radiation damage by FAR. The latter includes two mechanisms: (1) FAR between two randomly chosen atoms (FAR-a) which consists in exchanging the positions of two atoms on lattice sites, and (2) FAR between a randomly chosen atom and a PD (FAR-d) which consists in exchanging the positions between an atom and a vacancy (V) or a self-interstitial atom (SIA). We need to account for the removing and creation of PDs within a cascade. Here we consider that the PDs “disappear” during the heat spike, because the material is locally a liquid-like phase where there is no notion of PD. Later, during the process of quenching, only a small fraction of the PDs “reappears” somewhere in the cascade area, as if they had effectively jumped to another crystalline site during the cascade. The effective result of this process is modeled by FAR-d.

For recoil energy well above DTE producing displacement cascade, the overall effect of the mixing is modeled by FAR characterized by a relocation distance $r$. FAR occurs at a given frequency proportional to the radiation flux. First for FAR-a, we assume that the probability density function $p(r)$ of the relocation distance follows an exponential decay Enrique and Bellon (2000); R. A. Enrique, K. Nordlund, R. S. Averback and Bellon (2003); Demange et al. (2017):

[TABLE]

where $r_{\text{m}}$ is the mean relocation distance which is related to the size of the displacement cascade. Note that the latter depends on the material and the recoil energy of PKA. For example, the sizes of cascades generated in metals by fast neutrons or by heavy ions typically range between 10 and 100 Å English and Jenkins (1981); Phythian et al. (1995). At the atomic scale, the relocation distance $r$ is discrete and is equal to one of the $i$-th NN distances. We define the probability mass function $\mathcal{P}(i)$ so that the distribution $p(r)$ in the interval $[r_{i},r_{i+1}]$ is averaged to the $i$-NN point:

[TABLE]

where $r_{i}$ corresponds to the $i$-NN distance. In practice, we consider only a finite set of nearest neighbours, meaning that there is a cut-off relocation distance $L$-NN beyond which the probability is set to 0. In this case, we define the normalized probability mass function $\mathcal{P}_{L}(i)$ as:

[TABLE]

We introduce as well a simplified model associated with a single relocation distance $r_{m}$ because it gives access to an analytical solution.

Here we ignore FAR-d of SIA and this assumption is justified in Sec. II.3. Therefore, we consider only FAR-d of vacancy. We propose two categories of FAR-d models: either the same relocation model employed for FAR-a, or a model favoring the relocation sites close to the solute atoms in case of attractive binding energies between vacancy and solute atoms. The latter model makes sense because in the quench-like process at the end of the displacement cascade, the remaining PDs form preferentially where their formation energy is the lowest, that is in the vicinity of solute atoms.

In order to represent both categories of models, we introduce three models. Models 1 and 2 for the first category, and Model 3 for the second category including a thermodynamic effect on FAR-d. Model 1 includes a single relocation distance for both solute and vacancy, while Model 2 includes an exponential law for the relocation distance (Eq. (5)) for both species. Model 3 is similar to Model 2, the only difference is that when the relocated vacancy is located at a distance lower than a threshold value $R_{c}$ from the solute atom B, the vacancy is systematically exchanged with an atom randomly chosen among the 1NN atoms of B (chemically biased FAR-d).

For recoil energy below DTE, the effective result of the sub-threshold collision is modeled only by FAR-d. The model of FAR-d is the same in Model 1 while FAR-a is not performed. The relocation distance is set to 1NN distance $r_{1}$.

II.3 FAR frequencies

The FAR-d frequency is denoted $\Gamma^{\text{ad}}$, and the FAR-a frequency is denoted $\Gamma^{\text{aa}}$.

When the recoil energy is above DTE, the relocation frequency $\Gamma^{aa}$ can be deduced from the radiation dose rate $\phi$ based on the ion-beam mixing framework Haff and Switkowski (1977); Averback (1986). In our model, FAR-a reproduces the mixing of atoms in the displacement cascade, which is related to the number of PDs produced by the PKA. After the quenching phase, there is only a small fraction of surviving defects, which defines the unit of displacement per atom (dpa). Therefore, there is a factor $n_{\text{rel}}$ relating $\Gamma^{\text{aa}}$ and the radiation dose rate in unit of dpa/s (see Eq. (6)).

[TABLE]

From the literature, we set $n_{\text{rel}}=100$ Rivière and Dinhut (1983); Averback (1986); Müller et al. (1988). The latter number varies with the alloy thermodynamics due to the thermal effect on the atomic mixing rate in the cascade Workman et al. (1987). The frequencies of FAR-a and FAR-d depend on the number of cascades formed per unit of time as stated in Section II.2. Therefore, $\Gamma^{\text{aa}}$ and $\Gamma^{\text{ad}}$ are both proportional to the dose rate. Hence, they are proportionally related by:

[TABLE]

with $\gamma$ the proportionality constant. Note that $\gamma$ is set to 1 if not specified, i.e. $\Gamma^{ad}=\Gamma^{aa}=\Gamma$. Sensitivity studies concerning the value of $\gamma$ are shown in Sec. V.4.

When the recoil energy is below DTE, the sub-threshold irradiation does not induce FAR-a because no displacement cascade is produced, thus $\Gamma^{aa}=0$. In this case, the calculation of $\Gamma^{\text{ad}}$ is not related to $\Gamma^{aa}$ and is directly deduced from the recoil energy.

Note that the maximum dose rate under realistic irradiation condition is around 1 dpa/s Was (2007), thereby leading to a relocation frequency of about 100 s*-1*. The latter is still very small compared to the thermal jump frequency of SIA, even at low temperature. For instance, the SIA thermal jump frequency in pure nickel at 300 K is around $10^{10}$ s*-1* according to the atomic diffusion data given in Ref. Tucker et al. (2010). Therefore, we do not expect an important impact of FAR-d on the SIA-mediated diffusion properties. Therefore, we consider only FAR-d with vacancies, as stated in Sec. II.3. Yet, we emphasize that the extension of our framework to account for FAR-d of SIAs is straightforward.

II.4 Point defect concentration

The global concentration of PD varies under irradiation, mainly due to the production of Frenkel pairs, the mutual recombination between SIA and V, the elimination of PD at PD sinks such as grain boundaries and dislocations. The vacancy concentration at NESS $C_{V}^{\text{ness}}$ is estimated from a rate theory model Russell (1984); Schuler et al. (2017a):

[TABLE]

where $C_{V}^{\text{eq}}$ is the thermal vacancy concentration at equilibrium, $r_{c}$ is the SIA-V recombination radius usually assumed to be of the order of the lattice parameter $a_{0}$. $\Omega$ is the atomic volume. $\phi$ is the radiation dose rate. $k^{2}$ is the sink strength assumed to be constant with the radiation dose rate with typical values ranging from $10^{12}$ to $10^{19}$ m*-2* Soisson and Jourdan (2016) and $D_{V}$ is the vacancy diffusion coefficient. Note that the equilibrium concentration $C_{V}^{\text{eq}}$ is obtained from the vacancy formation enthalpy $H_{V}^{\text{f}}$ and entropy $S_{V}^{\text{f}}$ by:

[TABLE]

We may then replace the flux $\phi$ by its expression in terms of $\Gamma$ from Eq. (6) into Eq. (8), leading to a direct relationship between the vacancy concentration at NESS and $\Gamma$.

III Diffusion theory

For total exchange frequencies, we use the notation:

[TABLE]

where $W_{\textbf{n}\rightarrow\widetilde{\textbf{n}}}^{AV}=\omega_{\textbf{n}\rightarrow\widetilde{\textbf{n}}}^{AV}+\Gamma_{\textbf{n}\rightarrow\widetilde{\textbf{n}}}^{AV}$, $W_{\textbf{n}\rightarrow\widetilde{\textbf{n}}}^{BV}=\omega_{\textbf{n}\rightarrow\widetilde{\textbf{n}}}^{BV}+\Gamma_{\textbf{n}\rightarrow\widetilde{\textbf{n}}}^{BV}$ and $W_{\textbf{n}\rightarrow\widetilde{\textbf{n}}}^{AB}=\Gamma_{\textbf{n}\rightarrow\widetilde{\textbf{n}}}^{AB}$. Note that, although relocation frequencies do not depend on the configuration before and after the exchange, for the sake of clarity, we choose to follow the notation of the thermal jump frequencies.

We start with a Master Equation expressing the fact that the probability distribution of different configurations is controlled by the transition probabilities between two configurations:

[TABLE]

where $\bm{W}$ is a matrix with components $\bm{W}_{\textbf{n}\widetilde{\textbf{n}}}=W_{\widetilde{\textbf{n}}\rightarrow\textbf{n}}$ if $\textbf{n}\neq\widetilde{\textbf{n}}$ and $\bm{W}_{\textbf{n}\textbf{n}}=-\sum_{\widetilde{\textbf{n}}\neq\textbf{n}}W_{\textbf{n}\rightarrow\widetilde{\textbf{n}}}$. $\bm{P}=(P_{\textbf{n}})$ is a linear vector of probabilities of configurations (n).

For now, the recombination reactions between SIA and V are introduced at the upper scale, within the mean field rate theory model of the average PD concentrations (see Sec. II.4). These athermal events are not treated on the same foot as FAR because SIA and V are considered to be well-separated at the end of the displacement cascade under dilute approximation Nordlund et al. (2018). In this case, the SIA-V recombination requires long-range diffusion, thereby not incorporated in the microscopic Master Equation.

We explain below the method we use to determine the dynamical short range order (SRO) parameters at NESS and the diffusion properties from the Master Equation.

III.1 Dynamical short range order

Starting from the thermal equilibrium state, the mix of thermal jumps and FAR leads to NESS. The latter state is characterized by dynamical SRO parameters which depend on FAR frequencies and thermal jump frequencies. We define them from the configurational probabilities, deduced from a stationary condition applied to the Master Equation (Eq. (11)), also called the global detailed balance condition:

[TABLE]

The solution of Eq. (11) at NESS is noted $\bm{P^{\text{ness}}}=(P_{\textbf{n}}^{\text{ness}})$. The SRO parameter for configuration n is defined as the ratio between the configurational probability $P_{\textbf{n}}^{\text{ness}}$ and the one of the a reference configuration denoted $P_{0}^{\text{ness}}$.

III.2 Transport coefficients

The phenomenological transport coefficients ($\lambda_{\alpha\beta}$) are fundamental parameters to describe the diffusion of chemical species ($\alpha$, $\beta$) in alloy at the macroscopic scale. Fluxes of chemical species ($J_{\alpha}$) are proportional to these coefficients:

[TABLE]

with $\nabla\mu_{\beta}$ the established driving force deviating the system from equilibrium. Starting from NESS, we apply a small gradient of chemical potential and compute the resulting flux of atoms and vacancy. Here we extend the SCMF model to jump mechanisms not obeying the microscopic detailed balance. The SCMF theory was first proposed to study the diffusion process with atomic jumps following the principle of microscopic detailed balance Nastar (2005) but the latter is broken because of FAR. By following the nomenclature of Ref. Nastar (2005), the configuration is defined by a vector n. The latter consists in occupation numbers of all species on all sites i.e. {$n^{\text{A}}_{1}$,$n^{\text{B}}_{1}$,$n^{\text{V}}_{1}$; $n^{\text{A}}_{2}$,$n^{\text{B}}_{2}$,$n^{\text{V}}_{2}$; …}, with $n^{\alpha}_{i}$ equal to one if the site $i$ is occupied by species $\alpha$ and zero if not. The transition from configuration n to $\widetilde{\textbf{n}}$ is realized by thermally activated jumps or FAR, with a total frequency ${W}_{\textbf{n}\rightarrow\widetilde{\textbf{n}}}$. Within the standard SCMF theory in Ref. Nastar (2005), $P_{\textbf{n}}(t)$, the non-equilibrium distribution function of configuration n, is expressed as the product of the equilibrium probability $P^{\text{eq}}_{\textbf{n}}$ and a non-equilibrium contribution. Here we choose the reference state to be NESS, and replace $P^{\text{eq}}_{\textbf{n}}$ by the probability distribution function $P^{\text{ness}}_{\textbf{n}}$:

[TABLE]

The Master Equation (see Eq. (11)) is written for a certain configuration n as

[TABLE]

By applying the global detailed balance condition, i.e. Eq. (12), we obtain a reformulation of the Master Equation:

[TABLE]

Note that the standard SCMF theory relies on the microscopic detailed balance (${W}_{{\widetilde{\textbf{n}}\rightarrow\textbf{n}}}P^{\text{ness}}_{\widetilde{\textbf{n}}}={W}_{\textbf{n}\rightarrow\widetilde{\textbf{n}}}P^{\text{ness}}_{\textbf{n}}$). In that case, it is equivalent to consider the transition probabilities entering or exiting a given configuration. When the microscopic detailed balance is not satisfied the transition frequencies to be retained are the entering configurations. The derivation from the Master Equation (Eq. (16)) of the transport coefficients is similar to the standard SCMF theory in Ref. Nastar et al. (2000); Nastar (2005). It is summarized in Appendix.

III.3 SCMF theory under first shell approximation

Here we focus exclusively on the diffusion properties of a dilute binary model alloy A(B): a host matrix of atoms A containing a single solute atom of species B and a single vacancy. The crystallographic structure is chosen to be a fcc crystal. As explained in Sec. II.2, we consider the vacancy as the only type of PDs. Our purpose is to extend the SCMF theory to include athermal FAR mechanisms. In order to derive analytical transport coefficients, we start with a first shell approximation. This approximation consists in neglecting kinetic coupling and thermodynamic interactions between B and V if the distance between both species is beyond 1NN. FAR-a and FAR-d are restricted to exchanges between 1NN sites only. In such dilute alloy, there are five different atom-vacancy thermal exchange frequencies ($\omega_{i=0,1,2,3,4}$) which we designate after the Lidiard’s nomenclature Lidiard (1955) (see Fig. 1), to which we add 4 FAR-a frequencies ($\Gamma_{i=0,1,3,4}^{AB}$) and 5 FAR-d frequencies ($\Gamma_{i=0,1,2,3,4}^{AV}$ and $\Gamma_{2}^{BV}$). The total B-V exchange frequency is noted $W_{2}^{BV}$. The total A-V and A-B exchanges conserving the 1NN distance between B and V are respectively noted $W_{1}^{AV}$ and $W_{1}^{AB}$. The total A-V and A-B exchanges dissociating the B-V pair are respectively noted $W_{3}^{AV}$ and $W_{3}^{AB}$. The total A-V and A-B exchanges associating the B-V pair are respectively noted $W_{4}^{AV}$ and $W_{4}^{AB}$, and all the other A-V and A-B exchanges far from the solute atom B are respectively noted $W_{0}^{AV}$ and $W_{0}^{AB}$. Here we recall that:

[TABLE]

Within the first shell approximation, two configurational probabilities are considered: $P^{\text{ness}}_{1}$ for the configuration where B and V are located at 1NN and $P^{\text{ness}}_{0}$ for the dissociated configuration where B and V are beyond 1NN. The analytical expression of the 1NN-SRO is given by:

[TABLE]

where $E_{b}$ is the B-V 1NN binding energy which is deduced from the ratio of thermal frequencies, with $\text{exp}\left(E_{b}/k_{B}T\right)=\omega_{4}/\omega_{3}$. Note that $P^{\text{ness}}_{1}/P^{\text{ness}}_{0}$ is a SRO parameter revealing the binding tendency of B and V at the 1NN. As $\Gamma$ increases, the SRO parameter decreases towards $1$. Note that the decrease in the above-threshold situation is larger than what is expected in the sub-threshold situation, just because two relocation frequencies contribute to the decrease in the above-threshold case.

The expressions of the phenomenological coefficients $\lambda_{BV}$, $\lambda_{VB}$, $\lambda_{VV}$ and $\lambda_{BB}$ in a dilute binary fcc alloy are given by

[TABLE]

where $\Lambda=7W_{3}+2W_{1}+2W_{2}^{BV}$, $\Lambda_{3}^{\alpha}=3W_{3}^{A\alpha}-2W_{1}^{A\alpha}-W_{2}^{BV}$ and $\Lambda_{4}^{\alpha}=3W_{4}^{A\alpha}P^{\text{ness}}_{0}/P^{\text{ness}}_{1}-2W_{1}^{A\alpha}-W_{2}^{BV}$ $(\text{for}~{}\alpha=B,V)$, with $W_{i}=W_{i}^{AV}+W_{i}^{AB}$ for $i=0,1,3,4$. $C_{BV}^{p}$ is the concentration of B-V pair at 1NN distance and $C_{V}^{m}$ (resp. $C_{B}^{m}$) is the concentration of isolated V (resp. B). These concentrations can be deduced from the total concentrations of B and V (resp. $C_{B}$ and $C_{V}$) by a low temperature expansion formalism Sykes et al. (1973); Ducastelle (1991); Schuler et al. (2017b):

[TABLE]

with $C_{B}^{0}$, $C_{V}^{0}$ to be obtained by solving the following system of equations:

[TABLE]

where $z^{\text{ness}}=12P^{\text{ness}}_{1}/P^{\text{ness}}_{0}$ is the effective partition function at NESS and $z_{0}=12$.

Note that the term $\Lambda_{m=3,4}^{V}$ (resp. $\Lambda_{m=3,4}^{B}$) is related to the vacancy (resp. solute atom) mobility since it contains all the vacancy (resp. solute atom) jump mechanisms including A-V (resp. A-B) and B-V (resp. V-B) exchanges. At equilibrium, $\Lambda_{3}^{V}=\Lambda_{4}^{V}$ and $\Lambda_{3}^{B}=\Lambda_{4}^{B}$ due to the microscopic detailed balance. Hence the two off-diagonal equilibrium coefficients $\lambda_{VB}$ and $\lambda_{BV}$ are equal, according to the Onsager reciprocal relations. In addition, $\lambda_{VV}$ (resp. $\lambda_{BB}$) can be separated into two parts: $C_{V}^{m}W_{0}^{AV}$ (resp. $C_{B}^{m}W_{0}^{AB}$) and the rest. The latter represents the exchanges of the B-V pair at 1NN distance while the former represents the hops of the isolated V (resp. B).

In the case of sub-threshold irradiation for which there is no direct exchange between atoms (i.e. $\Gamma^{\text{aa}}=0$), the off-diagonal coefficients are equal and, from Eq. (19), (20) we get:

[TABLE]

Although the microscopic detailed balance is broken for the individual exchange frequencies $\omega_{3}$ and $\Gamma^{\text{ad}}$, it still holds for the sum of the latter frequencies, that is $W_{3}=W_{3}^{AV}=\omega_{3}+\Gamma^{\text{ad}}$ (see Eq. (18)).

By replacing the total frequencies by the corresponding thermally activated jump frequencies, and replacing the dynamical SRO at NESS by the equilibrium SRO, the transport coefficients turn out to be equivalent to the Onsager coefficients $L_{BV}$ of the five-frequency model within the first shell approximation Howard and Lidiard (1964).

The variation of $\lambda_{BV}$ with $\Gamma^{\text{ad}}$ depends on the full set values of the thermal-activated jump frequencies. When $\Gamma^{\text{ad}}$ is dominant before all $\omega_{i}$: $\lambda_{BV}\sim-C_{B}C_{V}\Gamma^{\text{ad}}$. Note that if $13\omega_{3}>2\omega_{1}$ ($L_{BV}<0$), then $\lambda_{BV}$ remains negative whatever the magnitude of the relocation frequencies. Otherwise, a change of sign of $\lambda_{BV}$ can be observed when $\Gamma^{\text{ad}}\simeq-(13\omega_{3}-2\omega_{1})/11$. Therefore, when a solute atom is dragged by a vacancy, FAR-d may change the sign of the solute-vacancy flux coupling and destroy the solute drag effect. In the opposite case, when $L_{BV}$ is negative, FAR-d does not change the sign of the solute-vacancy flux coupling.

We consider now the case of above-threshold irradiation. Then FAR have two contributions: FAR-a and FAR-d, with $\Gamma^{\text{aa}}=\Gamma$ and $\Gamma^{\text{ad}}=\gamma\Gamma$. The off-diagonal terms $\lambda_{BV}$ and $\lambda_{VB}$ are not equal and their difference $\Delta\lambda=\lambda_{VB}-\lambda_{BV}$ is given by:

[TABLE]

with

[TABLE]

Note that $\Delta\lambda=0$ in the two extreme cases when thermal jumps ($\omega$) are dominant (i.e. $\Gamma/\omega\rightarrow 0$) or negligible ($\Gamma/\omega\rightarrow\infty$). The sign of $\Delta\lambda$ is determined by the product $(1-\omega_{3}/\omega_{4})[\omega_{a}-(1-\gamma)\Gamma]$. If $\gamma=1$, this product involves thermal jump frequencies only. The first parenthesis is directly related to the equilibrium SRO parameter: $(1-\omega_{3}/\omega_{4})$ is positive if the vacancy and the solute atom attract each other and negative otherwise. The higher the thermodynamic attraction, the smaller the ratio $\omega_{3}/\omega_{4}$, and the larger the difference $\Delta\lambda$.

III.4 Extension of the KineCluE code

For a more precise calculation beyond the first shell approximation, we consider each configuration where V and B are located at a distance lower than the kinetic radius $R_{k}$, as B-V pair configuration. At distances larger than $R_{k}$, B and V are considered as isolated monomers. Therefore, 3 cluster contributions are included: monomer B, monomer V as well as B-V pair. Note that the calculation under first shell approximation performed at Sec. III.3 is a particular situation where the kinetic radius is set equal to the 1NN distance. The calculation of the cluster transport coefficients is performed using the KineCluE code Schuler et al. (2020). The latter accounts for all the kinetic paths within a pair cluster defined by radius $R_{k}$. Note that the kinetic radius can be set well beyond the 1NN distance in KineCluE. This allows us to perform converged calculation of cluster transport coefficients including long-distance FAR as well as long range kinetic correlations. In order to use NESS as reference state, a module is added to the code which calculates the NESS probability distribution by solving Eq. (11). Besides, the underlying principle of the microscopic detailed balance of the code is replaced by the global detailed balance condition (Eq. (12)). Detailed descriptions will be published elsewhere. Models 1, 2 and 3 presented in Section II.2 have been introduced into KineCluE. Note that the cluster radius $R_{c}$ in Model 3 is set equal to $R_{k}$ for simplicity.

III.5 Comparison between KineCluE results and Monte Carlo simulations

As mentioned in the introduction (Sec. I), as soon as one of the microscopic diffusion mechanisms ($W_{\textbf{n}\rightarrow\widetilde{\textbf{n}}}^{AV}$ and $W_{\textbf{n}\rightarrow\widetilde{\textbf{n}}}^{AB}$) does not obey the microscopic detailed balance, we cannot use the Allnatt formulae Allnatt (1965, 1982) to extract the phenomenological transport coefficients from atomistic kinetic Monte Carlo (AKMC) simulations. However, for a binary alloy with solute-point defect interactions restricted to 1NN pairwise interactions, we have shown in Sec. III.3 that detailed balance is fulfilled in the case of sub-threshold irradiation. Therefore, in this specific case, we may rely on the Allnatt formulae to obtain the Onsager matrix of the transport coefficients. As for the thermodynamic properties, we may apply AKMC to study the dynamical short range order characterizing a NESS from an average on the residential time, relying on the ergodic principle.

We choose here a model alloy with highly attractive vacancy-solute interactions because it emphasizes the effect of FAR on flux coupling. The migration barriers (in eV) are set to 0.95 for $\omega_{0}$ and $\omega_{3}$, 0.75 for $\omega_{1}$ and $\omega_{4}$, and 0.60 for $\omega_{2}$. The attempt frequency $\nu$ is chosen to be $10^{14}~{}\text{s}^{-1}$. As for the model of FAR, we choose Model 1, with $r_{m}$ equal to the 1NN distance $r_{1}$.

The AKMC simulation box is a fcc crystal of 2048 sites. It contains one single solute atom and one vacancy. We apply periodic boundary conditions and use a residence-time algorithm. At each Monte Carlo step, we propose the whole set of the thermal jumps and FAR. We select one exchange from the proposed mechanisms. After every exchange, we compute the residence time increment. From the fluctuations of atomic positions, we compute the transport coefficients. Note that the corresponding off-diagonal coefficients given by the AKMC method are by construction symmetric. As shown in Ref. Gallavotti (1996); Lau et al. (2007), they do not correspond to the transport coefficients $\lambda_{BV}$ and $\lambda_{VB}$ whenever one of the diffusion mechanism does not obey the detailed balance.

As for the KineCluE approach, the kinetic radius $R_{k}$ of the cluster B-V is set to $4a_{0}$.

Fig. 2 shows the evolution of the dynamical 1NN-SRO under sub- and above-threshold FAR. We obtain an excellent agreement between KineCluE and AKMC simulations on the SRO paramters. As expected, the dynamical SRO decreases with the relocation frequency with a higher rate in the case of an above-threshold irradiation.

Fig. 3 shows the variation of the transport coefficients with the frequency of FAR-d in the sub-threshold irradiation regime. Both KineCluE and AKMC methods give the same transport coefficients because the microscopic detailed balance holds for the total transition rates. However, when $\Gamma$ is small compared with thermal jump frequencies, we observe a slight discrepancy between the coefficients. Yet the size of the AKMC simulation box is comparable with $R_{k}$. The discrepancy is due to the difference in the applied boundary conditions between KineCluE and the AKMC method. In KineCluE, configurations of solute and vacancy located at a distance larger than the kinetic radius are not included in the calculation, while the AKMC method relies on periodic boundary conditions. In the latter, atoms or PDs exiting from the simulation box enter back through another side and add a kinetic correlation contribution to the transport coefficients.

In the case of an above-threshold irradiation, we observe in Fig. 4 a similar behaviour of the diagonal transport coefficient $\lambda_{BB}$, whereas the single off-diagonal coefficient measured in AKMC simulations does not correspond any more to the off-diagonal phenomenological transport coefficients obtained by KineCluE.

IV Results on diffusion properties

Here we focus on the above-threshold irradiation case. We consider a model alloy with relatively high migration barriers. Hence the alloy is potentially sensitive to FAR effects, just because the thermal jump frequencies are small with respect to the relocation frequency deduced from realistic dose rate. The energy interaction between B and V is restricted to a pairwise 1NN interaction. The migration barriers (in eV) are set to 1.10 for $\omega_{0}$, $\omega_{1}$ and $\omega_{3}$, 0.90 for $\omega_{4}$, and 0.80 for $\omega_{2}$. The attempt frequency $\nu$ is chosen to be equal to $5\times 10^{12}~{}\text{s}^{-1}$. The three relocation models indicated in the Section II.2 are considered. We use KineCluE to calculate the transport coefficients.

The parameter values that we set to estimate the vacancy concentration under irradiation are shown in Table 1. Here the mean relocation range $r_{m}$ and the cut-off distance $L$ for Models 2 and 3 are respectively set to 1NN ($\sqrt{1/2}\,a_{0}$) and 5NN ($\sqrt{5/2}\,a_{0}$) distances. The kinetic radius is set to $2a_{0}$.

IV.1 Dynamical short range order

Fig. 5 shows the profile of dynamical SRO as a function of $\Gamma$ for Models 1, 2 and 3. The probability for B and V at 1NN distance is reduced by FAR leading to an effective B-V interaction smaller than the thermodynamic one. The decrease of 1NN-SRO with the relocation frequency in Model 1 starts when $\Gamma$ is around $10^{-2}\,\text{s}^{-1}$. The decrease starts earlier in Models 2 and 3: respectively around $10^{-4}\,\text{s}^{-1}$ and $10^{-3}\,\text{s}^{-1}$. However, the 1NN-SRO of Model 3 converges towards non-zero value at large $\Gamma$. In Model 1, there is no interaction between B and V beyond the 1NN distance, whatever the relocation frequency. However, in Models 2 and 3, we observe that the effective B-V interaction extends beyond the range of the thermal one (i.e. beyond the 1NN). The effective interaction remains up to 5NN distance when $\Gamma$ is comparable to one of the thermal jump frequencies. This is due to the relatively long range FAR. In the extreme case when $\Gamma$ is dominant before the thermal jump frequencies, the B-V interactions are dropping in Models 1 and 2 whereas in Model 3 the 1NN attraction is slightly decreasing and the 2NN, 3NN, 4NN and 5NN are slightly increasing. The binding tendency of a vacancy around the solute atom is still high ($P_{1}^{\text{ness}}/P_{0}^{\text{ness}}\simeq 10^{2}$) due to the introduction of the biased FAR-d with the 1NN atoms of the solute atom in Model 3.

IV.2 Tracer diffusion coefficient

In the dilute limit, the tracer diffusion coefficient of solute B is written as

[TABLE]

Phenomenological models of diffusion under irradiation systematically rely on the assumption that the thermally activated diffusion and FAR take place in parallel Müller et al. (1988); Roussel and Bellon (2002). The tracer diffusion coefficient is then written as a sum of two diffusion coefficients:

[TABLE]

where $D_{B,\text{th}}^{*}$ is the thermal diffusion coefficient commonly deduced from diffusion experiments or atomic based diffusion models and $D_{B,\text{far}}^{*}$ is the diffusion coefficient of solute atom B resulting from FAR only. Note that both coefficients can be calculated by KineCluE. Unless one diffusion mechanism is dominant over the other, we expect a non-additive contribution to the solute tracer diffusion coefficient because of the kinetic correlations. In order to quantify the non-additive contribution, we define the parameter:

[TABLE]

Fig 6 shows the variation of the solute diffusion coefficient with the relocation frequency. We observe that the global tendencies of the diffusion coefficients obtained with the three models are similar. However, the three curves do not have the same asymptote at large $\Gamma$. The largest difference occurs when the correlation factor $f_{B}$ is increased by FAR. With Models 1 and 2, this factor tends to 1 when $\Gamma$ is dominant over the thermal jump frequencies, meaning that there are no kinetic correlations. However, in Model 3, the correlation factor tends to 0.46. The remaining kinetic correlations are due to the biased FAR-d. Besides, $\Delta D_{B}$ is high when $\Gamma$ is in the range of the thermal jump frequencies because then, there is a strong competition between the thermal mechanisms and FAR. In this example, $\Delta D_{B}$ spans from 100 % to 300 % depending on the relocation model.

IV.3 Flux coupling

We characterize the flux coupling between solute B and vacancy V by computing the wind factors Anthony (1969, 1970); Okamoto (1979)

[TABLE]

and

[TABLE]

Both wind factors describe the B-V flux coupling related to two different situations. The wind factor $\delta_{B\rightarrow V}$ gives the number of solute atoms following a vacancy under the driving force $\nabla\mu_{V}$ and the wind factor $\delta_{V\rightarrow B}$ indicates the number of vacancies dragged by a solute atom under the driving force $\nabla\mu_{B}$. If the wind factors are positive, a drag of B by V (or vice versa) may occur. As shown in Sec. IV.1, the interactions between the solute atom and the vacancy are reduced or even destroyed by FAR. Since the drag effect is highly related to this interaction, we study the effect of the relocation frequency $\Gamma$ on the wind factors.

Fig. 7 shows the variation of the wind factors with the relocation frequency. Whatever the relocation models, $\delta_{B\rightarrow V}$ and $\delta_{V\rightarrow B}$ globally decrease with $\Gamma$. However, $\delta_{V\rightarrow B}$ of Model 1 has a surprising non-monotonous behaviour: the drag effect is enhanced before being destroyed by FAR. $\delta_{B\rightarrow V}$ of Model 3 has also an atypical behaviour: it slightly increases and tends to a non-zero value at large $\Gamma$, meaning that the solute drag and vacancy drag effects are not totally destroyed. This is because the biased FAR-d maintains a flux coupling between B and V. This persistent flux coupling at high radiation flux should be very sensitive to the details of the relocation mechanism.

V Sensitivity study with respect to the model and alloy parameters

FAR models depend on the values of the mean relocation range $r_{m}$, the kinetic radius $R_{k}$ of the cluster B-V and the truncation distance $L$. However, the latter parameter is not a physical parameter. Since the relocation frequency exponentially decreases with the distance between B and V (see Eq. (5)), the value of $L$ does not affect the diffusion properties as long as it is large enough. Therefore, we focus here on the sensitivity of the results to the other two parameters: $r_{m}$ and $R_{k}$.

V.1 Kinetic radius

In general, the results given by KineCluE code converge with the kinetic radius $R_{k}$ Schuler et al. (2020). However, because $R_{c}=R_{k}$ in Model 3, the FAR-d models for a monomer vacancy and for a vacancy within the B-V pair are different. In this case, the results obtained with Model 3 may depend on the values of $R_{k}$. However, Fig. 8 shows that $D_{B}^{*}$, $\Delta D_{B}$ and $\delta_{B\rightarrow V}$ are not very sensitive to the change of the kinetic radius. Although, the decrease rate of $\delta_{V\rightarrow B}$ with $\Gamma$ is slower with $R_{k}=3a_{0}$ than $2a_{0}$. This is because the vacancy performs biased FAR-d with the 1NN atoms of the solute atom from longer distances.

V.2 Mean relocation range

Fig. 9 shows the effect of the mean relocation distance $r_{m}$ on the solute diffusion and flux coupling. First we focus on Model 1. Since the solute mobility is enhanced when increasing the relocation distance, the corresponding solute diffusion coefficient increases with $r_{m}$. Besides, according to the plot of $\Delta D_{B}$, the interaction between thermal jumps and FAR decreases with $r_{m}$. The thermally-activated jump distance and the thermal interaction between B and V are both restricted to 1NN. The larger the relocation distance, the smaller the B-V interaction. Thus B and V are more likely to diffuse as monomers, a kinetic regime where the diffusion properties related to the thermal jumps and FAR become additive. As for the flux coupling, the decreasing rate of $\delta_{B\rightarrow V}$ with $\Gamma$ increases with $r_{m}$. Thus the solute drag effect is destroyed more easily. Besides, the variation tendency of $\delta_{V\rightarrow B}$ with $\Gamma$ become qualitatively different when $r_{m}>$ 1NN. The vacancy drag effect is not enhanced when $r_{m}$ equals to 2NN and 3NN. This may be due to the same reason mentioned before: B and V have many more paths to escape from each other. As for the results obtained with Models 2 and 3, they have similar profiles as the ones in Model 1.

V.3 Atomic mixing rate

Since the previous section has shown that the effects of FAR are roughly the same in terms of the global tendency whatever the relocation model and the mean relocation distance, we choose the simplest model, Model 1 with $r_{m}=r_{1}$. As stated in Sec. II.3, the number of relocations per Frenkel pair created (i.e. $n_{\text{rel}}$) should be alloy specific due to the thermal effect on heat spike mixing. Fig. 10 shows the variation of the diffusion properties in function of radiation dose rate with different values of $n_{\text{rel}}$. The effect of FAR on the flux coupling and tracer diffusion occurs at a smaller dose rate when $n_{\text{rel}}$ increases. Moreover, we observe that $\Delta D_{B}$ decreases with $n_{\text{rel}}$. These results show the importance of $n_{\text{rel}}$ in the prediction of a critical dose rate when the effects of FAR on the flux coupling and tracer diffusion is paramount.

V.4 FAR-d frequencies

For the reason stated in Sec. V.3, Model 1 with $r_{m}=r_{1}$ is chosen for the sake of simplicity. Note that FAR-a and FAR-d are due to different phenomena in the displacement cascade: FAR-a describes the recoil mixing due to PKA while FAR-d models the lattice site change of PD during the quenching process. There is no guarantee that the frequencies of FAR-a ($\Gamma^{aa}$) and FAR-d ($\Gamma^{ad}$) are equal. Fig. 11 shows the plot of $D_{B}$, $\delta_{B\rightarrow V}$ and $\delta_{V\rightarrow B}$ as a function of relocation frequency $\Gamma^{aa}=\Gamma$ with different ratios $\gamma=\Gamma^{ad}/\Gamma^{aa}$. The global tendencies of the above quantities are not affected by the variation of the ratio $\gamma$. Besides, the tracer diffusion coefficient $D_{B}$ is not sensitive to the variation of the ratio $\gamma$. However, $\delta_{B\rightarrow V}$ decreases with $\gamma$ while the variation of $\delta_{V\rightarrow B}$ has the opposite tendency. For $\gamma\neq 1$, FAR-a and FAR-d effects on the solute atom and PD diffusion occurs at different dose rate. It respectively happens when $\Gamma^{aa}$ (i.e. $\Gamma$) and $\Gamma^{ad}$ (i.e. $\gamma\Gamma$) are of the same order of magnitude compared with thermal jump frequencies. In brief, the smaller the $\gamma$ value, the larger the difference between the frequencies for FAR-a and FAR-d, and the more important the strength of the flux coupling.

V.5 Thermal jump parameters

The effects of FAR depend on the radiation dose rate and the intrinsic thermal jump frequencies of the alloy. We use KineCluE to perform a sensitive study of the radiation kinetic properties with respect to the thermal jump frequencies. Fig. 12 shows the variation of $\Delta D_{B}$ and the wind factors $\delta_{B\rightarrow V}$, $\delta_{V\rightarrow B}$ with respect to $\Gamma$, for various values of $\omega_{4}$. The values of the other thermal jump frequencies are fixed. Model 1, with $r_{m}=r_{1}$, is chosen for the following discussion. The interactions between thermal jumps and FAR are emphasized in this case because the hop distances are both 1NN. The ratio $\omega_{4}/\omega_{3}$ directly affects the binding energy $E_{b}$ between solute atom and vacancy at 1NN. We observe that $\Delta D_{B}$ and wind factors increase with the binding energy. Besides, the larger the binding energy, the larger the enhancement of the wind factor $\delta_{V\rightarrow B}$ by FAR. This can be explained by noting that the solute atom and vacancy tend to be closer to each other with a larger binding energy. Therefore, the interaction between FAR and thermally activated diffusion of solute atom is more important, leading to a larger difference from what we would except with an additive model, i.e. Eq. (27). Moreover, the binding tendency of the vacancy and the solute atom increases, causing an enhancement of the wind factor $\delta_{V}$. As well, $\omega_{1}$ and $\omega_{2}$ have a non-negligible effect on the profile of $\Delta D_{B}$ and wind factors in function of $\Gamma$. Here we set $\omega_{4}$ to its initial value $2.3\times 10^{1}s^{-1}$ and we perform calculations with different values of $\omega_{1}$ and $\omega_{2}$. Fig. 13 shows that if $\omega_{2}$ is large compared to $\omega_{1}$ (more than 1 order of magnitude), $\Delta D_{B}$ and $\delta_{B\rightarrow V}$ increase with $\omega_{1}$ whereas the enhancement of $\delta_{V\rightarrow B}$ by FAR decreases with $\omega_{1}$. If the amplitudes of $\omega_{2}$ and $\omega_{1}$ are comparable (within 1 order of magnitude), the trends are opposite: $\Delta D_{B}$ and $\delta_{B\rightarrow V}$ decrease with $\omega_{1}$ whereas the enhancement of $\delta_{V\rightarrow B}$ by FAR increases with $\omega_{1}$. However, we observe that if the values of $\omega_{1}$ and $\omega_{2}$ are close (within 1 order of magnitude), the variations of $\Delta D_{B}$ and wind factors with $\Gamma$ are not sensitive to $\omega_{1}$.

VI Discussion and Summary

Neutron or ion irradiation in metals generates displacement cascades. We present a simplified model of this complex phenomenon by introducing FAR mechanisms, and an average creation rate of PD uniform in time and space. To calculate the energetic and kinetic properties, we write a Master Equation for the evolution of the distribution function which includes both the thermal jumps and FAR. We extend the SCMF theory to solve and compute the SRO parameters and the phenomenological transport coefficients at the NESS reached under irradiation. The main difficulty lies in the loss of the microscopic detailed balance when considering FAR mechanisms. Relying on Model 1 including FAR between 1NN sites only and a first shell approximation of the kinetic correlations, we derive analytical expressions of the phenomenological transport coefficients. We demonstrate that FAR does not produce a simple additive term to the transport coefficients. When the magnitude of the relocation frequency is in the range of the thermal frequencies, FAR interacts with the thermal diffusion mechanism, yielding non-symmetric off-diagonal transport coefficients and a solute tracer diffusion coefficient deviating from a direct sum of the contributions of thermal jumps and FAR. This deviation increases with the solute kinetic correlations. We use the automated code KineCluE to yield a more systematic study of the effect of the range and magnitude of FAR on the kinetic properties, including a sensitivity study with respect to the alloy thermodynamics and the models of relocation and PD production.

Due to the lack of data on the detailed mechanisms of FAR and PD production, we introduce Models 2 and 3 representing two extreme situations, expecting the real situation to be in-between. In Model 2, we assume that FAR is a fully random process while in Model 3, we introduce a biased FAR-d with the 1NN atoms of the solute atoms to reproduce the fact that the vacancy creation within a cascade is partially driven by the vacancy-solute thermodynamic attraction. As a result, part of the vacancy-solute SRO remains which in turns leads to a higher resistance of the vacancy-solute flux coupling to irradiation. Positive flux coupling is the result of strong kinetic correlations, which can be modified by introducing FAR mechanisms. Our sensitivity study shows that the magnitude of the surviving kinetic correlations strongly depends on the details of the biased FAR-d mechanism, while the reduction of correlations and flux coupling due to the randomizing processes are less sensitive to the details of the relocation events unless the distance of FAR is close to the thermodynamic range. A persistent vacancy-solute flux coupling at low temperature and high radiation flux may play an important role on the solute redistribution in irradiated materials. Therefore, the mechanism of PD production with respect to the solute atom spatial distribution within the displacement cascade should be analyzed more precisely.

Eventually, the effect of the interplay between thermal jumps and FAR on vacancy-solute positive flux coupling is important when the solute-vacancy thermodynamic attraction is large, the magnitude of the thermal jump frequencies compared with the relocation frequency and the range of thermodynamic interactions is close to the relocation distances. As for the tracer diffusion coefficients, their non-additivity property with respect to FAR and thermal jumps follows the same trend as the flux coupling phenomena in systems featuring positive flux coupling but may also arise in case of no positive flux coupling but strong correlated solute migration paths. For instance, the additive expression of Eq. (27) reproduces correctly the diffusion coefficient of Au in Al measured under irradiation Acker et al. (1974). This is because the vacancy-jump barrier in Al is around 0.58 eV Wu et al. (2016), hence the thermal jump frequencies are dominant over the relocation frequencies under realistic experimental conditions. However, we expect a non-negligible effect of FAR in Ni-based alloys because the vacancy-mediated migration barrier in pure Ni is high (around 1.09 eV Wu et al. (2016)).

Appendix A SCMF expressions at NESS

In this appendix, we introduce in short how to obtain the transport coefficients from the microscopic Master Equation by the standard SCMF formulation Nastar et al. (2000); Nastar (2005). Note that $\delta P_{\textbf{n}}(t)$ in Eq. (16) is a corrective term representing the modification of the effective distribution function $P^{\text{ness}}_{\textbf{n}}$ due to the presence of an applied driving force. It is written as

[TABLE]

where $\beta=1/k_{B}T$, $\delta\Omega$ is the normalization factor, $\delta\mu_{i}^{\alpha}$ is the deviation from the stationary-state chemical potential on site $i$ of the chemical species $\alpha$ compared to the bulk atom, and $h$ is the time-dependent effective Hamiltonian restricted to the pair interaction written as

[TABLE]

where $\nu_{ij}^{\alpha\beta}(t)$ is the time-dependent effective pair interactions. The latter can be determined by solving the kinetic equations deduced from the Master Equation. Here, Eq. (31) is linearized with respect to the terms $\beta\delta\mu_{i}^{\alpha}$ and $\beta h$ because we are close to NESS:

[TABLE]

Starting from the Master Equation Eq. (16), the time derivative of the ensemble average can be given by

[TABLE]

where $\widetilde{n}_{i}^{\alpha}$ is the occupation number of the configuration $\widetilde{\textbf{n}}$. By denoting $\langle\cdot\rangle^{\text{ness}}$ the ensemble average over the distribution function at NESS (i.e. $P_{\textbf{n}}^{\text{ness}}$), the derivative of the one-point average $\langle n_{i}^{\alpha}\rangle$ can be given by

[TABLE]

where $\delta\mu_{i}^{\alpha\gamma}=\delta\mu_{i}^{\alpha}-\delta\mu_{i}^{\gamma}$. By applying the continuity equation to the kinetic equation of the one-point average written as:

[TABLE]

we can deduce the expression of the flux of chemical species. The fluxes of solute atom and vacancy under first shell approximation are recognized to be:

[TABLE]

where $\hat{\nu}_{ij}^{BV}=\nu_{ij}^{BV}-\nu_{ij}^{BA}-\nu_{ij}^{AV}$.

Note that under the first shell approximation, the effective interactions are restricted to the pair B-V at 1NN. The term $\hat{\nu}^{BV}$ can be estimated from the stationary condition of the kinetic equation of the two-point average $\left<n_{i}^{B}n_{j}^{V}\right>$. As a result, $\hat{\nu}^{BV}$ can be expressed as a function of the chemical potential gradient. Therefore, the atomic fluxes of solute atom and vacancy are also functions of ${\nabla}\mu^{VA}$ and ${\nabla}\mu^{BA}$, allowing us to identify the transport coefficients and the expressions under first shell approximation are given in Eq. (19)–(III.3).