Helium bubble size effects on the surface morphological response of plasma-facing tungsten

We report a simulation study on the effects of helium (He) bubble size on the morphological evolution and pattern formation on the surface of tungsten used as a plasma-facing component (PFC) in nuclear fusion devices. We have carried out a systematic investigation based on self-consistent dynamical simulations of surface morphological evolution according to an atomistically-informed, 3D continuum-scale model that captures well the relevant length and time scales of surface nanostructure formation in PFC tungsten. The model accounts for PFC surface diffusion, driven by the biaxial compressive stress originating from the over-pressurized He bubbles in the near-surface region of PFC tungsten as a result of He plasma exposure, combined with the formation of self-interstitial atoms in tungsten that diffuse toward the PFC surface and the flux of surface adatoms generated as a result of surface vacancy-adatom pair formation upon He implantation; this transport of surface adatoms contributes to the anisotropic growth of surface nanostructural features due to the different rates of adatom diffusion along and across step edges of islands on the tungsten surface. Our detailed analysis reveals that varying the average He bubble size plays an important role in the PFC surface growth kinetics as well as the resulting surface topography. Specifically, we find that an increase in the He bubble size leads to a deceleration in the growth rate of the tungsten nanotendrils that emanate from the PFC surface. We also find that the separation distance between the resulting surface features increases with increasing He bubble size, as well as over time. This coarsening effect is a thermally activated process resulting in an accurate description of the temperature dependence of the average surface feature separation by an Arrhenius relation.


