Low temperature diffusivity of self-interstitial defects in tungsten

To accurately determine the migration barrier for a crowdion defect in tungsten, we performed ab initio calculations using the VASP package [25, 26], with the AM05 [27] pseudopotential including 12 valence electrons. The simulation supercell contained 129 atoms, with 5 × 5 × 5 gamma-point centered k-point grid and a plane-wave energy cutoff of 450 eV. The ion positions and simulation box shape of the 111 crowdion and 111 dumbbell structures were fully relaxed with energy and force convergence at ${10}^{-7}$ eV and ${10}^{-4}$ eV ${\mathring{\rm A} }^{-1}$. To determine the energy barrier the nudged elastic band method [28] calculations were performed using seven images between the ground state and saddle point. After completing the calculations, elastic corrections, as proposed and implemented by Varvenne et al [29], were applied to eliminate the effect of elastic periodic image interactions on the system energy. The results of these calculations are shown in figure 1(b), along with similar calculations [30] performed in a supercell of 2001 atoms using an embedded atom method (EAM) interatomic potential [31] developed by Marinica et al [18], which we also use in our molecular dynamics simulations. The EAM calculations, in good agreement with other available interatomic potentials [7], predict a very small crowdion migration barrier of 17.2 meV. This is consistent with classical molecular dynamics simulations, performed in this work and elsewhere [7, 8], which find, in agreement with experiment [14], that crowdions execute free Brownian motion above 150 K, as the migration barrier is overwhelmed by thermal agitations. However, as can be seen in figure 1(b), our DFT calculations yield an even smaller migration barrier of 2.07 meV, in agreement with previous ab initio based approaches to model this defect [32]. As we will see, this extremely small magnitude of the migration barrier has profound implications for the nature of the low temperature defect motion.

As mentioned above, below a crossover temperature ${{T}}_{\mathrm{cross}}$, non-classical phonon assisted tunnelling trajectories emerge as stationary action solutions in an imaginary time path integral evaluation of the defect density matrix [19]. ${{T}}_{{\rm{cross}}}$ is given by ${{T}}_{{\rm{cross}}}\simeq {\rm{\hslash }}| {\omega }_{0}^{s}| /2\pi {{k}}_{{\rm{B}}}$, where ${\omega }_{0}^{s}$ is the unstable frequency at the top of the migration barrier. Using ${\omega }_{0}^{s}$ extracted from our EAM calculations for a crowdion in tungsten gives ${{T}}_{{\rm{cross}}}=1.8$ K, whilst using the migration barrier found in DFT calculations gives ${{T}}_{{\rm{cross}}}=0.6$ K. This very low crossover temperature allows us to neglect all tunnelings phenomena above ∼1 K and focus only on the effect of quantum zero-point fluctuations on the thermally activated defect motion. In the supplementary material (available online at stacks.iop.org/NJP/19/073024/mmedia) we also provide an estimate of the decoherence temperature, finding ${{T}}_{\mathrm{decoh}}\simeq {10}^{-8}$ K, as expected due to tungsten's large atomic mass. As existing empirical potentials cannot capture the extremely low crowdion migration barrier found in our DFT simulations, the observed defect diffusivities found in our molecular dynamics simulations must be modified to account for the discrepancy. We return to this point below, when developing our theoretical treatment of the low temperature diffusion.

As tunnelling phenomena can be neglected for system temperatures above $\simeq 1K$ for the SIA clusters considered here, semi-classical simulations, which only retain quantum zero-point motion, are adequate to model the low temperature diffusion. To this end, we have implemented the colored noise quantum heat bath proposed by Barrat and Rodney [24] (see supplementary material) in the MD and spin-lattice dynamics program SPILADY [33]. The quantum heat bath technique is an approximate method to capture the quantum lattice vibrations of solids, which cannot treat systems with pathologically strong anharmonicity [34–36] as zero-point energy is erroneously distributed between all the vibrational modes. However, when parametrised with care, the method is in principle able to treat realistic potentials [24, 37], providing the heat bath friction constant is sufficiently high to prevent zero-point energy leakage but sufficiently weak to allow correct anharmonic behaviour. Whilst this is not in general possible for all systems [37], the strength of atomic bonding in tungsten allows this to be achieved with a parametrisation similar to that used in previous studies [24] (see supplementary material).

