Tungsten surface evolution by helium bubble nucleation, growth and rupture

Producing electricity from nuclear fusion is one of the greatest scientific and engineering challenges of the 21st century [1]. Magnetic confinement reactors, such as ITER and DEMO, require advances in both plasma physics technology and materials science, including structural and functional materials and the materials surfaces exposed to the plasma-material interface. Tungsten is a leading candidate to satisfy the varied demands of the plasma facing divertor due to its high melting point, low sputtering yield, high thermal conductivity, low transmutation probability, and low tritium retention [2–4].

Recent laboratory experiments performed in linear plasma devices indicate the possibility of substantial surface modification in tungsten exposed to low-energy mixed helium–hydrogen plasma, although the observed surface response is strongly temperature-dependent. Pitted surfaces are observed below ≈1000 K [5], whereas 'nanostructured' low-density 'fuzz' or 'coral' structures are observed between approximately 1000 and 2000 K [6–9], and micrometre-sized holes form above about 2000 K [10, 11]. The nanostructured 'fuzz' has recently been observed in the divertor regions of a tokamak device as well [12]. Such surface features could lead to changes in heat transfer, fuel (deuterium/tritium) retention [13], increased rates of erosion through both sputtering and dust formation [14], and embrittlement of the divertor, all of which can be detrimental to the plasma [15]. Transmission electron microscopy suggests that the nanoscale tendrils of the fuzz contain gas bubbles and/or cavities [16], which implies that bubble evolution is an important process in fuzz formation in tungsten.

Experiments at PISCES-B [17] and NAGDIS-II [18] suggest that the nanostructured fuzz forms when the surface temperature is between 1000 and 2000 K and incident ion energies are above 22 eV [19]. Holes form when the incident ion energy is below 22 eV at surface temperatures above 1200 K, or for temperatures above 2000 K regardless of incident ion energy [19]. This raises several questions about the mechanisms that control surface evolution, particularly with respect to the rate of gas bubble formation as a function of temperature and the mechanisms leading to severe surface roughening. Sharafat et al [20] suggested that near-surface stress gradients cause helium bubbles to be drawn to the surface in the 1000–2400 K temperature range. Krasheninnikov [21] has suggested a viscoelastic model of 'fuzz' growth, while a study by Kajita et al [19] found that 'pinholes' form on the surface before the 'fuzz' begins to grow and suggested that stresses created by bubble formation and growth were responsible for fuzz formation.

In this article, we present molecular dynamics (MD) simulations of the initial stages of helium clustering and bubble formation that demonstrate a number of sub-surface, atomic-scale processes that cause surface roughening and which are likely involved in the initial formation of nanoscale fuzz. The helium ion energy and temperature chosen for these simulations are representative of the conditions under which fuzz forms in linear plasma devices, as well as fusion tokamaks. The MD simulations show that surface roughening occurs by the addition of individual tungsten adatoms to the surface, as well as by the formation of adatom 'islands.' They also reveal that dislocation loop-punching from over-pressurized bubbles occurs not as a discrete process, but rather as a result of individual Frenkel pair production processes and the subsequent self-organization of tungsten self-interstitial atoms. The self-organization results in clustering of the individual interstitials into prismatic dislocation loops, which glide to the free surface where they form the adatom 'islands.' Furthermore, as the bubble size and/or pressure increases, the bubble can induce significant plastic deformation of the surface ligament, effectively 'bursting' and leaving a residual hole and surface pileup around the crater. While the time and length scales, and the helium implantation flux, in MD simulations are orders of magnitude different from those in the corresponding experiments; the mechanisms and physical processes of helium bubble formation, growth, and bursting that are observed should exist in any system for which conditions naturally favour the formation of over-pressurized gas bubbles, although the spatial distribution of bubble sizes will certainly vary with the implantation flux. The mechanisms presented herein are important and salient atomic-scale mechanisms, and are part of a larger multi-scale modelling effort under development to connect these phenomena with the conditions in plasma fusion devices such as ITER and DEMO.

