Strain Effects on the Diffusion Properties of Near-Surface Self-Interstitial Atoms and Adatoms in Tungsten

Tungsten (W) is a candidate for the plasma-facing components and divertor in future fusion applications. The material will be subject to a large particle influx (mainly helium and hydrogenic species) that will form bubbles. As bubbles grow, they compress the material, adding to thermal stresses, and eject self-interstitial atoms (SIAs—isolated or in clusters) to release internal pressure. These SIAs diffuse towards the surface in large stress/strain fields and on the surface are thought to act as precursors for nanotendril formation (also known as fuzz) that develops on the material surface modifying its morphology. In this work we analyze the effect of strain on the diffusion properties of both SIAs and adatoms. Relying on atomistic simulations, we compute the average time that a SIA created in the center of a tungsten slab takes to reach a (110) surface for different strains and temperatures. This time relates to the SIA diffusivity and allows to compute the activation energy and dipole tensor including surface effects. We observe a large strain effect that significantly modifies the propensity for SIAs to reach the surface and, hence, to cluster to form dislocation loops in the bulk crystal. Strain also alters the diffusivity of the adatom although to a lesser extent. Finally, we report on the resulting surface roughness evolution and its dependence on strain.

1 Introduction

In fusion environments, plasma-facing components (PFCs) are exposed to stringent conditions as imposed by the plasma, required for the nuclear reaction to occur. Such PFCs must withstand extreme conditions of temperature and particle fluxes that modify the material microstructure, potentially leading to failure. (Knaster et al., 2016) Understanding the materials response upon such extreme environments is therefore crucial for the development of clean, safe and abundant fusion energy generation systems. Tungsten (W) is the main candidate as PFC due to its favorable properties, such as low sputtering yield, high sputtering threshold, high thermal conductivity and high melting point. (Suzuki et al., 2003; Ihli et al., 2006).

Despite these advantages, certain drawbacks have also been observed, such as the formation of a fuzz-like nanostructure in the W surface due to the helium (He) ash implantation, which leads to the formation of high-pressure bubbles that eject self-interstitial atoms (SIAs) to relax such high stresses, consequently generating vacancies in the bubbles. The diffusion of these SIAs, either individually or in clusters, towards the surface in the presence of stresses generated thermally and by the bubble distribution, seems to be an important mechanism contributing to the observed surface morphological evolution and nanotendril formation (Takamura et al., 2006; Baldwin et al., 2009; Yang et al., 2015; Wang et al., 2017). This fuzz nanostructure modifies the erosion propensity and potentially causes embrittlement of the divertor component, (Wirth et al., 2015) leading to particle sputtering that may disrupt the plasma. Fuzz forms under certain plasma conditions, when W is exposed to He and He-deuterium (D) mixed plasmas but not for pure D plasmas (Baldwin et al., 2009), at exposure temperatures in the range of 900–2000 K and incident ion energies around 20–100 eV with incoming fluxes on the order of 10²⁴ m⁻² s⁻¹. (Kajita et al., 2009) Note that these incident energies are insufficient to generate radiation damage in the form of vacancies and self-interstitials, but lead to the formation of over-pressurized bubbles and subsequent trap mutation.

Many different models have been proposed in the literature to explain fuzz formation. (Krasheninnikov, 2011; Martynenko and Nagel’, 2012; Lasa et al., 2014; Fiflis et al., 2015; Ito et al., 2015; Ito et al., 2018; Dasgupta et al., 2019) Recently, Dasgupta et al. (Dasgupta et al., 2019; Dasgupta et al., 2020; Chen et al., 2021) have developed a continuum-scale model that tries to understand and predict the onset of fuzz formation depending on the exposure conditions. In this model, internal stresses due to the formation of highly over-pressurized bubbles is the main driving force for triggering surface instabilities leading to the formation of fuzz. Each bubble is modeled as an Eshelby inclusion acting as a source of hydrostatic compressive stress in the matrix. The model assumes a negligible contribution of particle sputtering to the surface morphology as the incident energies of the implanted He atoms are low. The model describes the change of a surface height function, h = h(x,y,t), in terms of transport of SIAs in one dimension toward the surface, as well as two-dimensional curvature-driven and stress-driven surface diffusion. In that model, the dependence of the process on strain/stress is included in the chemical potential of surface atoms. However, the fluxes also depend on the diffusion coefficient of the SIAs, both in the bulk and on the surface (after they reach the surface and become adatoms), and, in general, both of these diffusion coefficients are functions of the local strain fields. This dependence is still to be incorporated into larger-scale models. The purpose of this study is to better understand the effect of strain on the transport coefficients of SIAs in W in order to improve the predictive capabilities of the continuum-scale surface evolution model.