The simulation supercell contained 16001 atoms, with supercell vectors of $20a[001],20a[010]$ and $20a[001]$. The crowdion configuration was structurally relaxed with a conjugate gradients minimization and then thermalized with either a quantum or classical heat bath for 0.1 ns. Simulation trajectories were obtained at ten temperatures between 10 and 200 K, at intervals of 10 K between 10 and 80 K, with longer trajectories broken into four independently thermalized segments. The quantum and classical heat bath trajectories were all 4 ns in length, yielding robust statistical data on the defect diffusivity. The defect position was extracted by first time-averaging over a period of ∼0.05 ps, applying a high-pass energy filter to reveal the defect core, then taking the center of mass of the remaining atoms. At lower temperatures, the defect trajectory is characterized by isolated stochastic jumps; in this case, the diffusivity is given by $D={b}^{2}{N}_{\mathrm{jump}}/2{t}_{\mathrm{sim}}$, where ${N}_{\mathrm{jump}}$ is the number of defect jumps, ${t}_{\mathrm{sim}}$ the total simulation time and b the jump length between potential minima. At higher temperatures, typically above 150 K for the simulations presented here, the defect trajectory exhibits a continuous stochastic motion and the diffusivity $D=\langle {(\delta x)}^{2}\rangle /2\delta t$ is calculated through analysis of the mean squared displacement [8, 11].

The simulated defect diffusivities are displayed on Arrhenius plots $(1/{{k}}_{{\rm{B}}}{T},\mathrm{ln}D)$ in figure 2. The diffusivities extracted from classical heat bath molecular dynamics clearly exhibits the classical Arrhenius relationship $\mathrm{ln}| {D}_{\mathrm{cl}}({T})| =\mathrm{ln}| {D}_{0}| -{\rm{\Delta }}V/{{k}}_{{\rm{B}}}{T}$ below 100 K, with an activation energy ${\rm{\Delta }}V$ of 0.018 eV, in very good agreement with the value 0.0172 eV extracted from static calculations using the same EAM potential. From previous studies of crowdion motion in bcc metals, it is known that the Arrhenius behavior breaks down at temperatures above 150 K; at high temperature the diffusivity is that of a freely diffusing particle, $D={{k}}_{{\rm{B}}}{T}/\gamma$, where γ is the frictional coupling between the defect and thermal vibrations [7, 8, 16]. However, in contrast to purely classical simulations, the diffusivities extracted from quantum heat bath molecular dynamics deviate significantly from the classical law at low temperature, indicating the influence of quantum zero-point motion on defect migration.

Zoom In Zoom Out Reset image size

We have seen that due to the very low migration barriers found in static calculations, we are able to neglect all tunneling phenomena for system temperatures above 1 K, allowing the use of semiclassical theories which retain only quantum expectation values in the vibrational spectrum [23]. In previous work, a theoretical and numerical scheme was developed to analyze the thermal vibrations surrounding migrating crystal defects, initially applied to study the phonon drag force [8, 16]. The method uses a defect migration pathway as calculated from a barrier climbing calculations, defining the defect position coordinate r to be the reaction coordinate along this pathway. By considering possible structural variations along the migration pathway ($\delta r\in {\mathbb{R}}$) and perpendicular to the migration pathway (${\bf{x}}\in {{\mathbb{R}}}^{3{\rm{N}}-1}$) one can show that, as defect motion is not in general an eigenmode of the host crystal [16, 38], the potential energy V can be expanded as