MD simulations were performed using the Large-scale Atomic/Molecular Massively Parallel Simulator (LAMMPS) [22]. The interactions between tungsten atoms were described by the N-body potential of Finnis and Sinclair [23] as modified by Ackland and Thetford [24] and subsequently modified at short range by Juslin and Wirth [25]. Helium interactions were obtained from the tungsten–helium potential of Juslin and Wirth [25] and the helium potential of Beck [26, 27] as modified at short distances by Morishita et al [28].

The simulation cell was either a 20a × 20a × 30a or a 30a × 30a × 40a block of pure BCC tungsten, where a, the lattice parameter, is 0.3170 nm, 0.3189 nm, or 0.3209 nm at 500 K, 1200 K, or 2000 K, respectively. These temperatures were chosen to coincide with the three regimes of tungsten surface response observed experimentally [19]. Periodic boundary conditions were used in the x- and y-directions, while the z-direction was a free surface. All three dimensions of the box remained constant in length throughout the simulations. To avoid net drift of the atoms in the vertical direction due to the escape of helium atoms, the coordinate origin was relocated periodically so that the position of the centre of mass of the bottom two layers of tungsten atoms remained constant throughout the simulation.

Two different surfaces have been considered: one single crystal that terminates along a (1 0 0) plane and one containing a ∑5 grain boundary that also terminates along (1 0 0) planes. The ∑5 symmetric tilt grain boundary was introduced following the procedure of Tschopp and McDowell [29].

Helium atoms were regularly added to the tungsten bulk throughout each simulation in two different manners. For simulating helium retention, diffusion, and clustering; each helium atom had an x- and y-coordinate distributed uniformly over the entire x–y plane, while the z-coordinate was selected randomly from the depth distribution predicted for an incident helium ion energy of 60 eV by SRIM [30] with a mean depth of 1.35 nm (see the online supporting information (stacks.iop.org/NF/53/073015/mmedia)). This energy was chosen to be in the middle to high end of ion energies expected in ITER and within the energy range of ion energies that have been found to produce fuzz in laboratory plasma devices [7, 19]. To obtain statistics regarding helium retention and cluster formation, 20 independent simulations were performed at each temperature for each surface orientation, and the results were averaged.

These simulations were designed to predict the response of tungsten surfaces to the implantation of helium corresponding to an incident ion energy of 60 eV. However, to avoid the computational cost associated with (a) the dissipation of kinetic energy of an incident helium atom with too little energy to displace or sputter a tungsten atom and (b) analysing the probability of reflection, the helium atoms were inserted below the tungsten surface without kinetic energy according to the depth profile predicted from SRIM [30] calculations. We note that MD simulations at 1200 K with 60 eV helium atoms normally incident on a (1 0 0) surface (see the online supplementary data) show reflection rates on the order of 70% to 90%. Performing such simulations also requires the use of much smaller time steps until the helium dissipates its kinetic energy. Therefore, MD simulations to completely follow the helium implantation and energy dissipation processes would require that the vast majority of the simulation time be spent on helium atoms that would have little observable effect on the longer term defect—gas evolution of the solid, and complicated temperature control techniques would be required to avoid overheating within the simulation.

The total amount of helium implanted in each simulation consisted of 800 helium atoms, at a rate such that 50 000 time steps occurred between each helium atom insertion. An adaptable time step was used to account for possible high-energy interactions due to the random placement of helium while minimizing simulation wall clock time. Using an adaptable time step, each time step lasts between 10⁻⁴ ps and 10⁻³ ps, and the time between two helium insertions varies around a mean of about 10 ps. The resulting flux is about 2.5 × 10²⁷ He m⁻² s⁻¹. It should be noted that the flux expected at the ITER divertor strike point is about 10²⁴ He m⁻² s⁻¹ [31], which itself is a bit higher than the flux used in the linear plasma device experiments that observe tungsten fuzz growth [32]. Again, although there is a significant difference in incident gas flux, which can be expected to influence the quantitative details of bubble number density and spatial distribution, it is believed that MD simulations do provide practical insight into relevant mechanisms of bubble formation and growth, along with corresponding surface topology changes. Future studies will further evaluate the implications of the gas insertion rate to the gas bubble populations and surface topology changes.

