A modified formula for non-Arrhenius diffusion of helium in metals

TL;DR

This paper introduces a modified diffusion formula for helium in metals, improving the accuracy of kinetic Monte Carlo simulations by accounting for low-energy barriers in helium diffusion.

Contribution

A new stochastic model based on Brownian motion and cosine potential is derived to better predict helium diffusion rates in metals, surpassing traditional Arrhenius law limitations.

Findings

Model predictions align with molecular dynamics simulations.

The approach improves accuracy for low-energy barrier diffusion.

Enhances KMC simulation of helium-vacancy cluster growth.

Abstract

Helium diffusion in metals is the basic requirement of nucleation and growth of bubble, which gives rise to adverse degradation effects on mechanical properties of structural materials in reactors under irradiation. Lattice based Kinetic Monte Carlo approach is widely adopted to study the evolution of helium-vacancy clustering. However, the imple-mentation of Arrhenius law to prediction the event rate of single interstitial helium solute diffusion in metal is not always appropriate due to low-energy barrier. Based on a stochastic model, a modified formula is derived from the Brownian motion upon a cosine-type potential. Using the parameters obtained from molecular dynamics simulation for the diffusivity of single helium solute in BCC W, the prediction of our model is consistent with the results from dynamical simulation and previous model. This work would help to develop a more accurate KMC scheme for the growth of helium-vacancy clusters, as well as other low-energy reactions in materials science.

Figures

Equations

$$k\text{He}+l\text{V}+{\text{He}}_m{\text{V}}_n\rightleftharpoons {\text{He}}_{m+k}{\text{V}}_{n+l}$$

$$\nu ={\nu }_0\exp (-E_m/k_{\text{B}}T)$$

$$D=g{\lambda }^2\nu =D_0\exp (-E_m/k_{\text{B}}T)$$

$$D=\mu k_{\text{B}}T=\frac{k_{\text{B}}T}{m^{\ast }\gamma }\propto T$$

$$m^{\ast }\ddot{X}=F-m^{\ast }\gamma \dot{X}+f(t)$$

$$D=\frac{⟨{\left[X(t)-X(0)\right]}^2⟩}{2t}=\frac{m^{\ast }}{k_{\text{B}}T}{\int }_0^{\infty }⟨\dot{X}(t)\dot{X}(0)⟩\text{d}t$$