Introduction
In the International Thermonuclear Experimental Reactor (ITER) and other future magnetically confined devices, the plasma-facing components (PFCs), such as the divertor, will experience extreme heat and particle fluxes during steady-state operation. These extreme conditions limit the choice for the PFC material which can withstand cumulative incident particle fluence and energy to facilitate long uninterrupted operation and even further extreme heat flux up to 1 GW m −2 under plasma transients, such as type I Edge Localized Modes (ELMs) [1]. Tungsten is the most reliable material to be used as a PFC due to its outstanding thermomechanical properties [2][3][4][5]. However, as a divertor PFC material, tungsten poses a unique challenge under low-energy (35 eV or above) helium (He) plasma exposure [6][7][8] and over a temperature range from 900 K to 2000 K, which are typical conditions expected in the ITER divertor. Under such conditions, tungsten readily forms a fragile crystalline nanotendril-like structure, commonly called 'fuzz' [6][7][8][9][10][11][12]; formation of nanotendrils emanating from the PFC surface constitutes a precursor to the growth of fuzz, a surface layer consisting of such interconnected nanotendrils with complex morphological features. Over the last two decades, a number of theoretical [13][14][15], simulation [16][17][18][19][20][21][22][23], and experimental [6,10,24,25] studies have been conducted to decipher the intriguing physics/mechanism of fuzz formation because, if it remains unmitigated, the fragile fuzz of the high-Z tungsten (W) can be exfoliated and released into the plasma and, thus, severely compromise the performance of the reactor [9][10][11][12]. Helium has a very low solubility in tungsten; as a result, upon implantation, He atoms readily self-cluster and nucleate helium bubbles, which subsequently grow by absorbing helium and emitting tungsten Frenkel pairs, a process known as the 'trap-mutation' reaction [24,26]. Helium bubbles in PFC tungsten are at an over-pressurized state and, consequently, they exert an equibiaxial compressive stress [22] on the tungsten matrix in the near-surface region of the PFC material.
In our previous studies, we have postulated that the competition between the elastic strain energy in the near-surface PFC tungsten region and the PFC surface free energy is the primary driving force for the early-stage growth of nanostructures on the plasma-exposed tungsten surface, acting as a precursor to fuzz formation [22]. We have developed a continuum-scale model of PFC tungsten surface morphological response based on the aforementioned hypothesis, benchmarked the model by comparing its predictions with experimental results [22], and further characterized the early stage of fuzz growth as a surface morphological instability [22,27]. However, evolving from a precursor to fuzz to the actual fuzz formation involves a much more complex process, where dislocation loop punching [18,26,28], near-surface helium bubble evolution dynamics [29][30][31][32][33], helium bubble bursting, and, as a result of that, surface crater formation [18,[33][34][35] together play an important role in the formation of PFC surface nanostructure. To model such complex phenomena, governed by disparate spatiotemporal scales, and implement such a model to simulate PFC surface morphological response over experimentally relevant time (hours) and length (micrometers) scales [13,14,[19][20][21][36][37][38][39], we have followed a hierarchical multiscale modeling strategy for model development. Using our model, we have systematically studied the effects of elastic softening due to the presence of helium bubbles in the PFC near-surface region [40,41], surface hole formation due to bubble bursting [42], change in He bubble pressure due to increasing the tungsten surface temperature [43], and surface vacancy-adatom pair formation [44] on the PFC surface morphological response. In these studies, the average size of the helium bubble in the PFC near-surface region plays an important role in determining the PFC surface morphological evolution; in all our previous studies [22,27,[41][42][43][44], based on a very limited set of experimental observations available in the literature [12,45], we considered the maximum value reached by the average He bubble radius when the He content in the nearsurface region reaches a steady state to be in the range of 0.75 nm to 1.5 nm.
However, a recent experimental study [46] has reported a broader size range of He bubbles in the nearsurface region of PFC tungsten, with the average bubble radii extending up to 3 nm. The average He bubble size in the PFC near-surface region is dependent on the irradiation conditions [12,47] and an exact dependence has not yet been established. Therefore, in the present study, we treat the average He bubble radius in the PFC nearsurface region at steady-state He content as a model parameter. Through such a parametric study, in this article, we have thoroughly investigated the combined effects of He bubble size and surface temperature on the surface morphological response of PFC tungsten. We find that a change in the He bubble size has a direct impact on the growth kinetics of the PFC surface nanostructure as well as on the characteristics of the surface topography during the early stage of the fuzz formation process.

Morphological evolution model for PFC tungsten surface
In our continuum-scale model for the morphological evolution of a PFC tungsten surface exposed to a He plasma with an incident flux of 2.7 × 10 21 ions m −2 s −1 and a He ion incident energy of 75 eV, (which is within the range of incident energies expected for the divertor region in ITER and has been used in experimental studies including that of [22]), we divide the tungsten near-surface region into two layers. The top layer consists of a homogeneous distribution of nanometer-scale spherical He bubbles of uniform size and number density, with an isotropic He/W interfacial free energy for the bubble/matrix interface, and has an average thickness on the order of tens of nanometers, while the bottom layer is practically free of He ions and is considered as perfect and pristine bulk tungsten; the number density of He bubbles is set, as detailed in [22], by the He concentration in the nanobubble region and the amount of He in each bubble determined by the He/vacancy ratio in the bubble [22]. The continuity equation for the conservation of tungsten atoms that governs the surface morphological evolution of PFC tungsten is expressed as where h = h(x, y, t) is the surface height function with x, y, and z being the usual Cartesian coordinates, , ∇ s is the surface gradient operator, and Ω is the atomic volume of tungsten at temperature T. The tungsten surface atomic (mass) flux, J s , the self-interstitial atomic flux toward the surface along the z-direction, J I , the step-edge diffusive current, J E [48][49][50], and the terrace diffusive current, J T [48][49][50], together provide the driving forces that facilitate the surface feature formation process on the PFC tungsten surface. J s is proportional to the gradient of the chemical potential, μ s , of the tungsten atoms present on the surface. This chemical potential, μ s , consists of contributions from the free energy of the curved tungsten surface and the elastic strain energy, generated due to the over-pressurized He bubbles in the top layer (or nanobubble region) of the near-surface region of PFC tungsten. We have approximated the dislocation loop formation and loop punching process from helium bubbles with a flux, J I , of self-interstitial atoms (SIAs), which is proportional to the gradient of the chemical potential, μ I , of the SIAs along the z-direction. J E and J T are generated as a result of the surface vacancy-adatom pair formation effect and are described and discussed in full detail in [44]. We have neglected the effects of tungsten sputtering from the PFC surface because the He incident ion energy is less than the sputtering threshold (∼150 eV) for tungsten [22,29]. The surface atomic flux, step-edge diffusive current, and terrace diffusive current are strongly dependent on the surface, step-edge, and terrace diffusivities, respectively, which along with the surface adatom density and surface free energy per unit area are functions of the crystallographic orientation of the plasma-exposed surface. The above surface evolution model and its parameterization is described in more detail in [22,27,[41][42][43][44]. We emphasize that helium bubble growth in the governing equation, equation (1), is considered by the direct dependence on the He bubble radius of the biaxial stress in the near-surface nanobubble region of PFC tungsten, from the beginning of plasma exposure to He content saturation where the over-pressurized He bubbles reach their steady-state size; this radius-dependent stress level determines the contribution of the elastic strain energy to the surface chemical potential, μ s , which has been detailed in the description of the model in [22].

Numerical simulation techniques
We performed dynamical simulations using a 20 μm × 20 μm tungsten surface supercell to represent the boundary of the nanobubble region exposed to the plasma. The PFC tungsten surface (approximated by the above supercell) was discretized using a 512 × 512 grid, based on a convergence test, to optimize the computational cost for given numerical accuracy [41]. The surface height at each grid point (i,j) represents the local surface morphology obtained from the dynamical simulations in the form of the surface height function h (i, j, t). The thickness of the nanobubble region was set to be h 0 = 30 nm [51] and periodic boundary conditions were applied in the x-and y-directions. We computed the stress field in the nanobubble layer numerically by solving the elastostatic boundary value problem (BVP) using a spectral collocation method and discrete fast Fourier transforms (FFTs). Coupled self-consistently with the elastostatic BVP, we used an operator-splitting semi-implicit front tracking method with adaptive time stepping to record the PFC surface morphological evolution [22,[52][53][54]. The initial tungsten surface configuration was prepared by randomly perturbing the surface from its planar morphological state at a root mean square (RMS) surface roughness of SR = 5 nm, which is consistent with the tungsten samples used in the He plasma exposure experiments in [22]. The 3D simulations of the surface morphological evolution of PFC tungsten were conducted at various temperatures, over the range from T = 700 K to T = 1100 K, and the average He bubble radius, at He saturation (i.e., when the He content in the nanobubble region reaches its steady state), was varied from r b,s = 1 nm to r b,s = 3 nm. Consistent with our previous work [44], we have used two types of surface morphological characterization metrics, namely, the average nanotendril height, h t tendril ( ), and the RMS surface roughness, SR(t). To determine the average nanotendril height, h t tendril ( ), we have identified local maxima of the surface height function as the heights of individual nanotendrils and taken the average value of the identified nanotendril heights over the entire simulation supercell.

Simulated PFC tungsten surface morphologies
We have conducted 3D dynamical simulations to predict the morphological response of a W(110) surface under helium-plasma irradiation with an incident He ion energy of 75 eV and He ion incident flux of 2.7 × 10 21 ions m −2 s −1 and characterized in detail the kinetics of the surface morphological evolution and the feature patterns forming on the surface over the temperature range from T = 700 K to T = 1100 K. Representative simulation results are reported in figures 1 and 2. The simulated plasma exposure conditions are consistent with those in the experiments used to benchmark our initial model [22] and the He incident flux is well within the particle flux range expected in the ITER divertor. In our simulations, helium plasma exposure is continued until the evolving tungsten surface reaches an RMS surface roughness, SR, of 50 nm, i.e., one order of magnitude higher than the initial surface roughness of the tungsten sample exposed to the plasma. The total helium concentration in the near-surface region of plasma-irradiated tungsten, which is referred to as He content, is approximated based on a first-order He accumulation kinetics and it reaches a saturation (steady-state) value of n He,s = 16 % [41] within a few minutes of plasma exposure, consistent with the experimental observations [51,55]. Figure 1 shows the top views of the PFC tungsten surface morphologies at the temperatures of ((a)-(d)) T = 700 K, ((e)-(h)) T = 900 K, and ((i)-(l)) T = 1100 K. The average He bubble radii used in the simulations are r b,s = 1nm, r b,s = 1.5 nm, r b,s = 2 nm, and r b,s = 3 nm as indicated in the surface morphologies depicted in figure 1. Figure 2 shows the corresponding cross-sectional surface profiles, h(x), along the x-direction at locations marked with the white solid lines in figures 1(a), (e), and (i). In figure 1, consistent with the findings of our recent study [44], we observe the formation of a surface pattern that consists of extended ridge-like surface features oriented along the vertical (y) direction, aligned with the [001] crystallographic direction, which we attribute to the difference in the (terrace step) edge-diffusion barriers along different surface crystallographic directions. We also observe that the separation distance between these y-aligned surface features increases with increasing temperature at given He bubble size and also increases with increasing He bubble size at given temperature. At low temperature, we observe a shorter average spacing between nanostructures formed on the surface because the vacancy-adatom pair formation effect, which gives rise to terrace and step-edge currents, is stronger than at higher temperature compared to the biaxial compressive stress-driven surface diffusion [44], which dominates at high temperature and small bubble sizes. Also, an increase in the average He bubble size causes a decrease in the bubble pressure and, thus, decreases the overall level of biaxial compressive stress, σ 0,s , exerted on the near-surface nanobubble region of PFC tungsten at He saturation, which is the main driving force for surface feature formation, resulting in coarser surface features. Specifically, as detailed in section 3.  In figure 3, we report the growth kinetics of the tendril-like surface features (the formation of which corresponds to the early stage of fuzz growth on the PFC surface) by tracking the evolution of the average nanotendril height for all the average helium bubble sizes at He saturation and tungsten surface temperatures examined in this study. The results establish that the growth rate of the nanotendrils that emanate from the PFC tungsten surface decreases with increasing He bubble size. We attribute this clear trend to the decrease with increasing bubble size in the level of biaxial compressive stress (σ 0,s ) acting on the nanobubble region, one of the key driving forces behind the PFC surface growth kinetics. At relatively higher temperature, the effect of the average He bubble size on the growth kinetics is much stronger than the effect at a lower temperature because, at lower temperature, surface vacancy-adatom pair formation has a more prominent effect on the surface feature pattern formation [44], as discussed above.

Analysis of surface morphological evolution of PFC tungsten
In all the simulations under different temperature conditions and for different average He bubble sizes, we have tracked the helium-implanted surface layer thickness measured as the sum of the average nanotendril height, h t tendril ( ), calculated with respect to the lowest point on the PFC surface, and the initial thickness h 0 of the nanobubble region. Figure 4 shows the simulated surface layer thickness evolution as a function of fluence, Φ. In all plots, black open circles represent the simulated surface layer thickness and red and blue solid lines represent the scaling law at an early stage of surface morphological evolution where the layer thickness is proportional to Φ 0.25 , and optimal fits to the simulation results after the early stage of surface morphological evolution, respectively. Figures 4((a)-(d)) show that, at the low temperature of T = 700 K, the simulated surface layer thickness evolves according to a sublinear growth mode after a minimum He plasma exposure duration of 100 minutes. Since fluence Φ is proportional to time t, we have chosen time (as commonly done in surface growth studies) as the preferred variable instead of fluence to explain the scaling laws that characterize the evolution of layer thickness growth. In this sublinear growth mode, as observed in figures 4((a)-(d)), the time dependence of the surface layer thickness evolution scales as t 0.5 , as demonstrated by the fitted blue solid lines, which is consistent with experimental measurements [8]. At T = 900 K and a He bubble radius of r b,s = 1nm, shown in figure 4(e), the surface layer thickness evolution follows a superlinear (exponential) growth mode after crossing the early stage of surface layer evolution at a He fluence of Φ ∼10 25 m −2 . However, as the average He bubble size increases, the surface layer thickness growth mode in the later stage of surface evolution shifts from superlinear (at the small bubble radius of r b,s = 1 nm) to sublinear following a time dependence that scales as t 0.5 , as seen in figures 4(f), (g), and (h). For surface evolution at T = 1100 K and an average He bubble size of 1 nm, shown in figure 4(i), the surface layer thickness exhibits a strong superlinear (exponential) growth mode, which is described very well by the blue solid line that represents exponential growth as discussed in detail in [27] at this (smallest) bubble size (r b,s = 1 nm). However, as the average He bubble size increases at this high surface temperature, the growth rate in the later stage of surface evolution shifts from exponential to sublinear following a time dependence of t 0.85 (at r b,s = 1.5 nm), where the exponent of 0.85 has no special physical meaning but is simply an optimal fitting parameter to this portion of the respective dataset, and t 0.5 (at r b,s = 2 and 3 nm), as seen in figures 4(j), (k), and (l), respectively. In almost all of the above cases (with the exception of that at the highest surface temperature and smallest He bubble size where the growth mode is exponential throughout the recorded surface evolution), at the very early stage of surface morphological evolution, the layer thickness growth is observed to have a sublinear behavior which scales approximately with t 0.25 , as shown by the red solid lines. This slower growth kinetics at the very early stage of PFC surface evolution is observed because the surface layer evolution is driven mainly by the terrace and step-edge fluxes, as described in detail in [44] and, thus, it is not affected by the He bubble size that controls the driving force of stress due to the over-pressurized He bubbles. However, at the later stages of surface evolution, the layer thickness growth kinetics is faster because it is dominated by the equibiaxial compressive stress induced by the over-pressurized He bubbles and, thus, it is sensitive to the average He bubble size. Therefore, an increase in the He bubble radius decreases the stress in the nanobubble region, which slows down the growth kinetics as seen in figures 4(j), (k), and (l). We emphasize that these PFC surface growth kinetics results, exhibiting a surface layer thickness evolution that scales with t 0.5 or, equivalently, Φ 0.5 are in good agreement with the experimental observations reported in [8]. Figure 5 shows the surface feature coarsening kinetics represented by the time evolution of the average surface feature separation along the x-direction of the W(110) surface, λ x . This surface feature separation, λ x , is calculated by performing a 2D FFT analysis following the approach of [42,44]. Missing data points within the plotting range of figures 5(a) and (b), for r b,s = 1 and 1.5 nm, respectively, correspond to cases (small bubble sizes and high temperature) where the surface morphological metrics grow beyond the predictive capabilities of our model, which is limited to the onset of nanotendril formation; this is because small bubble radii lead to the development of high biaxial compressive stress levels in the nanobubble region, which constitute the main drivers of plasma-exposed surface growth at the later stage of plasma exposure. In all cases, we observe that the surface features coarsen over time and (except for simulations at the highest temperatures and smallest He bubble radius examined) their separations reach steady-state values, which increase with increasing temperature. To characterize this temperature dependence, we plot the obtained λ x values at different time instants (He plasma exposure durations) over the temperature range examined, from T = 700 K to T = 1100 K, and average He bubble sizes of r b,s = 1 nm, 1.5 nm, 2 nm, and 3 nm in figures 6(a), (b), (c), and (d), respectively. In all cases, λ x decreases with increasing inverse temperature (i.e., increases with increasing temperature), while it increases as the duration of the PFC surface exposure to a He plasma is increased, which further clarifies the coarsening effect of the surface features. For a more detailed analysis of the dependence of λ x on temperature, we demonstrate that the simulation results are described very well by optimal fits according to an Arrhenius relation, exp where ΔE is a thermal activation barrier for the coarsening process. A consistent activation barrier of approximately 0.2 eV across all the simulation datasets confirms that the coarsening effect is a thermally activated process, and this is consistent with experimental findings reported in [55]. The thermal activation barriers calculated by fitting the Arrhenius relation to the simulation data, shown in figure 6, for different average He bubble sizes and durations of He plasma exposure, are listed in table 1.

Morphological stability analysis
To fundamentally understand the origin of the mechanism responsible for the surface pattern formation on a PFC tungsten surface exposed to a He plasma (specifically for this study, to a He ion incident energy of 75 eV and a He incident flux of 2.7 × 10 21 ions m −2 s −1 ), we performed a linear stability analysis by disturbing the planar tungsten surface morphology through a low-amplitude plane-wave perturbation according to the approach described in Refs. [44,52,54]. From this linear stability analysis, we obtain the maximally unstable wavelength, max l , of the PFC tungsten surface, which provides an upper bound for the surface feature separation distance, λ x , predicted from our dynamical simulations. The maximally unstable wavelength is expressed in dimensionless form, where lengths are made dimensionless with the characteristic surface length scale l [43], as  In all plots, the colored solid lines denote the optimal fits to the corresponding datasets according to an Arrhenius relation. The error bars are from statistical analysis of five independent simulations for each parameter set, corresponding to different initial surface configurations. Table 1. Thermal activation barrier values in eV corresponding to optimal fits of the datasets of figure 6 according to Arrhenius relations. The last column ( max l ) represents the derived barrier from optimal fitting of the LST-predicted values for the maximally unstable wavelength, max l , according to Arrhenius relations. ], and l is the characteristic surface length scale, expressed as l M s 0 0, 2 g s º [43]; in the above equations, δ s /Ω is the area density of surface atoms, M 0 is the biaxial modulus of tungsten, γ is the tungsten surface free energy per unit area, and σ 0,s is the equibiaxial compressive stress level in the nanobubble region at He saturation. These results for max l are plotted in figure 6 and are in good agreement with the simulation results (indicating comparable activation barriers while providing an upper bound to λ x ) and our previous findings [22,27,41,42,44]. By using leading-order approximations to weaker (than Arrhenius) temperature dependences of material properties and employing the binomial expansion and its proper truncation at low values of the dimensionless parameter 9 8 T x X ( ), we have expressed max l as the function of temperature (polynomial of a Boltzmann factor) In equation (3), exp = -, E s is the activation barrier for surface diffusion, E t is the activation barrier for terrace diffusion, and A, B, C, and D are dimensionless constants. At low temperatures and large He bubble sizes, ε = 1 and equation (3) can be simplified to the single Arrhenius equation , with a thermal activation barrier of E A /2 ; 0.40 eV, which is comparable to the barrier obtained from fitting according to an Arrhenius equation the linear stability theory (LST) predictions for max l at the low-T and large-r b,s asymptotic limit (0.35 eV for this limit in figure 6(d)). This analysis provides theoretical support for the temperature dependence of the feature separation distance λ x according to an Arrhenius relation, as demonstrated by the simulation results of figure 6.

Summary and conclusions
In summary, we have conducted a detailed computational and theoretical analysis of the effects of He bubble size on the surface morphological evolution of plasma-facing tungsten focusing on the W(110) surface orientation. We have found that varying the average He bubble radius in the near-surface nanobubble region of PFC tungsten has a direct impact on the kinetics of the surface morphological response of the PFC material and the resulting surface topography. Specifically, increasing the He bubble size leads to a decrease in the biaxial compressive stress in the nanobubble region, which is the main driving force for nanotendril formation and growth on the PFC surface. Thus, the growth rate of the nanotendrils emanating from the PFC tungsten surface decreases with increasing He bubble size. Over most of the range of He bubble sizes and temperatures examined, the growth kinetics of the PFC surface layer thickness exhibits a time dependence that is proportional to the square root of time (or incident He fluence from the plasma exposure), which is consistent with experimental data. We also found that the separation distance between the resulting surface features forming on the PFC surface upon He plasma exposure increases with increasing average He bubble size and duration of He plasma exposure. We have analyzed this coarsening effect systematically and found that it is a thermally activated process with the average surface feature separation depending on temperature according to an Arrhenius relation. Moreover, the simulation results are in good agreement with theoretical predictions based on linear stability analysis of the PFC surface morphological response, which provides analytical support and a fundamental interpretation to the simulation data.
The continuum-scale model used in this study to examine the surface morphological evolution of PFC tungsten accounts for many important phenomena that occur in PFC tungsten such as softening of the elastic moduli of the PFC material, generation of equibiaxial stress due to over-pressurized He bubbles in the nearsurface nanobubble region of PFC tungsten, as well as anisotropic diffusion of surface adatoms generated through surface vacancy-adatom pair formation on the plasma-exposed tungsten surface. However, the fuzz formation process on the PFC tungsten surface is also mediated by other important dynamical phenomena, related to He bubbles formed in the near-surface region of the PFC material, such as dynamical bubble formation, coalescence of He bubbles, and He bubble bursting. Consequently, a fundamental understanding of the effects of He bubble size on the morphological evolution of the PFC tungsten surface will help us incorporate these additional complex physics (bubble dynamics) into our atomistically-informed continuum-scale model toward enabling more accurate predictive capabilities through a more direct bubble-physics-based modeling, for which efforts are currently in progress. Finally, we mention that in this study we focused on the W(110) surface crystallographic orientation for consistency with our previous studies on surface nanotendril formation and using the W(110) surface as prototypical given its lowest surface free energy per unit area, and, thus, highest thermodynamic stability, among low-index tungsten surfaces. We are currently working on the morphological evolution of other tungsten surfaces with different low-index crystallographic orientations, such as W(001), which will help us evaluate the impact of grain orientation on the surface growth kinetics of plasma-exposed polycrystalline tungsten.