In this work we rely on atomistic simulations to systematically study the effect of strain on the diffusivities of SIAs and surface adatoms in W. We first describe the methodology that has been followed in this work in Section 2. We present the results in Section 3 where we show the significant effect of strain in the transport of SIAs in bulk W, with a less significant contribution to surface adatom diffusion. In Section 4 we discuss the implications of these results in the study of fuzz formation. Lastly, we draw the main conclusions of our study in Section 5.

2 Methodology

We have used molecular dynamics (MD) simulations as implemented in the LAMMPS code (Plimpton, 1995) with an Embedded Atom Method (EAM) interatomic potential developed by Marinica et al. (Marinica et al., 2013) to model the behavior of body-centered cubic (bcc) W; this interatomic potential was designed specifically targeting defect properties in bcc W. We have used three different sets of simulations, with samples modeled by slab supercells oriented in the $\left[11\overline{2}\right]$, $\left[111\right]$, and $\left[1\overline{1}0\right]$ directions, with free surfaces with z normals and a time step of 2 × 10⁻³ ps.

The first set employs a sample of dimensions 61.5 × 76.1 × 86.6 Å³, with two free surfaces in the z direction and $\mathrm{26,880}$ atoms. This orientation has been chosen since the $\left(1\overline{1}0\right)$ surface is the lowest energy surface in W. (Marinica et al., 2013) A SIA is added at the middle of the sample and the time to reach the surface is computed as a function of temperature and equibiaxial strain. We have used three different temperatures, 700, 850, and 1,000 K, enforced using a Langevin thermostat with relaxation time equal to 1 ps, and strains ranging from −0.007 to 0.007. The biaxial strain in the x and y directions was applied after steady-state thermal expansion was reached at each temperature running at 0 pressure for 200 ps. Twenty independent simulations were run for each case.

The second set of simulations intended to compute the effect of strain on the diffusivity of adatoms on a W$\left(1\overline{1}0\right)$ surface. In this case we run simulations at five different temperatures, from 1,000 K to 1,400 K, five different equibiaxial strains, in the range from −0.008 and 0.008, and carry out five independent runs for each set of conditions. Again, the strain was applied after steady-state thermal expansion was reached at each temperature. These dynamical simulations were complemented performing nudged-elastic band (NEB) (Henkelman et al., 2000; Henkelman and Jónsson, 2000) computations to obtain the minimum energy path (MEP) between two neighboring equilibrium adatom sites at the various levels of applied equibiaxial strain. These NEB simulations were performed in the samples with the volume obtained at 1,000 K for the different strains.

The last simulation set aimed at analyzing the effect of biaxial strain on the morphological evolution of a W$\left(1\overline{1}0\right)$ surface. Five strains, from −0.008 to 0.008 were applied following the same methodology as above at a temperature of 1,000 K and one SIA was inserted randomly into the system every 200 ps, mimicking the ejection of SIAs from He bubbles. The initial sample size was 123.0 × 152.2 × 86.6 Å³ with 107,520 atoms. These conditions result in a flux of SIAs of 2.67 × 10²⁵ m⁻²s⁻¹, overestimating experimental He implantation fluxes that range from 10²¹ to 10²³ m⁻² s⁻¹. Such overestimation highlights one of the main limitations of MD models in terms of the total time that can be simulated. Nevertheless, these simulations provide valuable information about the basic mechanisms occurring in the W slab as SIAs are added and diffuse towards and on the surface and the strain effect on such processes.

To track the position of the SIA, we performed a common neighbor analysis (CNA) (Honeycutt and Andersen, 1987) to obtain the atoms that are not located in bcc lattice positions. For the first passage time simulations, a SIA is specified to reach the surface when all the atoms in the bulk crystal are arranged in bcc lattice sites.

3 Results

3.1 Self-Interstitial Atom Diffusion in Bulk Tungsten

According to the first passage time model, (Redner, 2001) the most probable time for a diffusing particle to reach a point in space from a given origin in 1D is given by

${\tau }_n=\frac{\text{Δ}{x}^2}{6D_{SIA}}(1)$

where Δx is the distance traveled by the walker (in this case a SIA) and D_SIA is the SIA diffusivity. We have estimated this time running the first set of simulations described in Section 2. We have employed these results to fit a diffusion coefficient following transition state theory (TST) in the form

$D_{SIA}=D_0\text{exp}\left\{\frac{-\text{Δ}H}{k_{B}T}\right\}=D_0\text{exp}\left\{\frac{-(\text{Δ}E+{\sigma }_a:\epsilon \text{Ω})}{k_{B}T}\right\}(2)$