$$\begin{eqnarray}&&V(r+\delta r,{\bf{x}})={V}_{0}(r+\delta r)+\sum {c}_{k}{x}_{k}\delta r+\displaystyle \sum _{k}m{\omega }_{k}^{2}\displaystyle \frac{{x}_{k}^{2}}{2},\end{eqnarray} \tag{ 1 }$$

where ${V}_{0}(r)$ is the periodic migration potential shown in figure 1(b). The vibrational frequencies $\{{\omega }_{k}\}$ and coupling constants $\{{c}_{k}\}$ are related to the Hessian matrix of second derivatives [8, 16]. The coupling constants $\{{c}_{k}\}$ allow the defect system to be thermalized through interaction with the surrounding thermal vibrations, giving rise (in a classical picture) to a dissipative drag force $-\gamma \dot{r}$ and a fluctuating thermal force of square magnitude $\gamma {{\rm{k}}}_{{\rm{B}}}{T}$, where the frictional coupling γ can be expressed [8, 16] as a function of the $\{{\omega }_{k}\}$ and $\{{c}_{k}\}$. We have evaluated the $\{{c}_{k}\}$ for various defects in bcc metals [8, 16], including the crowdion defect considered here; explicit expressions for $\{{\omega }_{k}\}$, $\{{c}_{k}\}$ and γ are given in the supplementary material.

To derive an effective migration rate in the temperature regime of interest, we note that equation (1) is the widely studied 'oscillator bath' Hamiltonian used in the theory of OQS [23], which aims to integrate away the vibrational degrees of freedom ${\bf{x}}$ to leave an effective distribution function for the defect density matrix, which then yields a Fokker–Planck equation for the Wigner quasi-probability function $W(\delta r,{p}_{r})$ [19]. Expanding the potential around r = 0 to obtain $V{(\delta r)=V(0)+m{\omega }_{0}^{2}(\delta r)}^{2}/2$, a standard result is that $W(\delta r,{p}_{r})$ follows the approximate evolution equation [23]

$$\begin{eqnarray}&&\displaystyle \frac{\partial W}{\partial t}=\left(\displaystyle \frac{\partial }{\partial {p}_{r}}(m{\omega }_{0}^{2}\delta r+\gamma {p}_{r})-\displaystyle \frac{{p}_{r}}{m}\displaystyle \frac{\partial }{\partial r}+D\displaystyle \frac{{\partial }^{2}}{\partial {p}_{r}^{2}}\right)W,\end{eqnarray} \tag{ 2 }$$

where explicit expressions for D and $\gamma$ are given in the supplementary material. The evolution equation (2) has a steady state distribution ${W}_{\infty }(\delta r,{p}_{r})={Z}^{-1}\exp (-[V(\delta r)+{p}_{r}^{2}/2m]\gamma /D)$, that clearly resembles a classical Gibbs distribution with an effective temperature (see supplementary material)

$$\begin{eqnarray}&&{{T}}_{\mathrm{OQS}}\equiv \displaystyle \frac{D}{{{k}}_{{\rm{B}}}\gamma }=\displaystyle \frac{{\rm{\hslash }}{\displaystyle \sum }_{{k}}{c}_{{k}}^{2}{\omega }_{{k}}^{-2}\sin ({\omega }_{{k}}{\tau }_{C})\coth \left(\tfrac{{\rm{\hslash }}{\omega }_{{k}}}{2{{k}}_{{\rm{B}}}{T}}\right)}{2{{k}}_{{\rm{B}}}{\displaystyle \sum }_{{k}}{c}_{{k}}^{2}{\omega }_{{k}}^{-3}\sin ({\omega }_{{k}}{\tau }_{C})},\end{eqnarray} \tag{ 3 }$$