A slightly different methodology was used to identify the gas bubble growth and loop-punching mechanism. In this approach, helium was introduced at an even higher rate (one every picosecond) following equilibration of the tungsten block at 1200 K. Each helium atom was introduced one at a time at the centre of a growing bubble, which was placed ten lattice units (3.19 nm) below the surface, followed by a short minimization run in which only the most recently inserted atom was allowed to move to ensure that the newly inserted helium would not overlap with pre-existing atoms. No initial kinetic energy was given to the inserted helium, but the helium atoms as a group were temperature-controlled to 1200 K using a Nosé–Hoover thermostat [33–35]. This approach allowed us to follow the growth and pressure relief process of a single bubble without artificially pre-loading the bubble with an extremely high gas pressure.

The first set of simulations evaluates the effects of temperature and surface morphology on gas atom clustering processes resulting in helium bubble formation. In each simulation, helium atoms diffuse and either form clusters or leave through the free surface following insertion. Since the mean depth for inserting a 60 eV helium atom is 1.35 nm, we expect that a significant fraction will diffuse to and be released from the surface, and that the fraction of the helium retained will be temperature dependent.

Figure 1 shows the average helium retention observed over twenty MD simulations as a function of temperature and surface/boundary orientation. For the (1 0 0) single crystal surface at 500 K, helium retention is quite high and appears to saturate around 85%. With increasing temperature, the helium retention decreases substantially, with values of 72% at 1200 K and 65% at 2000 K. For the simulations with a ∑5 grain boundary, significantly higher levels of helium retention at 1200 and 2000 K were obtained than for the single crystal specimens.

Zoom In Zoom Out Reset image size

Figure 2 presents cluster size distributions in terms of the percentage of helium in a given size cluster. The cluster sizes are divided into groups: single helium clusters, which are the most mobile; clusters containing 2–7 helium atoms, which are also mobile; 'medium' clusters containing 8–50 helium atoms; and finally, 'large' helium clusters with more than 50 atoms. The transition from small to medium clusters was selected according to cluster mobility. At a size of seven to eight helium atoms in a cluster, the pressure is sufficiently large to 'punch' a tungsten atom off its lattice site, creating a vacancy and a self-interstitial (Frenkel pair) to reduce the pressure. This 'critical size' varies with temperature due to the activation energy of self-interstitial formation and the corresponding pressure of the gas. We define this process as 'trap mutation,' and the resulting helium plus vacant lattice site cluster is essentially immobile over MD timescales. The smaller clusters, with sizes from two to about seven helium atoms, are still relatively mobile and are observed to undergo Brownian-type motion during the simulations, whereas all larger clusters that have experienced trap mutation are effectively stationary. For convenience, we have also sub-divided the medium and large helium/vacancy clusters into smaller bins based on the number of helium atoms contained in the bubble.

Zoom In Zoom Out Reset image size

The helium cluster size distributions as functions of temperature and surface/boundary orientation after insertion of 800 helium atoms, corresponding to a helium exposure of 2 × 10¹⁹ He m⁻² (not accounting for reflections), are shown in figure 2. The cluster size distribution for the single crystal tungsten (figures 2(a)) spans a broad range, with numerous small- and medium-sized clusters. At higher temperatures, the cluster distribution appears to shift to larger sizes, whereas at 500 K, there are more numerous, smaller clusters.

However, the inclusion of a ∑5 grain boundary has a pronounced effect on the cluster size distribution, shown in figure 2(b). While the helium retention is similar at 500 K to the single crystal, it is higher for the simulations with a grain boundary at elevated temperatures than comparable single crystal simulations. However, the size distribution of the clusters is substantially smaller with the grain boundary present. Indeed, the grain boundary effectively serves as a getter for the helium, presumably as a result of the added free volume provided by the grain boundary in which to trap helium. Notably, this behaviour is opposite to what would be expected from the viewpoint of grain boundaries serving as fast diffusion pathways. The higher density of helium near the grain boundary facilitates cluster nucleation, which produces a much higher density of small clusters. Thus, within the time frame of these simulations, helium cluster nucleation is much more pronounced with the grain boundary present, whereas cluster growth produces larger helium bubble sizes in the single crystals.