with D₀ a diffusivity prefactor, ΔH the activation enthalpy, ΔE the activation energy, σ_a an activation stress tensor related to the difference of the stress tensor at the saddle point and at the initial configuration, ε the applied strain tensor and Ω the atomic volume as a reference constant volume for the system, which for W is Ω = 15.8 Å³. It should be mentioned that σ_a and ε in Eq. 2 denote the respective stress and strain tensors and the corresponding product is a double dot product (carried out through the contraction operator), which, multiplied by the atomic volume, gives a scalar quantity with units of energy. Fitting simultaneously the values gathered from MD we obtain D₀ = 12.9 Å² ps⁻¹, ΔE = 0.19 eV, and σ_a = 110.6 GPa, considering tensile strains as positive. The large value of the activation stress highlights the strong effect of the biaxial strain on the SIA diffusivity in the bulk crystal. Furthermore, the positive sign of σ_a implies that the compressive strain leads to a decrease in the activation enthalpy, increasing the diffusivity, and hence, decreasing the first passage time. Note that ΔE is an effective activation energy in the $\left[1\overline{1}0\right]$ direction and that it contains the effect of the free surface.

To better understand the role of the free surface, we have computed the diffusivity of a SIA in bulk W at ε = 0 and 700 K, 850 K and 1,000 K, respectively. We computed the mean-squared displacement tracking the center of mass of the non-bcc atoms obtained from CNA and used an Einstein relation to calculate the diffusivity. The effective activation energy following from this analysis resulted in a value of around 0.06 eV, which is significantly lower than the value obtained in the presence of the surface (0.19 eV). This value is fully consistent with recent experiments (Heikinheimo et al., 2019) that provide an upper bound for the SIA migration energy of 0.1 eV. To rationalize this result, we have computed the strain induced by the free surface relaxation, which generates a strain in the sample ε_s = 0.006 in the direction normal to the free surface at 1,000 K. Using σ_a = 110.6 GPa calculated from the first passage time and the same Ω as above, the increase in activation energy is ≈0.066 eV, resulting in an activation enthalpy of 0.126 eV, much closer to the value obtained from the first passage time computations.

Figure 1 shows the results from both MD (solid square symbols) and the theoretical model (solid lines). Figure 1A shows the dependence of the estimated first passage time on the applied strain for the three different temperatures studied, while Figure 1B presents the estimated time as a function of temperature. Fitting the MD results according to Eq. 2 clearly reproduces the trends computed from MD, where we observe a significant effect of strain on the first passage time, which may result in differences between the estimated times by almost two orders of magnitude. As we biaxially compress the material the time decreases as the diffusivity increases. The temperature effect seems to follow the usual Arrhenius relation for the cases studied in this work.

FIGURE 1

FIGURE 1. First passage time for a SIA to reach the free surface traveling a distance $\text{Δ}x=4.4$ nm as obtained with MD and the theoretical model as a function of (A) applied equibiaxial strain, and (B) inverse temperature.

3.2 Adatom Diffusion on a $\left(1\overline{1}0\right)$ Tungsten Surface

3.2.1 Dynamics

Running MD simulations of the dynamics of the adatom on the W$\left(1\overline{1}0\right)$ surface, we computed the mean-squared displacement ($\langle R^2\rangle$) and the diffusivity following Einstein’s relation in 2D, i.e., $D_{ad}=\langle R^2\rangle /4t$, where t is the time. We have again fitted the diffusivity following TST as in Eq. 2. In this case, the fitting procedure gives $D_0=44.35{\text{Å}}^2$ ps⁻¹, $\text{Δ}E=0.8$ eV, and ${\sigma }_a=17.75$ GPa. We observe again that the compressive biaxial strain increases the diffusivity, although the effect is weaker than in the bulk, being within one order of magnitude, as shown in Figures 2A,B. This activation stress is an average value over the x and y directions. Although the individual values might not be exactly the same, the most likely strain state on a PFC is equibiaxial (Dasgupta et al., 2019) and that is our focus here. We note that the model follows closely the MD data, capturing the effect of both strain and temperature (see Figure 2). We observe that the adatom hops follow mostly one of the two $\langle 111\rangle$ directions on the $\left(1\overline{1}0\right)$ plane, leading to a fairly isotropic diffusion. This is in agreement with previous results in the literature, (Chen and Ghoniem, 2013) showing that the energy of an adatom hopping transition is lower than those of other mechanisms such as the exchange mechanism or the crowdion mechanism on the $\left(1\overline{1}0\right)$ surface, as is the fact that the dependence of the surface adatom diffusivity on strain is weak.

FIGURE 2

FIGURE 2. Diffusion coefficient of an adatom in a W $\left(1\overline{1}0\right)$ surface as a function of (A) equibiaxial strain, and (B) inverse temperature.

3.2.2. Minimum Energy Path for Adatom Diffusion