where ${\tau }_{C}=0.05$ ps is the calculated correlation time of the thermal force acting on the defect [8]. It is straightforward to show that the OQS temperature ${{T}}_{\mathrm{OQS}}\to {T}$ at high temperature, but at low temperature ${{T}}_{\mathrm{OQS}}$ tends to a constant value due to zero-point fluctuations. In the inset of figure 2 we compare ${{T}}_{\mathrm{OQS}}$ to the effective 'classical' temperature of a solid state system with normal modes $\{{\omega }_{k}\}$ [21]

$$\begin{eqnarray}&&{{T}}_{{\rm{C}}}=\displaystyle \frac{1}{3({N}-1){{k}}_{{\rm{B}}}}\displaystyle \sum _{{k}=1}^{3{N}-3}\displaystyle \frac{{\rm{\hslash }}{\omega }_{{k}}}{2}\coth \left(\displaystyle \frac{{\rm{\hslash }}{\omega }_{{k}}}{2{{k}}_{{\rm{B}}}{T}}\right),\end{eqnarray} \tag{ 4 }$$

which can also be evaluated numerically through Maxwellian analysis of the atomic velocities to give ${{T}}_{{\rm{C}}}=m\langle {V}^{2}\rangle /(3({N}-1){{k}}_{{\rm{B}}})$. We see that both effective temperatures show a very similar behavior, with ${{T}}_{\mathrm{OQS}}$ tending to a slightly lower temperature than ${{T}}_{{\rm{C}}}$. Equations (2) and (3) show that, in the approximation of a harmonic potential well, the coupling of a crystal defect to the surrounding vibrational modes results in an effectively classical Fokker–Planck equation (2), with all quantum mechanical effects contained in the ${\rm{\hslash }}$ dependent effective system temperature ${{T}}_{\mathrm{OQS}}$. Although no exact result exists for the general nonlinear migration potential in the open quantum systems theory setting, an approximate diffusivity D_OQS that accounts for quantum zero-point motion can be produced by using the effective temperature ${{T}}_{{\rm{OQS}}}$ in the classical expression for one dimensional motion in a periodic potential [39, 40]

$$\begin{eqnarray}&&{D}_{\mathrm{OQS}}=\displaystyle \frac{{{k}}_{{\rm{B}}}{{T}}_{\mathrm{OQS}}{b}^{2}/\gamma }{{\int }_{0}^{b}\exp \left(-\tfrac{F(r)}{{{k}}_{{\rm{B}}}{{T}}_{\mathrm{OQS}}}\right){\rm{d}}r{\int }_{0}^{b}\exp \left(\tfrac{F(r)}{{{k}}_{{\rm{B}}}{{T}}_{\mathrm{OQS}}}\right){\rm{d}}r},\end{eqnarray} \tag{ 5 }$$

where $F(r)$ is the free energy of migration, given in the harmonic approximation by

$$\begin{eqnarray}&&F(r)={V}_{0}(r)+{{k}}_{{\rm{B}}}{{T}}_{\mathrm{OQS}}\displaystyle \sum _{{k}}\mathrm{ln}| {\omega }_{{k}}| +{F}_{0},\end{eqnarray} \tag{ 6 }$$

and where F₀ is a constant term that does not affect our results. When ${{{k}}_{{\rm{B}}}{T}}_{\mathrm{OQS}}$ is much smaller than the free energy barrier ${\rm{\Delta }}F$, ${D}_{\mathrm{OQS}}$ is given by the well known expression from the Kramers transition state theory [39]

$$\begin{eqnarray}&&\mathop{\mathrm{lim}}\limits_{{{k}}_{{\rm{B}}}{{T}}_{\mathrm{OQS}}\ll {\rm{\Delta }}F}{D}_{\mathrm{OQS}}=\displaystyle \frac{m{\omega }_{0}^{g}| {\omega }_{0}^{s}| }{\gamma }\left(\displaystyle \prod _{{k}=1}^{3N-4}\displaystyle \frac{{\omega }_{{k}}^{g}}{{\omega }_{{k}}^{s}}\right)\times \exp \left(-\displaystyle \frac{{\rm{\Delta }}V}{{{k}}_{{\rm{B}}}{{T}}_{\mathrm{OQS}}}\right),\end{eqnarray} \tag{ 7 }$$