As previously mentioned, 'trap mutation,' or 'loop-punching,' processes of the growing helium clusters as they make the transition to bubbles produces mobile self-interstitial tungsten atoms that can migrate to the surface. Figure 3 shows the surface topology and sub-surface helium cluster distributions that are observed at a helium fluence of 2 × 10¹⁹ He m⁻² (not accounting for reflections) showing the effects of temperature, and single crystal versus grain boundary with a (1 0 0) surface orientation. The observed surface roughening shown in figure 3 indicates that individual tungsten surface adatom formation has occurred along with the formation of entire 'islands' of adatoms. The tungsten adatom and/or 'island' formation is observed at all three temperatures.

Zoom In Zoom Out Reset image size

Regions of surface islands appear regardless of surface orientation or temperature, as seen in figure 3. It is unclear the extent to which this island formation is impacted by the helium cluster size distribution or the presence of a grain boundary, since figure 3(b) also shows surface island formation, although the density of the islands appears lower than for the single crystals. However, as noted previously, the presence of the grain boundary significantly reduced the average helium cluster size.

Additional MD simulations were performed to investigate the mechanisms of island formation and to follow helium bubble growth and pressure relief to larger bubble sizes, using the modified helium insertion method described in section 2. This process, while non-physical, represents the gradual diffusion of helium atoms towards the bubble, accelerating the implantation and diffusion processes in a tractable MD simulation. Figure 4 and the associated movie (see the online supplementary material (stacks.iop.org/NF/53/073015/mmedia)) show the growing helium bubble and resulting tungsten self-interstitial atoms (SIAs) and surface adatoms at 1200 K. It is important to note that this simulation is considerably accelerated, but still provides sufficient time for pressure relief to occur through trap-mutation processes that lead to the formation and glide to the free surface of prismatic dislocation loops.

Zoom In Zoom Out Reset image size

The helium cluster, or bubble, demonstrates the same loop-punching processes described above. This first occurs when the cluster consists of seven helium atoms (figure 4(a)). The resulting SIA, along with the subsequent SIAs produced as the cluster continues to grow, is initially bound to the over-pressurized bubble as individual SIAs, and does not immediately form a SIA cluster or migrate to the free surface. Figure 4(b) shows the region around the nascent helium bubble after 22 helium insertions, in which the helium atom cluster is surrounded by five SIAs bound to the bubble.

After about 40 helium atoms have been introduced, and seven SIAs have been produced by individual trap mutations, the SIAs self-assemble into a 〈1 1 1〉 prismatic loop (figure 4(c)). As more helium atoms are added, this loop gathers more SIAs and is drawn to the surface, presumably in response to image dislocation stresses, at which point the loop detaches from the bubble and glides to the surface (figure 4(d)), where it annihilates in the form of a 28 adatom 'island' (visible in figure 4(e)). A cyclic process of repeated trap mutation, SIA trapping or binding at the bubble periphery, and ultimate self-organization into a prismatic 〈1 1 1〉 dislocation loop is observed as the bubble grows. Once a prismatic 〈1 1 1〉 dislocation loop has self-assembled from the individual SIAs at the bubble periphery, the loop continues to grow until the image force causes it to migrate to the surface, at which point it rapidly annihilates to form an adatom 'island.' This is shown in figure 4(e), in which we see a second island forms on the surface after approximately 80 helium atoms. After about 250 helium atoms, the bubble occupies roughly 100 lattice sites and another dislocation loop consisting of 53 interstitials glides to the surface, creating an adatom island which displaces a previously formed adatom island up to a second layer (figure 4(f)). It is not clear from this simulation whether the formation of the second layer was influenced by the presence of the first layer, only that such multiple-layer islands are possible.