To rationalize the results obtained from the computation of the mean-squared displacement, we have performed NEB simulations to compute the MEPs and the respective activation barriers at various applied strain levels. Figure 3 shows the results, where we observe that indeed the compressive equibiaxial strain leads to a decrease in the activation enthalpy. The value of $\text{Δ}E$ at zero strain is 0.68 eV, compared to the 0.8 eV value obtained from the dynamics runs. We attribute this difference of 0.12 eV in the values of $\text{Δ}E$ computed by the two methods to statistical errors in the fitting procedure, small anharmonic effects, but also the fact that other hopping mechanisms with higher activation barriers, such as the exchange and crowdion mechanisms, that are not examined in the NEB analysis may be occurring during the MD simulations. The maximum difference in the energy barrier is 0.08 eV between 0.008 and $-0.008$ strain. The obtained value for the activation energy ($\text{Δ}E=0.68$ eV) is somehow smaller than that obtained in previous studies, (Chen and Ghoniem, 2013) where the authors reported a value of 0.87 eV for the hopping transition from density functional theory NEB simulations. This suggests that the interatomic potential used in this study overestimates the adatom diffusivity on the $\left\{1\overline{1}0\right\}$ free surfaces in W. Moreover, we have obtained a monotonic relation between strain and activation enthalpy with tensile strain, although DFT results (Chen and Ghoniem, 2013) show a non-monotonic relation, with strain always increasing the activation enthalpy. In both cases, DFT and empirical potential, the dependence of the enthalpy barrier on strain is weak.

FIGURE 3

FIGURE 3. Minimum energy path for an adatom hopping along a $\langle 111\rangle$ direction on a W$\left(1\overline{1}0\right)$ surface at various applied equibiaxial strain levels.

3.3 Surface Roughness Evolution

The third set of MD simulations examines the evolution of the surface morphology as added SIAs diffuse to the surface and migrate on it to form clusters and aggregates that lead to an increase in surface roughness depending on the applied equibiaxial strain. As shown in Section 3. One, strain has a profound effect on the SIA transport in the bulk material. The higher the tensile strain the longer the time the SIA spends inside the material, increasing the defect concentration, and hence, the propensity to react with incoming SIAs. Indeed, we see that the formation of SIA clusters and dislocation loops is more probable as the tensile strain increases. We have used the root mean squared surface roughness (SR), which is defined consistently with Ref. (Dasgupta et al., 2019)

$SR\left(t\right)=\sqrt{\frac{1}{N^2}{\displaystyle \sum }_{i=1}^N{\displaystyle \sum }_{j=1}^N{\left(h\left(i,j,t\right)-\overline{h}\left(t\right)\right)}^2}(3)$

with $\overline{h}\left(t\right)$ the average height

$\overline{h}\left(t\right)=\frac{1}{N^2}{\displaystyle \sum }_{i=1}^N{\displaystyle \sum }_{j=1}^{N}h\left(i,j,t\right)(4)$

We have run the MD simulations for 2 μs, adding a total of 10,000 SIAs to the W in the slab supercell. To compute the SR values, we have discretized the domain in x and y in small cells, with indices i and j, and found the atoms belonging to each of the surfaces. The size of the cells was chosen such that the minimum number of atoms per cell is one, with cell sizes ranging between 9 and 25 Å². For each surface, the average value of the z coordinates of the atoms in each cell was computed and stored as h(i,j,t). The average SR of both surfaces was then calculated and it is shown in Figure 4, where we present a five-point average in time of the SR evolution.

FIGURE 4

FIGURE 4. Surface roughness evolution at various levels of applied equibiaxial strain.

Initially, we see that equibiaxial compression leads to a slightly larger SR since the SIAs travel fast to the surface. However, as we add more atoms, the SR at compressive strains plateaus. At this point, extra atoms create new bcc layers under compression and do not increase the SR, while under tension, the SR continues to grow, creating a crossover at around 1 μs. The same trend is followed by the number of atoms at the surface, the surface area and the distance between free surfaces as obtained by the alpha-shape algorithm (Stukowski, 2014) as implemented in the Ovito software (Stukowski, 2010) using 3 Å as the sphere probe radius. The latter is shown in Figure 5 after 2 μs, where we note that the average distance between free surfaces decays fairly linearly with the applied equibiaxial strain level. This indicates that the surface grows faster under compressive rather than tensile strain. It is also worth mentioning that the 2 μs reached by our MD simulations are too short compared with the experimental observations in which exposure times reach hundreds of minutes (i.e., several hours) and the surface roughness may reach values on the order of tens of nm.

FIGURE 5

FIGURE 5. Average distance between free surfaces located by the alpha-shape algorithm (Stukowski, 2014) implemented in Ovito (Stukowski, 2010) using a sphere probe of 3 Å radius as a function of the applied equibiaxial strain level.