where ${\omega }_{0}^{g}$ and $| {\omega }_{0}^{s}|$ are the frequencies at the ground and saddle point of the static migration potential $V(x)$, taken here either from EAM or DFT calculations, whilst the ${\{{\omega }_{k}^{o,s}\}}_{1}^{3N-4}$ are the eigenfrequencies at the ground state and saddle state [39], which here we take only from EAM calculations. At effective temperatures such that ${{{k}}_{{\rm{B}}}{T}}_{\mathrm{OQS}}\gg {\rm{\Delta }}F$ we obtain the free diffusion coefficient

$$\begin{eqnarray}&&\mathop{\mathrm{lim}}\limits_{{{{k}}_{{\rm{B}}}{T}}_{\mathrm{OQS}}\gg {\rm{\Delta }}F}{D}_{\mathrm{OQS}}=\displaystyle \frac{{{k}}_{{\rm{B}}}{{T}}_{\mathrm{OQS}}}{\gamma }.\end{eqnarray} \tag{ 8 }$$

Equation (5) constitutes our main theoretical result, an expression for the migration rate of a crowdion in a heavy atom system, valid down to temperatures of at least 1 K. In figure 2 equation (5) is illustrated using the migration potential ${V}_{0}(x)$ found in molecular statics and DFT calculations. We see that using the migration barrier of 17.2 meV derived from the empirical potential, ${D}_{\mathrm{OQS}}$ is in good agreement with the diffusivity extracted from our quantum heat bath molecular dynamics simulations. When using the much lower migration barrier of 2.07 meV found in our DFT calculations, as such a small migration barrier is easily overcome by all effective classical temperatures, the diffusivity is essentially identical at all system temperatures to the free diffusion limiting form (8).

Whilst we have developed an accurate model for single SIA crowdions, irradiated materials will in practice have a large distribution of SIA cluster sizes [5]. Migration barriers of larger SIA loops have not yet been calculated in DFT due to computational restrictions on the system size; however, from previous classical heat bath molecular dynamics studies of SIA cluster motion in bcc metals, it is known that for larger SIA clusters the migration barrier is approximately constant [11, 15], with a prefactor that scales as the inverse square root of the cluster size N [15]. In addition, above around 200 K the cluster exhibits free Brownian motion [16], where $D={{k}}_{{\rm{B}}}{T}\,{\gamma }^{-1}$, with a friction coefficient γ that also scales as $1/\sqrt{N}$. We have performed classical molecular dynamics simulations of a 7-SIA cluster, using trajectories of up to 94 ns at low temperature to obtain statistically robust measurements of the diffusion. We find a migration barrier of around 70 meV at low temperature and a friction coefficient of 13 meV ${{\rm{\mathring{\rm A} }}}^{-2}$ ps⁻¹ at high temperature. However, as our DFT calculations revealed that the crowdion migration barrier was 12% of the EAM value; accounting for this discrepancy yields an SIA cluster migration barrier of 8 meV. The effects of such a small migration barrier will be overcome at effective classical temperatures of $50\,{\rm{K}}$, allowing us to use the high temperature approximation (8) for the loop mobility. Collating the information from ab initio calculations and quantum heat bath molecular dynamics simulations thus gives our main result, the diffusivity of an N-SIA cluster in tungsten

$$\begin{eqnarray}&&{D}_{\mathrm{OQS}}({T},N)=\displaystyle \frac{{{k}}_{{\rm{B}}}{{T}}_{\mathrm{OQS}}({T})}{2\sqrt{\pi N}{\gamma }_{0}},\end{eqnarray} \tag{ 9 }$$