As the bubble continues to grow, the trap-mutation and self-organization of prismatic loops continues, followed by glide to the surface of the prismatic dislocation loops, continues to occur, and as mentioned previously, this process sometimes results in multiple layers of adatom islands. In all cases, islands can be traced back along the close-packed 〈1 1 1〉 directions to the bubble. The formation and motion of the prismatic loops produces plastic deformation that, combined with the bubble growth, thins the surface ligament above the bubble. In this simulation, the surface ligament above the bubble has been reduced in thickness to only two to three layers of tungsten atoms after about 800 helium atoms are contained in the bubble. The bubble, which was initially located 3.19 nm below the surface, has now grown to approximately 1.6 nm in diameter and has created over 400 surface adatoms. The slightly aspherical bubble then ruptures, or bursts through the surface (figure 4(g)). Subsequently, helium is released from the bubble, and a cavity of roughly 450 vacancies is left below the surface with some additional surface roughness occurring as a result of the cratering and the rupture process of the tungsten surface ligament (figures 4(h) and (i)). About 450 tungsten atoms remain above the original surface after the bubble bursts.

One clear observation from this MD simulation is the manner in which bubble expansion occurs via a trap mutation or loop-punching process to relieve over-pressurization. This process occurs, as shown in figure 4 and the associated video, through many individual Frenkel pair formation events, with the resulting SIAs initially remaining trapped at the periphery of the bubble. Only upon the creation of sufficient SIA density around the bubble do these interstitials self-assemble into a prismatic loop that can glide to the surface and annihilate in reaction to image stresses, leading to the formation of adatom islands on the surface. While the helium implantation rate is tremendously high in this simulation, this rate should not change the loop-punching mechanism from that of single interstitial formation through repeated trap mutation to one involving the simultaneous creation, or punching, of a larger-sized prismatic loop. We believe, therefore, that the combination of trap-mutation processes and the self-assembly of tungsten interstitials into prismatic loops is an important part of the mechanism driving surface topology evolution. Combined with the diffusion of individual tungsten interstitials to the surface, surface diffusion, and bubble rupture processes, such mechanisms can initiate surface roughness leading to the formation of tungsten 'fuzz.' The tungsten fuzz observed experimentally forms on much longer time and length scales, most likely by extensive repetition of all of these events.

We have performed molecular dynamics simulations that demonstrate several atomic-scale, sub-surface phenomena that result in surface roughening, which is likely involved in initial formation of nanoscale 'fuzz' observed in tungsten specimens exposed to low-energy helium plasmas. Implanted helium atoms are relatively mobile in the metal and have a strong driving force to cluster. Small helium clusters remain mobile until the cluster size reaches approximately seven to eight helium atoms, at which size the mobility decreases dramatically through a trap-mutation process that created a vacant lattice site and tungsten self-interstitial (Frenkel pair). The presence of grain boundaries increases helium retention by providing a sub-surface helium sink in which the higher density of helium favours cluster nucleation that produces smaller clusters at the same helium fluence. Helium retention rates are strongly influenced by temperature, and the presence of grain boundaries.

Once the helium clusters reach a critical size, the trap-mutation process and resulting vacancy effectively immobilizes the helium cluster, while self-interstitial tungsten atoms can initially remain loosely bound to the cluster or migrate to the nearby free surface. As the density of loosely bound SIAs at a cluster increases, the interstitials organize into prismatic dislocation loops. These prismatic loops migrate along the 〈1 1 1〉 glide cylinder directions, and are drawn by image stresses to the surface, where adatom 'islands' form. As helium bubbles grow through these processes, the surface ligament immediately above the bubble eventually thins to the point of rupture. This combination of island formation, single self-interstitial migration, and bubble rupture provides a mechanism that is likely active during the initial stages of fuzz growth in tungsten.

Further work will evaluate the effects of the helium insertion rate, or flux, on bubble formation, size and number density, as well as surface topology changes; and attempt to extend these simulations to time and length scales comparable to experimental results. Additional work is also in progress to refine our understanding of the influence of temperature and surface orientation on surface topology evolution, including the rate of adatom formation, the effects of surface roughness, and the effect of temperature and bubble depth on bubble rupture. The mechanisms presented herein are important and salient atomic-scale mechanisms, and are part of a larger multi-scale model being developed to connect these phenomena with the conditions in plasma fusion devices such as ITER and DEMO.

The authors thank the other investigators in the Plasma–Surface Interactions Science Center at the Massachusetts Institute of Technology and the University of California, San Diego. They are are grateful to the United States Department of Energy for financial support through grant DE–SC00–02060.