Figure 6 shows the time evolution of the W slab in the simulation supercell when an equibiaxial tensile strain $\epsilon =0.008$ is applied. Only atomic arrangements with structures different from that of a bcc crystalline lattice following a common neighbor analysis (Stukowski, 2010) are presented and colored according to their potential energy. Dislocation lines are shown in light green as obtained by the dislocation-extraction algorithm, (Stukowski and Albe, 2010) that identifies the generated dislocation as edge in character with Burgers vector $\frac{1}{2}\langle 111\rangle$. As explained above, tensile strain increases the concentration of SIAs in bulk and enhances the propensity for cluster formation, eventually leading to the formation of edge dislocations as presented in Figure 6. Figure 6A displays the structure after 1 μs, where we observe a large density of dislocations and SIA clusters. Dislocations are close to the surface of the slab, and at 1.032 μs (Figure 6B) one dislocation segment reaches the surface, generating a depletion or groove on the surface, that increases the surface roughness. Surface diffusion leads to the flattening of the depletion and, hence, grooves are unstable at this stage. As the structural evolution proceeds, dislocations annihilate at the surface leading to the configuration shown in Figure 6C after 1.6 μs, where two dislocation segments impinge on the surface. As the top dislocation is absorbed (Figure 6D) at a time of 1.8 μs, some extra depletions are formed. The figure shows the evolution of the surface morphology as SIAs and dislocations reach the free surfaces.

FIGURE 6

FIGURE 6. Evolution at different simulated times denoted in each subfigure, of a W slab with $\left(1\overline{1}0\right)$ surfaces under equibiaxial tensile strain ${\epsilon }_{11}={\epsilon }_{22}=0.008$ showing the presence of edge dislocations (green lines) generated as aggregates of SIAs and annihilating at the surface, forming grooves that dissolve due to diffusive processes.

4 Discussion

The SIA concentration profile at steady state in the kinetic regime dominated by annihilation at the free surfaces is obtained from the diffusion equation $D_{SIA}\left(\epsilon \right)[{\partial }^2{c}_{SIA}/\partial z^2]=-G$, where the recombination term is neglected. $D_{SIA}\left(\epsilon \right)$ is the SIA diffusion coefficient that depends on the applied strain as derived above and G is a source term for SIA production related to the He implantation (taken as $G={10}^{-11}1/\text{s}/{\text{A}}^3$). Imposing boundary conditions $c_{SIA}\left(0\right)=0$ and $\partial c_{SIA}/\partial z\left(L/2\right)=0$ in a 1D system of length L, we obtain

$c_{SIA}\left(z\right)=\frac{G}{2D_{SIA}}z\left(L-z\right)(5)$

The former boundary condition implies that the free surface at z = 0 acts as a sink for SIAs while the latter one is a symmetry condition in the middle of the W slab. The steady-state concentration depends on strain as the diffusivity $D_{SIA}$ does. Figure 7 shows the results for five different applied equibiaxial strains ranging from −0.02 to 0.02. We observe that as the diffusivity varies with varying applied strain, including changes of the strain state from tensile to compressive, the SIA concentration in the slab increases, reaching a difference of almost five orders of magnitude for the strains examined in this study. This SIA concentration difference reaches three orders of magnitude compared to the case where strain is not considered $\left(\epsilon =0\right)$.

FIGURE 7

FIGURE 7. Steady-state concentration of SIAs in a material slab of 90 Å in thickness obtained by the solution to the diffusion equation given by Eq. 5 at different levels of biaxial strain.

It is important to note that under reactor operating conditions, a PFC is expected to be subject to compressive equibiaxial strains induced by the presence of over-pressurized He bubbles in the PFC near-surface region. However, the generation of a significant density of edge dislocations as a result of aggregation of SIAs or directly emitted by the bubbles, through the loop punching process, may lead to more complex scenarios, with local strain fields that can vary in character. Moreover, the strain state is likely to also change binding energies, which might modify the source terms, that we have however kept constant throughout this study, as we focus on the transport properties. Furthermore, surface growth and fuzz formation under biaxial compressive stress and promotion of surface roughness under biaxial tensile stress are not in contradiction with each other. We know that, due to the formation of over-pressurized helium bubbles, the near-surface region of PFC tungsten is under an equibiaxial compressive state of stress and, according to our surface morphological evolution models (see Refs. (Dasgupta et al., 2019; Dasgupta et al., 2020; Chen et al., 2021)), this is consistent with the formation and growth of nanotendrils emanating from the PFC surface, which is known to be the precursor to fuzz formation in PFC tungsten. In other words, the development of biaxial compressive stress is expected to promote fuzz formation. This is not at odds with the promotion of surface roughness by tensile stress. For example, it is well known that tensile stress (as well as compressive stress) can induce surface morphological instabilities, such as the Asaro-Tiller-Grinfeld instability (Maroudas, 2011) that promote surface roughness and may lead to surface cracking. The results presented here intend to convey a general picture of the effects of strain on the transport of SIAs and its effect on the morphological evolution of the free surface, in conjunction with providing activation parameters that can be readily used in larger-scale models, such as in cluster dynamics or reduced-order frameworks. (Blondel et al., 2018; Dasgupta et al., 2019)