where ${\gamma }_{0}=1.38$ meV ${{\rm{\mathring{\rm A} }}}^{-2}$ ps⁻¹ and the effective temperature ${{T}}_{\mathrm{OQS}}$ is formally given by equation (3). A simple form of (9), appropriate for implementation in numerical simulations, can be obtained after two simplifying approximations. Firstly, as ${\gamma }_{0}$ is proportional to the magnitude of the Hessian matrix of second derivatives (see supplementary material and [8]) we conclude that ${\gamma }_{0}$ is proportional to μ, the shear modulus. In addition, we have found that ${{T}}_{\mathrm{OQS}}$ can be well approximated by the simple function $\sqrt{{{{T}}_{{\rm{D}}}}^{2}/{a}^{2}+{{T}}^{2}}$, where a is a fitting parameter. Under these approximations we obtain the readily usable form

$$\begin{eqnarray}&&{D}_{\mathrm{OQS}}(N)\simeq \displaystyle \frac{7.04\times {10}^{4}}{\mu \sqrt{N}}\sqrt{{\left(\displaystyle \frac{{{T}}_{{\rm{D}}}}{3.6}\right)}^{2}+{{T}}^{2}}\quad \mu {{\rm{m}}}^{2}\,{{\rm{s}}}^{-1},\end{eqnarray} \tag{ 10 }$$

where $\mu =400$ GPa and ${{T}}_{{\rm{D}}}=310\,{\rm{K}}$ are the shear modulus and Debye temperature of tungsten, giving a numerical value of ${D}_{\mathrm{OQS}}(N)\simeq 176\sqrt{{85}^{2}+{{T}}^{2}}/\sqrt{N}\mu {{\rm{m}}}^{2}\,{{\rm{s}}}^{-1}$. We note that the approximation in equation (10) predicts a breakdown in classical behaviour once $T\le {T}_{{\rm{D}}}/3.6$, in agreement with a high temperature expansion of the Bose–Einstein distribution [41].

Our main result, equation (9), provides a predictive calculation of the diffusivity of self-interstitial defects in pure tungsten, valid down to cryogenic temperature, which are presently used in a wide variety of experiments investigating post irradiation microstructure evolution. Due to the influence of zero-point motion, we find that the diffusivity remains very high, of order ${10}^{3}\,\mu {{\rm{m}}}^{2}\,{{\rm{s}}}^{-1}$ for a 100-SIA cluster. This high value, which has never been measured or estimated previously, has clear implications for transmission microscopy experiments conducted in $1\mbox{--}3\,\mu {\rm{m}}$ in thick foils at low temperatures [4, 6]. Using a diffusivity of ${10}^{3}\,\mu {{\rm{m}}}^{2}\,{{\rm{s}}}^{-1}$ we expect clusters to diffuse across typical foil thicknesses in around 0.03 s, leading to a likely depreciation in the cluster population even at cryogenic temperatures in the absence of impurities or other imperfections. Furthermore, as the diffusivity remains constant below 50 K, the effective activation energy is zero and therefore it is unlikely that self-interstitial defect migration will feature in the measured spectrum of energies derived from resistivity recovery experiments.

We gratefully acknowledge many stimulating discussions with K. Arakawa and M.-C. Marinica concerning the nature of low temperature defect diffusion. This work was stimulated by the preliminary experimental electron microscope observations of low temperature migration of defects reported in [6] and elsewhere. While our work involves a number of approximations, particularly about the absence of impurities or other defects, which in any experimental observation likely have a profound effect on the measured transport properties, we believe that the general approach described above will help with the development of a comprehensive theoretical framework within which experimental data can be interpreted. Density functional calculations were carried out using the EUROfusion Marconi supercomputer facility in Italy. This work has been carried out within the framework of the EUROfusion Consortium and has received funding from the Euratom research and training programme 2014–2018 under grant agreement No 633053 and from the RCUK Energy Programme [grant number EP/P012450/1]. To obtain further information on the data and models underlying this paper please contact PublicationsManager@ukaea.uk. The views and opinions expressed herein do not necessarily reflect those of the European Commission.