$$D=\{\begin{array}{rlr}& D_0{e}^{-E_m/k_{\text{B}}T},& \text{if}{E}_m\gg k_{\text{B}}T\\ & \mu k_{\text{B}}T,& \text{if}{E}_m\ll k_{\text{B}}T\end{array}$$

$$D=\mu k_{\text{B}}T{\left[{\int }_0^{\lambda }{e}^{\frac{E(X)}{k_{\text{B}}T}}\text{d}X{\int }_0^{\lambda }{e}^{-\frac{E(X)}{k_{\text{B}}T}}\text{d}X\right]}^{-1},\text{if}\gamma \gg {\nu }_0$$

$$D=\frac{\pi k_{\text{B}}T}{\gamma E_m}\exp (-E_m/k_{\text{B}}T),\text{if}\gamma \ll {\nu }_0$$

$$E(X)=\frac{E_m}{2}\left[1-\cos (2\pi X/\lambda )\right]$$

$$4{\pi }^2{\nu }_0^2=\frac{1}{m^{\ast }}{\left(\frac{{\partial }^{2}E}{\partial X^2}\right)}_{X=0}\Rightarrow 2\pi {\nu }_0=\frac{\pi E_m}{m^{\ast }{\nu }_0{\lambda }^2}=\frac{\pi g{E}_m}{m^{\ast }{D}_0}$$

$$b=g\frac{\gamma }{2\pi {\nu }_0}=\frac{g}{m^{\ast }\mu }\frac{m^{\ast }{D}_0}{\pi g{E}_m}=\frac{D_0}{\pi \mu E_m}$$

$$D=\mu k_{\text{B}}T{\left[J_0\left(\frac{E_m}{2k_{\text{B}}T}\right)\right]}^{-2}$$

$$D=\{\begin{array}{rlr}& \pi \mu E_m{e}^{-E_m/k_{\text{B}}T},& \text{if}{E}_m\gg k_{\text{B}}T\\ & \mu k_{\text{B}}T,& \text{if}{E}_m\ll k_{\text{B}}T\end{array}$$

$$D$$

$$=D_0\frac{k_{\text{B}}T}{\pi E_m}{\left\{{\left[J_0\left(\frac{E_m}{2k_{\text{B}}T}\right)-1\right]}^2+\frac{D_0}{\pi \mu E_m}\right\}}^{-1}$$

$$D=\{\begin{array}{rlr}& D_0{e}^{-E_m/k_{\text{B}}T},& \text{if}{E}_m\gg k_{\text{B}}T\\ & \mu k_{\text{B}}T,& \text{if}{E}_m\ll k_{\text{B}}T\end{array}$$

$$\nu =\frac{D}{g{\lambda }^2}={\nu }_0\frac{k_{\text{B}}T}{\pi E_m}{\left\{{\left[J_0\left(\frac{E_m}{2k_{\text{B}}T}\right)-1\right]}^2+\frac{D_0}{\pi \mu E_m}\right\}}^{-1}$$

$$\epsilon (T)=\frac{\left|D_A(T)-D_E(T)\right|}{D_E(T)}\times 100\%$$

Full text

1 Introduction

Helium (He) is one of the most common productions in structural materials of fission and fusion reactors under irradiation. Due to the insolubility of He atom in metals, it would be easily trapped into sinks, such as vacancies and grain-boundary, and would form the helium-bubble in a long-term evolution process. The formation and accumulation of helium plays important role in the high-temperature kinetics of microstructure evolution of nuclear materials under irradiation, leading to the adverse ageing effects, e.g. high-temperature helium embrittlement [1, 2]. Therefore, the mechanistic understanding of bubble evolution is considered as one of the key issues in nuclear materials science and engineering [3]. In this paper, we concentrate on the non-Arrhenius diffusion behavior of helium in metals.

As long-term phenomena, the growth of helium bubbles in metals are made up by numerous atomic activation processes, including the migration of defects, the combination of helium and vacancies, and the dissociation of helium-vacancy clusters [4]. In theoretical research, a well-developed lattice based kinetic Monte Carlo (short for “KMC” in the rest of paper) approach [5, 6] is widely adopted, by treating the evolution of helium-vacancy clusters as a chemical-reaction-like process, e.g.,

[TABLE]

as shown in Fig. 1. Further, each individual atomic process is treated as a classical reaction event, whose occurring frequency $\nu$, i.e., event rate, is associated with the reaction-path and the corresponding energy barrier $E_{m}$, usually described by Arrhenius equation in KMC [7], as

[TABLE]

with $\nu_{0}$ the attempt frequency, $k_{\textrm{B}}$ the Boltzmann constant, and $T$ the absolute temperature.

Note that, because of the repulsive interaction between helium and metal atoms and a small migratory energy barrier, the forward reaction process in Eq. (1) seems to be diffusion-controlled at finite temperatures (See Fig. 1 of Ref. [8] or the schematics in Fig. 1), that the growth rate of a helium-vacancy cluster is mainly determined by the mass-flux of helium towards the trap [9], equivalent to the diffusion behavior. In the present KMC scheme, the diffusivity $D$ is usually described by an Arrhenius equation, as

[TABLE]

where $D_{0}=g\nu_{0}\lambda^{2}$ is a pre-factor, with $g$ a geometrical factor and $\lambda$ the distance between two adjunct stable sites in real-space.

For the thermal assistant diffusion event, the prerequisite of using Arrhenius equation in Eq. (3) to estimate the diffusivity is that the diffusion should be treated as a thermodynamic reversible process, where the energy barrier between two adjunct stable sites along the reaction path should be far larger than the thermal energy, i.e., $E_{m}\gg k_{\textrm{B}}T$. Equivalently, the mean time for the defect escaping out of the trap should be far larger than the characteristic relaxation time for the defect-system from saddle point configuration towards equilibrium state. For the helium-vacancy clusters or mono-vacancy, the migration energy is in the order of $\sim$1eV, which satisfies this prerequisite in a wide range of temperature, e.g., $T<10^{4}$K. In this case, the diffusion behavior could be well-described by Arrhenius equation [3].

However, calculations based on density functional theory (DFT) [10] indicated that $E_{m}\sim$ 0.1eV of single interstitial helium in metals, therefore, the migratory energy barrier is comparable with (even larger than), the thermal energy at $T>1000$K, breaking down the condition of $E_{m}\gg k_{\textrm{B}}T$. In this scenario, it will reveal anomalous diffusion behavior beyond the Arrhenius description [11], as demonstrated by the recent molecular dynamics (MD) simulations [12, 13, 14] of helium diffusion in BCC W ($E_{m}\sim 0.15$eV). More precisely, the real-space trajectory of interstitial helium at 2600K (See Fig. 2(b) in Ref. [14]) indicated that it seems not a hopping random walker with a long resident time between its stable sites, but a thermal assisted Brownian motion in a stochastic media. Therefore, the high-temperature diffusivity of single interstitial helium in BCC W does not obey Arrhenius equation but the Einstein-Smoluchowski equation, i.e.,

[TABLE]

where $\mu=1/m^{*}\gamma$ is the classical mobility, with $m^{*}$ the effective mass of the diffusing object, and $\gamma$ the friction coefficient provided by the host media for the diffusion.

In fact, the diffusion behavior of defects having a small transition barrier, e.g., $E_{m}\leqslant 0.1$eV, is categorized to be the low-energy reactions [15], which is found to be a common phenomenon in other disciplines, such as adatom migration upon crystal surface [16, 17, 18], self-interstitial diffusion [11] and dislocation motion [19] in metals. Following Kramers’ theory [20], theoretical studies [21, 22] discussed the physical picture of these low-energy reactions: (1) the low-temperature behavior is governed by Arrhenius equation at $E_{m}\gg k_{\textrm{B}}T$, which the quasi-equilibrium properties, i.e., $E_{m}$ and $\nu_{0}$ shown in Eq. (2), are used to described the reversible nature; (2) the high-temperature behavior is governed by Einstein-Smoluchowski equation at $E_{m}\ll k_{\textrm{B}}T$, which the non-equilibrium properties, like $\mu$ shown in Eq. (4) to determine the heat dissipation, are used to describe the irreversible nature [23, 24, 25]. On the basis of this physical picture, attention have been paid on the low-energy reaction behavior at $E_{m}\sim k_{\textrm{B}}T$, the intermediate stage between the low- and high-temperature limits, by integrating both the reversible and irreversible characteristics into a unique framework [21, 22], which is believed to be the key issue for the theoretical description of the nature of the low-energy reaction.

For the issue of the nucleation and growth of helium-vacancy clusters, most of the multi-scale modeling studies based on KMC focus on the accuracy of input parameters calculated from DFT or MD simulations [4], rather than the more important issue, i.e., to introduce the low-energy reaction mechanism into the KMC scheme[15]. One of the recent endeavors is the work of Ref. [12] (short for “saw-tooth model” in the following), who used Zwanzig [26] and Mori [27] projection operator approach to construct a stochastic model of Brownian motion upon a simplified saw-tooth potential and propose a unique equation to describe the ‘anomalous’ diffusion behavior for single helium in BCC W in the whole temperature region (See Fig. 4 in Ref. [12]). However, the saw-tooth model is derived on the basis of Brownian motion, where the diffusion behavior at the intermediate temperature region is described by the simple interpolation of the behaviors at low- (i.e., Eq. (3)) and high-temperature limits (i.e., Eq. (4)). Further, the thermal fluctuation information at the basin of potential well for helium diffusion is missed due to the usage of simple saw-tooth potential.

To overcome the limitations in saw-tooth model [12], we construct an alternative stochastic model of hopping random walk upon a cosine potential (short for “cosine model” in the following), then derive a modified formula of single helium diffusion in metals, which is suitable for the KMC scheme in the issue of growth of helium-vacancy clusters. As presented in Sec. 2, cosine model will start with an assumption that the stochastic motion of single helium in metals could still be treated act as a step-by-step random walker hopping among its stable sites, i.e., the basins of a cosine potential provided by atoms of the host metal, and the motion of helium around its potential minimum is regarded as a damped harmonic motion, upon which an effective quality factor is introduced to describe the ‘anomalous’ diffusion behavior at the intermediate temperatures. The comparison of prediction of cosine model and the results obtained from molecular dynamics simulations will be then placed in Sec. 3, as well as the discussion. We hope this modified formula could help to develop a more accurate KMC scheme for the long-term evolution of helium-vacancy clusters and bubbles, as well as other low-energy reactions.

2 A modified analytical formula

2.1 Stochastic dynamics model

In BCC metals, a single helium atom occupying the tetrahedral interstitial sites generally exerts a biased periodic crystal potential and results in the local distortion and resonance vibrational modes, as well as the scattering center for phonon modes of the host matrix. With the thermal fluctuations and interaction provided by host atoms, helium atom travels upon the crystal potential until it is trapped by sinks, as schematic as in Fig. 1. Considering the case that helium atom is far away from the sinks in BCC W, the migration energy $E_{m}$ of interstitial helium in BCC W is found as small as $\sim$0.145eV [12], whose diffusion could be regarded as a low-energy reaction. According to saw-tooth model [12], the trajectory of interstitial helium solute inside the metals $X=X(t)$ is governed by a stochastic Langevin equation as

[TABLE]

where $X$ is the instantaneous position of helium atom, with $m^{*}$ the effective mass; $F=-\partial_{X}E(X)$ is the restoring force, with $E(X)$ the periodic crystal potential provided by host matrix, i.e., migratory barrier; $-m^{*}\gamma\dot{X}$ and $f(t)$ represent respectively the actions of dissipation and fluctuation provided by host atoms, with $\gamma$ the phonon-drag friction coefficient and $f(t)$ a Gaussian-type random force, where $\left\langle f(t)\right\rangle=0$, and $\left\langle f(t)f(t^{\prime})\right\rangle=2m^{*}\gamma k_{\textrm{B}}T\delta(t-t^{\prime})$.

The diffusivity $D$ could be derived from the trajectory $X(t)$ governed by Eq. (5), as

[TABLE]

where $\left\langle\cdots\right\rangle$ represents the running-time average. In particular,

[TABLE]

at low- and high-temperature limits, where $D_{0}=g\nu_{0}\lambda^{2}=\nu_{0}\lambda^{2}$ with $g=2/3$ for helium diffusion in BCC metals. The problem is that it is very difficult to obtain the analytical expression of $D$ at the intermediate temperature, i.e., $E_{m}\sim k_{\textrm{B}}T$, for an arbitrary potential $E(X)$ in Eq. (5). Substantial progresses were achieved just in the limiting cases. The representative works are respectively the Lifson-Jackson formula [28] in the large friction limit

[TABLE]

with $\mu=1/m^{*}\gamma$ the classical mobility, and Risken’s expression [29] based on a cosine-type potential in the low friction limit

[TABLE]

In underdamped limit, i.e., $\gamma\gg\nu_{0}$, the defect trajectory shows the long tracks (\gg$$\lambda) inside the metals, otherwise the overdamped condition, i.e., $\gamma\ll\nu_{0}$, leads to a typical short steps (\sim$$\lambda) of hopping between the potential minima [21]. The former case is suitable to describe the high-temperature behavior of helium diffusion [13], but the form of the latter case is more applicable in KMC scheme. Therefore, we would like to derive our analytical expression on the basis of Eq. (8) and introduce an effective quality factor $b=g\cdot\left(\gamma/2\pi\nu_{0}\right)$, to account for the damping feature for helium diffusion upon a cosine potential in the following.

2.2 Derivation of cosine model

Assuming a cosine-type force field of $E(X)$ is exerted by a single interstitial helium in metals as

[TABLE]

The attempt frequency $\nu_{0}$ could be then derived as

[TABLE]

Therefore, the quality factor $b$ is written as

[TABLE]

Meanwhile, the diffusivity could be derived by substituting Eq. (10) into Eq. (8), as [22]

[TABLE]

where $J_{0}(x)$ is the modified Bessel function of the first kind, with $J_{0}(x)$=1 at $x\ll 1$ and $J_{0}(x)$=$\sqrt{(2\pi x)}e^{x}$ at $x\gg 1$. Therefore,

[TABLE]

which satisfies the high-temperature limiting behavior but not the one at low-temperature limit. To solve this problem, we modify Eq. (13) by introducing the quality factor $b$, in form of

[TABLE]

In this regard, we have

[TABLE]

satisfying the low- and high-temperature limiting conditions. The $\gamma^{-1}$ dependence is could be found in Eq. (15), which is consistent with the prediction of a more rigorous theory proposed in Ref. [21]. The occurring rate $\nu$ of single helium hopping event conveniently used in KMC is then derived as

[TABLE]

so that $\nu=\nu_{0}$ at $E_{m}\gg k_{\textrm{B}}T$ and $\nu=\mu k_{\textrm{B}}T/\lambda^{2}$ at $E_{m}\ll k_{\textrm{B}}T$, which satisfies the Arrhenius and Einstein-Smoluchowski equations, respectively.

Eq. (15) is the derived analytical expression of the diffusivity for a single helium solute in metals, which bridges the characteristics of quasi-equilibrium (i.e., the pre-factor $D_{0}$ and migration energy barrier $E_{m}$) and non-equilibrium (i.e., classical mobility $\mu=1/m^{*}\gamma$) in low-energy reaction with an effective quality factor $b$.

3 Numerical example of helium in BCC W

Fig. 2 shows the migratory energy profile of single interstitial helium solute in BCC W hopping between two neighboring tetrahedral sites, which is calculated using modified conjugated gradient (MCG) [30] method in Ref. [12]. Accordingly, a cosine-type potential, i.e., Eq. (10), is used to describe the detail information with $E_{m}=0.145$eV. It can be seen the cosine potential is in general agreement with the migratory profile, including the curvature at the basin of potential well and the intermediate part. In addition, the saddle-point state, i.e., the top of migratory profile, is a meta-stable state but with a tiny barrier, $<0.01$eV, which could not be reproduced using Eq. (10).

Note that, parameters in Eq. (15) are all of many-body nature, which can be obtained from either atomistic modeling or experimental measurements. Using the parameters estimated in Ref. [12] for single helium diffusion in BCC W, i.e., $D_{0}=4.3$ Å/ps, $E_{m}=0.12$eV, and $\mu=8.4\textrm{\AA}^{2}/\textrm{ps}/\textrm{eV}$, getting $b=1.36$, the temperature dependence of diffusivity in Eq. (15) could be achieved. Here, this input migration energy $E_{m}=0.12$eV, obtained by Arrhenius-fitting the dynamic simulation results at low-temperature limit, is a slightly smaller than that obtained in molecular static calculation (i.e., $E_{m}=0.145$eV from MCG method), which might be arising from the entropic term, or the temperature dependence of migration energy. The results are plotted in Fig. 3, as well as the comparison with MD simulation results [12, 13] and the prediction using saw-tooth model [12] and the Arrhenius equation (Eq. (3)) with the same input parameters. It can be seen that our current model, i.e., cosine model in Fig. (3), could reproduce the feature of low-energy reaction at various temperatures ranging from $\sim$0K to 3000K. A slight difference in results between cosine model Eq. (15) and saw-tooth model appears in the low-temperature region, i.e., 300K $<T<500$K, where the Arrhenius behavior is revealed, because of the usage of different potential to describe the characteristics of the equilibrium state for helium migration. On the other hand, the Arrhenius equation give rise to a relative large overestimation for diffusivity at high-temperature region, e.g., $T>600$K.

In a word, in the case of $E_{m}\ll k_{\textrm{B}}T$, Einstein equation is more appropriate to describe the diffusion behavior of helium in metals, due to its low migration energy barrier. A relative error $\varepsilon(T)$ could be defined as

[TABLE]

to check the error of the application of Arrhenius equation in the high-temperature region, e.g., $T>600$K in BCC W. Here, $D_{A}(T)$ and $D_{E}(T)$ are the diffusivity predicted using Arrhenius equation (i.e., Eq. (3)) and Einstein equation (i.e., Eq. (4)), respectively. Fig. (4) presents the temperature dependence of $\varepsilon$. It can be seen $\varepsilon$ could be as large as 50% at 1000K$<T<2000$K, resulting in a significant error of growth rate of helium-vacancy clusters or bubbles, thus the kinetics of microstructural evolution in nuclear materials.

4 Conclusion

For interstitial helium solute in metals, its non-Arrhenius diffusion behavior at finite temperatures could be completely described by the quasi-equilibrium diffusion parameters, e.g., pre-factor $D_{0}$ and migration energy $E_{m}$, in the framework of Arrhenius equation. The non-equilibrium diffusion parameter, e.g., classical mobility $\mu$, should be involved to describe the diffusion behavior at high-temperatures. In this paper, we proposed a stochastic model and derived a modified formula, with those three parameters involved, to describe the non-Arrhenius diffusion behavior of single interstitial helium in metals, which could help to develop a more accurate lattice based kinetic Monte Carlo scheme for the long-term evolution of helium-vacancy clusters and bubble. Our model is built up on the basis of a stochastic model of Brownian motion upon a cosine potential. We assumed that the helium migration could be treated as a step-by-step random hopping event so that the diffusivity satisfy the Lifson-Jackson’s formula [28] under overdamped condition. Further, by introducing an effective quality factor to account for the damping feature, a modified formula for diffusivity (i.e., Eq. (15)) or jump-frequency (i.e., Eq. (17)) of single helium solute is achieved. Of course, the story is not ending, because the related diffusion parameters still have to be obtained empirically or from other calculations. Compared to the saw-tooth model [12], the cosine potential used here is more comfortable to estimate the migratory information near the basin of potential well self-consistently, such as the attempt frequency, which could be presented in our further paper. It has to be noted that, our modified formula is not restricted in the case of helium diffusion, but could be promoted to the general cases of low-energy reaction in materials science.

Acknowledgment

This work is initiated and funded by the Guangzhou Science and Technology Project (No. 201707020002) and NSFC (No. 11672339, No. 11602310, No 11602311), to which the authors are thankful. Y. Zheng also thanks support from the Special Program for Applied Research on Super Computation of the NSFC-Guangdong Joint Fund (the second phase), Fok Ying Tung Foundation, Guangdong Natural Science Funds for Distinguished Young Scholar and China Scholarship Council.

Reference