As a matter of fact, the reduced-order (continuum-scale) model proposed by Dasgupta et al. (Dasgupta et al., 2019) describes the surface evolution in terms of the SIA fluxes towards the PFC surfaces and the adatom surface diffusion

$\frac{\partial h}{\partial t}=H\text{′Ω}{\nabla }_s\cdot {\mathbf{J}}_s+\text{Ω}{J}_{SIA}(6)$

with $H\prime ={[1+{\left(\frac{\partial h}{\partial x}\right)}^2+{\left(\frac{\partial h}{\partial y}\right)}^2]}^{\frac{1}{2}}$, ${\nabla }_s$ the surface gradient operator, $\text{Ω}$ the atomic volume for W, ${\mathbf{J}}_s$ the surface mass flux, and $J_{SIA}$ the SIA flux from the bulk to the surface, which can be written, following classical phenomenology (Nernst-Einstein equation), as

$J_{SIA}=-\frac{D_{SIA}{c}_{SIA}}{k_{B}T}\frac{\partial {\mu }_{SIA}}{\partial z}(7)$

At steady state, $D_{SIA}{c}_{SIA}$ is constant as all the SIAs generated reach the surface, which can be easily seen substituting Eq. 5 into Eq. 7. However, as previously shown, the value for the concentration varies significantly leading to different clustering propensities and the formation of dislocation loops that will modify the microstructure and the surface morphology. This is a direct consequence of the fact that, in a mean-field approach, the clustering probability is proportional to the square of the concentration of the diffusing species. Also note that the He content in the nanobubble region reaches a steady state (Woller et al., 2015), but the surface morphology does not and keeps evolving under the action of the biaxial compressive strain.

The surface flux can be written as ${\mathbf{J}}_s=\frac{D_s{\delta }_s}{\text{Ω}{k}_{B}T}{\nabla }_s\left(\mu \right)$ with μ the chemical potential that depends on the surface curvature and elastic strain energy, $k_{B}T$ the Boltzmann factor, and ${\delta }_s/\text{Ω}$ the number of surface atoms per unit area. $D_s$ depends again on the local strain as described in Section 3.2, although the strain effect is expected to be significantly lower than for the bulk case. This difference in the relevance of the strain effect on the two diffusion processes is attributed to the significantly lower activation barrier for SIA diffusion as compared to that for adatom diffusion on the surface. However, the effect on adatom diffusion is far from negligible, as highlighted in Figures 2, 3, and should be accounted for in coarse-grained models of surface morphological evolution.

5 Conclusions

In summary, we have analyzed the effect of equibiaxial strain on the transport of self-interstitial atoms towards a $\left(1\overline{1}0\right)$ surface and of an adatom at this surface in tungsten relying on a molecular dynamics methodology. We observe a strong effect of strain on the bulk transport properties, with a large activation stress that significantly modifies the concentration of SIAs. Such a change in SIA concentration leads to strong variation on the clustering propensity and, hence, the formation of dislocation loops. For tensile equibiaxial strain the probability of loop formation is higher than under compression, which leads to the formation of edge dislocation loops. These dislocation segments eventually reach the free surface and generate large atomic depletions or grooves that are slowly filled in and dissipate from the surface through surface diffusion. The strain effect on adatom surface diffusion is considerably weaker than in the bulk, with a monotonic relation between strain and activation enthalpy. Compressive and tensile equibiaxial strains lead to faster and slower transport, respectively. We have also computed the evolution of the surface roughness, showing a crossover at around 1 μs, with compression and tension leading to a rougher surface morphology before and after that time instant, respectively.

Data Availability Statement

The raw data supporting the conclusion of this article will be made available by the authors, without undue reservation.

Author Contributions

EM: Conceptualization, Methodology, Software, Writing-Reviewing and Editing. BS: Methodology, Reviewing and Editing. BW: Reviewing and Editing. DM: Methodology, Writing-Reviewing and Editing.

Conflict of Interest

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Acknowledgments

BS and EM thank support from the startup package at Clemson University. Clemson University is acknowledged for generous allotment of compute time on Palmetto cluster. EM (partial support), DM and BW acknowledge support from the U.S. Department of Energy, Office of Science, Office of Fusion Energy Sciences and Office of Advanced Scientific Computing Research through the Scientific Discovery through Advanced Computing (SciDAC) project on Plasma Surface Interactions under Award Nos. DE-SC0008875 and DE-SC0018421.

References

Baldwin, M. J., Doerner, R. P., Nishijima, D., Tokunaga, K., and Ueda, Y. (2009). The Effects of High Fluence Mixed-Species (Deuterium, Helium, Beryllium) Plasma Interactions with Tungsten. J. Nucl. Mater. 390-391, 886–890. doi:10.1016/j.jnucmat.2009.01.247

CrossRef Full Text | Google Scholar

Blondel, S., Bernholdt, D. E., Hammond, K. D., and Wirth, B. D. (2018). Continuum-scale Modeling of Helium Bubble Bursting under Plasma-Exposed Tungsten Surfaces. Nucl. Fusion 58 (12), 126034. doi:10.1088/1741-4326/aae8ef

CrossRef Full Text | Google Scholar

Chen, C.-S., Dasgupta, D., Weerasinghe, A., Wirth, B. D., and Maroudas, D. (2021). Effects of Elastic Softening and Helium Accumulation Kinetics on Surface Morphological Evolution of Plasma-Facing Tungsten. Nucl. Fusion 61, 016016. doi:10.1088/1741-4326/abbf64

CrossRef Full Text | Google Scholar

Chen, Z., and Ghoniem, N. (2013). Biaxial Strain Effects on Adatom Surface Diffusion on Tungsten from First Principles. Phys. Rev. B 88 (3), 035415. doi:10.1103/physrevb.88.035415

CrossRef Full Text | Google Scholar

Dasgupta, D., Kolasinski, R. D., Friddle, R. W., Du, L., Maroudas, D., and Wirth, B. D. (2019). On the Origin of ‘fuzz’ Formation in Plasma-Facing Materials. Nucl. Fusion 59, 086057. doi:10.1088/1741-4326/ab22cb

CrossRef Full Text | Google Scholar

Dasgupta, D., Maroudas, D., and Wirth, B. D. (2020). Prediction of Temperature Range for the Onset of Fuzz Formation in Helium-Plasma-Implanted Tungsten. Surf. Sci. 698, 121614. doi:10.1016/j.susc.2020.121614

CrossRef Full Text | Google Scholar

Fiflis, P., Curreli, D., and Ruzic, D. N. (2015). Direct Time-Resolved Observation of Tungsten Nanostructured Growth Due to Helium Plasma Exposure. Nucl. Fusion 55 (3), 033020. doi:10.1088/0029-5515/55/3/033020

CrossRef Full Text | Google Scholar

Heikinheimo, J., Mizohata, K., Räisänen, J., Ahlgren, T., Jalkanen, P., Lahtinen, A., et al. (2019). Direct Observation of Mono-Vacancy and Self-Interstitial Recovery in Tungsten. APL Mater. 7 (2), 021103. doi:10.1063/1.5082150

CrossRef Full Text | Google Scholar

Henkelman, G., and Jónsson, H. (2000). Improved tangent Estimate in the Nudged Elastic Band Method for Finding Minimum Energy Paths and Saddle Points. J. Chem. Phys. 113 (22), 9978–9985. doi:10.1063/1.1323224

CrossRef Full Text | Google Scholar

Henkelman, G., Uberuaga, B. P., and Jónsson, H. (2000). A Climbing Image Nudged Elastic Band Method for Finding Saddle Points and Minimum Energy Paths. J. Chem. Phys. 113 (22), 9901–9904. doi:10.1063/1.1329672

CrossRef Full Text | Google Scholar

Honeycutt, J. D., and Andersen, H. C. (1987). Molecular Dynamics Study of Melting and Freezing of Small Lennard-Jones Clusters. J. Phys. Chem. 91, 4950–4963. doi:10.1021/j100303a014

CrossRef Full Text | Google Scholar

Ihli, T., Hermsmeyer, S., Köhly, C., and Norajitra, P. (2006). Integration of an Advanced He-Cooled Divertor in a DEMO-Relevant Tokamak Geometry. Fusion Eng. Des. 81, 121–126. doi:10.1016/j.fusengdes.2005.05.002

CrossRef Full Text | Google Scholar

Ito, A. M., Takayama, A., and Nakamura, H. (2018). Triple Hybrid Simulation Method for Tungsten Fuzzy Nanostructure Formation. Plasma Fusion Res. 13 (0), 3403061. doi:10.1585/pfr.13.3403061

CrossRef Full Text | Google Scholar

Ito, A. M., Takayama, A., Oda, Y., Tamura, T., Kobayashi, R., Hattori, T., et al. (2015). Molecular Dynamics and Monte Carlo Hybrid Simulation for Fuzzy Tungsten Nanostructure Formation. Nucl. Fusion 55 (7), 073013. doi:10.1088/0029-5515/55/7/073013

CrossRef Full Text | Google Scholar

Kajita, S., Sakaguchi, W., Ohno, N., Yoshida, N., and Saeki, T. (2009). Formation Process of Tungsten Nanostructure by the Exposure to Helium Plasma under Fusion Relevant Plasma Conditions. Nucl. Fusion 49 (9), 095005. doi:10.1088/0029-5515/49/9/095005

CrossRef Full Text | Google Scholar

Knaster, J., Moeslang, A., and Muroga, T. (2016). Materials Research for Fusion. Nat. Phys 12 (5), 424–434. doi:10.1038/nphys3735

CrossRef Full Text | Google Scholar

Krasheninnikov, S. I. (2011). Viscoelastic Model of Tungsten ‘fuzz’ Growth. Phys. Scr. T145, 014040. doi:10.1088/0031-8949/2011/t145/014040

CrossRef Full Text | Google Scholar

Lasa, A., Tähtinen, S. K., and Nordlund, K. (2014). Loop Punching and Bubble Rupture Causing Surface Roughening -A Model for W Fuzz Growth. Epl 105 (2), 25002. doi:10.1209/0295-5075/105/25002

CrossRef Full Text | Google Scholar

Marinica, M.-C., Ventelon, L., Gilbert, M. R., Proville, L., Dudarev, S. L., Marian, J., et al. (2013). Interatomic Potentials for Modelling Radiation Defects and Dislocations in Tungsten. J. Phys. Condens. Matter 25 (39), 395502. doi:10.1088/0953-8984/25/39/395502

PubMed Abstract | CrossRef Full Text | Google Scholar

Maroudas, D. (2011). Surface Morphological Response of Crystalline Solids to Mechanical Stresses and Electric Fields. Surf. Sci. Rep. 66, 299–346. doi:10.1016/j.surfrep.2011.05.001

CrossRef Full Text | Google Scholar

Martynenko, Y. V., and Nagel’, M. Y. (2012). Model of Fuzz Formation on a Tungsten Surface. Plasma Phys. Rep. 38 (12), 996–999. doi:10.1134/s1063780x12110074

CrossRef Full Text | Google Scholar

Plimpton, S. (1995). Fast Parallel Algorithms for Short-Range Molecular Dynamics. J. Comput. Phys. 117, 1–19. doi:10.1006/jcph.1995.1039

CrossRef Full Text | Google Scholar

Redner, S. (2001). A Guide to First-Passage Proccesses. Cambridge: Cambridge University Press.

Stukowski, A., and Albe, K. (2010). Extracting Dislocations and Non-dislocation crystal Defects from Atomistic Simulation Data. Model. Simul. Mater. Sci. Eng. 18 (8), 085001. doi:10.1088/0965-0393/18/8/085001

CrossRef Full Text | Google Scholar

Stukowski, A. (2014). Computational Analysis Methods in Atomistic Modeling of Crystals. JOM 66 (3), 399–407. doi:10.1007/s11837-013-0827-5

CrossRef Full Text | Google Scholar

Stukowski, A. (2010). Visualization and Analysis of Atomistic Simulation Data with OVITO-The Open Visualization Tool. Model. Simul. Mater. Sci. Eng. 18 (1), 015012. doi:10.1088/0965-0393/18/1/015012

CrossRef Full Text | Google Scholar

Suzuki, S., Ueda, Y., Tokunaga, K., Sato, K., and Akiba, M. (2003). Present Research Status on Divertor and Plasma Facing Components for Fusion Power Plants. Fusion Sci. Technol. 44, 41–48. doi:10.13182/fst03-a308

CrossRef Full Text | Google Scholar

Takamura, S., Ohno, N., Nishijima, D., and Kajita, S. (2006). Formation of Nanostructured Tungsten with Arborescent Shape Due to Helium Plasma Irradiation. Plasma Fusion Res. 1, 051. doi:10.1585/pfr.1.051

CrossRef Full Text | Google Scholar

Wang, K., Doerner, R. P., Baldwin, M. J., Meyer, F. W., Bannister, M. E., Darbal, A., et al. (2017). Morphologies of Tungsten Nanotendrils Grown under Helium Exposure. Sci. Rep. 7, 42315. doi:10.1038/srep42315

PubMed Abstract | CrossRef Full Text | Google Scholar

Wirth, B. D., Hammond, K. D., Krasheninnikov, S. I., and Maroudas, D. (2015). Challenges and Opportunities of Modeling Plasma-Surface Interactions in Tungsten Using High-Performance Computing. J. Nucl. Mater. 463, 30–38. doi:10.1016/j.jnucmat.2014.11.072

CrossRef Full Text | Google Scholar

Woller, K. B., Whyte, D. G., and Wright, G. M. (2015). Dynamic Measurement of the Helium Concentration of Evolving Tungsten Nanostructures Using Elastic Recoil Detection during Plasma Exposure. J. Nucl. Mater. 463, 289–293. doi:10.1016/j.jnucmat.2014.11.126

CrossRef Full Text | Google Scholar

Yang, Q., You, Y.-W., Liu, L., Fan, H., Ni, W., Liu, D., et al. (2015). Nanostructured Fuzz Growth on Tungsten under Low-Energy and High-Flux He Irradiation. Sci. Rep. 5, 10959. doi:10.1038/srep10959

PubMed Abstract | CrossRef Full Text | Google Scholar