New helium bubble growth mode at a symmetric grain-boundary in tungsten: accelerated molecular dynamics study

ABSTRACT This work, with an emphasis on helium irradiation rates appropriate for fusion-plasma conditions, advances the understanding of helium evolution at grain boundaries in W, an important consideration in the understanding of W as a plasma-facing component. Using accelerated molecular dynamics, helium bubble nucleation and growth at a symmetric Σ5[100](310) tilt grain-boundary in W is studied. The simulations reveal that the growth mode associated with bubble growth at the grain-boundary leads to a suppression of the helium supply to the bubble and hence to arrested growth. Such an unconventional bubble growth mode may dominate in materials with a high density of sinks. GRAPHICAL ABSTRACT IMPACT STATEMENT Simulations of helium bubble nucleation and growth at a tungsten grain-boundary approaching helium-irradiation rates for fusion first-walls reveal a novel growth mode.

The design of materials that withstand the extreme radiative environment present in fusion reactors is an outstanding challenge. For example, tungsten (W) is a leading material choice for plasma-facing material (PFM), due to its excellent material properties including high melting temperature and low sputtering yield [1]. As a divertor component, however, W will be implanted with a high flux of low energy ( ∼ 50 eV [2,3]) helium (He) ions. This He irradiation has been shown to cause morphological changes at the nanoscale, such as the formation of low-density 'fuzz' or 'coral' at surfaces [1]. This is a key problem since these morphology changes can severely disrupt the material's performance, ultimately leading to contamination of the plasma. As a result, many modeling and simulation studies at the atomistic scale have been undertaken to understand the various aspects of He behavior in W, including W surface evolution [4], the mobility of small He clusters [5][6][7][8], He bubble morphology in the bulk of W or near W surfaces [9][10][11][12], He segregation at surfaces and grain-boundaries (GBs) [5,13,14], and the dynamics of small He clusters near GBs [15]. When used as a PFM, W is typically polycrystalline. Experimentally, ultrafine-grained and nano-grained W have been actively studied to understand the role of grain boundaries under high-flux, low-energy He irradiation [16][17][18]. On the modeling and simulation side, it has been shown that GBs are generally strong sinks for He atoms or even small He clusters [5,13,15]. Once a He atom segregates to a GB, its migration is greatly retarded [13,19]. However, He bubble growth at GBs in W has not received significant attention in the literature (although several studies of He bubbles in bcc Fe have been carried out [20][21][22][23][24][25][26]), which suggests that this important aspect of the effect of GBs in the evolution of He requires further study. This work, using accelerated molecular dynamics (AMD), begins to address this issue, with an emphasis on He irradiation rates appropriate for fusion-plasma conditions, thereby advancing our understanding of He evolution at GBs.
The AMD method used in this work is the parallel replica dynamics (ParRep) method [27,28] as implemented in the LAMMPS simulation package [29]. Par-Rep addresses the timescale limitation of conventional molecular dynamics (MD) for systems with rare events by using multiple replicas of the system to achieve a parallel speedup in the time domain [28]. By deploying ParRep on the Titan high-performance computing (HPC) clusters at Oak Ridge National Laboratory, the simulation timescales reported here can extend into the μs. The W grain-boundary chosen in the current study is a wellstudied symmetric 5[100](310) tilt GB [30][31][32]. This GB is chosen for its well-ordered GB structure at the temperature of interest, i.e. 1000 K, typical of a fusion environment. The well-ordered GB structure helps maintain the parallel efficiency of ParRep, which is typically on the order of 80% for the simulations presented here. To describe the interatomic interactions, we employed the He-W, He-He, and W-W interatomic potential set developed by Juslin and Wirth [33]. In Ref. [33], the original W-W interaction derived by Ackland and Thetford [34] was modified at distances smaller than 0.26 nm as the original potential did not describe short-range repulsion accurately. This potential is well tested and has been used successfully before. Finally, in all simulations, events involving the motion of He in He clusters containing four or more He atoms were ignored. In these clusters, the He is effectively in a gaseous state and the rate of motion is very high. We ignore this motion to increase the efficiency of the simulation. The validity of this choice is discussed near the end of the paper.
The ground state structure of the 5[001](310) tilt GB in W is shown in Figure 1(a-c). The simulation supercell in the ParRep simulations is 6.6 nm in x = [310] grain1 /   is in the y direction. Periodic boundary conditions are applied in the y and z directions. In the x direction, which is perpendicular to the GB plane, slab surfaces with a vacuum thickness of 0.8 nm on each side of the two surfaces are imposed (Figure 1(c)). The relatively large size of the supercell in the GB plane alleviates possible spurious image interactions due to the presence of the large defect clusters that are generated during the ParRep simulations. Figure 1(a) shows the typical 'kitelike' structural units in the GB region (viewing down the tilt axis). In Figure 1(b), a side view of the supercell (along [130] grain1 /[310] grain2 , the z direction) is shown. Note that the two grains shift with respect to one another along the [001] direction, therefore breaking the mirror symmetry across the GB plane. This structure is in excellent agreement with earlier results in W using firstprinciples density functional theory (DFT) calculations [30,32]. Such an asymmetric atomic structure for the tilt  [001](310) GB has also been observed experimentally by Mo and Ta [31]. This structure thus provides a realistic starting structure for our ParRep simulations.
The simulations begin by placing 8 He atoms inside a pre-existing W vacancy, at the center of the GB plane. These initial He atoms are introduced into a spherical region 0.14 nm in radius, at randomly generated coordinates [9]. This serves as the nucleus of a He bubble. The simulations were then carried out at a fusion relevant temperature of 1000 K, maintained using a Langevin thermostat. The supercells were scaled to match the predicted lattice constant for W at 1000 K of 0.3183 nm. A time-step of 0.5 fs was used in the ParRep simulations to account for the small atomic mass of He. After an initial thermal equilibration time of 10 ps, the structure was then used as the starting point of the ParRep runs (with typically 40-100 replicas, each using 32-128 CPU processors to parallelize the force call) with a He insertion rate of 1 He atom/10 ns. This insertion rate is close to the upper bound of what can be expected at a nearsurface GB in ITER [35]. 1 The slow growth rate used here provides ample time for defects generated during the simulations to migrate. Helium is inserted directly into the center of the existing bubble, with care taken to ensure atoms are not placed too close to one another (with the condition that no current atoms in the simulation are within a distance of 0.06 nm of the added helium atom). While this simulation procedure was used in our previous work examining He bubble growth in the bulk [9], it only approximates the real He arrival through a diffusive process. We will comment on the appropriateness of this condition below.
In Figure 2(a-f), the temporal evolution of the He bubble growth process revealed by the ParRep simulations is shown. The He bubble exhibits a first 'trap mutation' (TM) event (in which a W atom is pushed out of the surface of the bubble, creating an additional vacancy) at 75 ps. Following this process, W atoms are pushed out of the bubble surface into interstitial positions, leaving behind vacancies that allow the bubble to grow. In this case, the W interstitials are created within the GB plane, initially in the center of the kite structure (see Figure 1(a)). Continued TM events create more isolated W interstitials and some of these W interstitials escape away from the (at that time) small He bubble. One snapshot of this early growth stage is captured in Figure 2(a), after 108 ns in the ParRep simulation, where a migrating isolated W interstitial residing far from the He bubble is highlighted by the circle. The overall process shares qualitative similarities (such as TM events during growth induced by the high He pressure) with the growth of He bubbles in bulk W [9], but key differences are also observed. Due to the strong tendency for interstitials to remain trapped in the GB plane, the self-interstitial-atom (SIA) loops observed in bulk W do not form. Rather, at the GB, the W interstitials gradually start to aggregate around the He bubble, forming a very ordered structure. Such aggregation starts at one side of the He bubble, as shown in Figure 2(a), and gradually expands to the other side of the bubble in Figure 2(b) (at 363 ns in the ParRep simulation).
As the simulation continues, the two branches of the W interstitial aggregate connect with each other, as shown in Figure 2(c) (at 850 ns in the ParRep simulation). Gradually, the W interstitials develop into a halo around the He bubble, surrounding it with a full ring of excess W (see Figure 2(d,e), at 1254 ns and at 1821 ns of the ParRep simulation, respectively). We have identified that the edge lines of the W interstitial halo regions are in the , [311], or [001] directions within the GB plane (see Figure 2(e)) [36]. These also correspond to the high atomic density directions at the GB plane.
In Figure 2(f), a side view of the structure presented in Figure 2(e) is shown (viewed down the tilt axis of the GB). In Figure 2(f), there is an asymmetry in the strains around the bubble, which is possibly due to the stochastic nature of the growth process. It is important to emphasize that, in contrast to the bulk, in which loops eventually break free from the bubble and can thus reach some far away sink (such as a surface), this does not happen at the GB. Rather, the interstitials remain trapped at the GB in the vicinity of the bubble. Also, bubbles growing at the GB are seen to (at least initially) grow as 2D platelets, not as the more-or-less spherical objects observed in the bulk. This indicates that creating new W Frenkel pairs is less energetically costly in the GB than in the bulk. As the bubble grows, however, the high interfacial energy cost of platelet geometries could presumably eventually lead to a transition to 3D growth. How this transition will affect the growth mode of bubbles will be the subject of a future investigation.
As mentioned above, for the simulation presented in Figure 2, He was introduced directly into the bubble. While this might be a reasonable approximation of the initial arrival process, as the interstitial halo grows, the manner in which He arrives at the bubble might be perturbed. In order to understand the effect of the W interstitial halo on He arrival, a new ParRep simulation was performed. Beginning with the final structure from the previous simulation (specifically, the structure shown in Figure 2(e)), we introduced a new He atom in the GB plane away from the He bubble. This initial configuration is shown in Figure 3(a). The additional He atom is quickly attracted to the periphery of the W interstitial halo, as shown in Figure 3(b). The newly placed He atom is able to break away from the W interstitial halo region, migrating randomly in the nearby region, before it is attracted back to the W interstitial halo. After about 0.47 μs of ParRep simulation time (Figure 3(c)), the extra He atom comes to rest at a nearby edge or step in the W interstitial halo. Thus, while the extra He atom is attracted to the halo, the halo prevents it from reaching the bubble. This second simulation ( Figure 3) suggests that the interstitial halo structure found from the first ParRep simulation (Figure 2) may partially hinder subsequent He atom arrival into the existing He bubble, thereby possibly limiting its size, or at least its growth rate. That is, if, as the bubble grows, the trapped W interstitials impede the arrival of subsequent He, the bubble will stop growing once the halo is complete.
To provide a more quantitative understanding to this effect, we designed a special configuration that captures the salient features of the interstitially-loaded GB. This configuration has 1/4 monolayer of extra W (or 'interstitial') atoms at the GB plane. The extra atoms at the GB plane were placed according to the observed positions of the trap-mutated W interstitials in the first ParRep simulation. As shown in Figure 4, the extra 1/4 monolayer of interstitial atoms was introduced in a square patch at the GB. After the configuration was thermally equilibrated at  300 K for 1 ps and relaxed, climbing-image nudged elastic band (CI-NEB) [37,38] calculations were carried out to calculate migration barriers for He in various regions of the GB plane, corresponding to different regions around the interstitial halo of Figure 2.
Three migration paths are shown in Figure 4. The migration path for He in the GB plane far away from the W interstitials is the A2 to A3 path in Figure 4. This is equivalent to He migration in the pristine GB. The migration barrier for this path is 1.19 eV, an energy barrier much larger than in bulk W (0.15 eV), indicating much slower migration of He in the GB than in the bulk, consistent with previous work [13,19]. He migration from the pristine region to the reconstructed halo zone can occur via various processes in multiple directions. For example, the migration path from B1 to A1 (see Figure 4) is a path from the reconstructed region to the pristine structure along a [310] direction. The migration barrier is 3.25 eV while the reverse migration barrier is 4.16 eV. The latter value represents the barrier for a single He atom to penetrate the halo along the [310] direction. Obviously, this is a very high energy barrier to overcome, even at 1000 K. An alternative path, along a [001] direction, corresponds to the A4 to B3 via B2 path (see Figure 4). This path has an effective migration path of 2.39 eV, passing through the shallow intermediate state B2. Again, this is a very high barrier at 1000 K. These results indicate that the interstitial halos that surround the He bubble can block the arrival of subsequent He atoms, which implies an unconventional constrained He bubble growth mode at the GB compared to that encountered in bulk W. The barrier for He migration from state B2 to B3, ∼ 1 eV, represents He migration within the halo structure itself. This suggests that, if He can penetrate the halo or is directly 'deposited' into the halo from the bulk, migration will occur faster than in the pristine boundary.
Finally, to ensure that, even in the pristine boundary where the migration energy of He is 1.19 eV, He would have sufficient mobility to nucleate and grow bubbles (as opposed to becoming trapped where it first 'sticks' to the boundary), we also simulated the He atom bubble nucleation process itself. In this scenario, He atoms are introduced randomly into the GB plane (that is, at a random position within the GB plane) with an insertion rate of ∼ 1 He atom/100 ns, which is in the range of He irradiation rates under plasma fusion conditions. The insertion of He atoms at the GB plane is justified by the fact that, once a He atom is introduced into bulk W, the He atom migrates quickly to the GB plane, with a segregation energy of 2.3 eV (the energy gained by placing a He atom at the GB as compared to placing it in the bulk). Such a high segregation energy practically eliminates any out-of-plane emission of He atoms once they have arrived at the GB, at least on the timescale of bubble growth.
During this last ParRep simulation, when two He atoms encounter each other and form a pair, the pair is mobile and can also dissociate back into isolated He atoms. However, the He dimer can also undergo selftrapping through a TM event, at which point it becomes immobile and serves as the seed for further growth of a bubble. This self-trapping TM event occurs at lower interstitial He content as compared to the bulk. In Figure 5(a), by 719 ns when a total of 8 He atoms had been introduced into the simulation cell, a cluster comprised of 5 He atoms and 3 W vacancies had formed through TM processes. By 806 ns, another cluster containing 3 He atoms and 1 W vacancy formed ( Figure 5(b)). At 933 ns, a third cluster, this time containing 2 He atoms and 1 W vacancy, formed ( Figure 5(c)). Two of the clusters grew, such that, by a time of 1166 ns ( Figure 5(d)), clusters containing 4 and 6 He atoms had been formed. These simulations indicate that, even though He migration at the GB is slower than in the bulk, there is sufficient mobility at the GB for He atoms to migrate and nucleate and grow bubbles.
In order to improve the parallel efficiency of the Par-Rep simulation, the dynamics of any clusters containing four or more He atoms were ignored in the ParRep transition event detection. This assumes that once a cluster of four He atoms forms and self-traps, it remains selftrapped and does not dissolve back into the W matrix. To test the validity of this assumption, we have calculated the energetics of 1 He atom joining with a 3-atom He cluster occupying a W vacancy at the GB plane to form a 4-atom He cluster. We found that the dissociation barrier for the cluster to break up or emit a He atom is high, about 2 eV, indicating that the assumption of the stability of these clusters is valid and that ignoring He events for such clusters will not corrupt the dynamics of the system.
Together, the ParRep simulations described above provide a picture of He bubble growth at GBs that is significantly different than that found in the bulk. In both the bulk and at the GB, a He atom can migrate, either three-or two-dimensionally, respectively, encountering other He atoms and nucleating He bubbles ( Figure 5). In both cases, as a bubble grows and the pressure increases, it emits interstitials, typically forming some kind of dislocation loop structure. In the bulk, these loops can escape the bubble. However, at the GB, they are trapped and decorate the bubble (Figure 2). Because these interstitials surround the bubble in the plane of the GB, the same plane in which the He is diffusing, once the bubble is completely surrounded, its supply of He is effectively cut off. This is a consequence of the high energy barrier for He to cross through the interstitial halo around the bubble ( Figure 4). Further, as He encounters the halo, it is likely to become trapped at features along the halo, such as steps (Figure 3), which then likely becomes the nucleation site of a new bubble. Thus, at the GB, we expect a more dispersed set of smaller bubbles, compared to equivalent conditions in the bulk, at least over some time regime. However, once the GB plane is decorated by bubbles that are separated by regions of interstitial-loaded GB (the halos), new He atoms arriving from the bulk can join the GB in the reconstructed region, where they will be able to diffuse towards the bubble and reinitiate the bubble growth. We thus expect that the bubble growth dynamics will initially rise relatively quickly, temporarily saturate, and then resume growth again.
Similar behavior can occur in the bulk [39]. However, in that case, the mechanism is different. In the bulk, if the He arrival rate is high enough, He can pin on dislocations that were punched out from the bubble but did not yet escape. This leads to the nucleation of new bubbles. At the GB, the timescale for the emitted interstitials to escape seems to be extremely high, and thus would not be nearly as sensitive to the He arrival rate.
We have only examined He bubble growth at one special grain boundary. We have seen that the structural evolution of the bubble is related to the kite structure of the boundary and the ability of the boundary to reconstruct as interstitials are emitted from the bubble. Such accommodation of interstitials has been observed in previous work and shown to lead to significant changes in GB properties [40][41][42]. Thus, the extent to which this behavior is general would depend on the ability of the GB to accommodate interstitials in a very ordered and compact way. If interstitials would rather disperse, as they might at other types of GBs [43], the behavior might differ from what was observed here. Further work is needed to elucidate the relationship between the growth kinetics and the boundary structure.
That said, the bubble growth behavior observed here has interesting consequences for the evolution of plasmaexposed tungsten. As bubbles grow at GBs, they may not lead to as significant changes in surface morphology as bubbles in the bulk, since the interstitials cannot escape to the surface as readily. Thus, the introduction of a high density of sinks for He may help delay, if not mitigate, the surface evolution responsible for the formation of fuzz. This is similar in spirit to the idea that a high density of sinks can mitigate bubble growth in materials such as oxide dispersion strengthened (ODS) steels [44], though the origin of the effect is different here in that, instead of dispersing He in small bubbles everywhere, the sinks instead trap the interstitials and keep them from escaping to the surface.
To conclude, using ParRep simulations, we have found that the kinetics associated with bubble growth at GBs in W is qualitatively different from that in the bulk. Interstitial loops, though punched out from the bubble just as in the bulk, are not able to escape due to the thermodynamic interactions with the boundary itself. This leads to an arrested growth in which the He supply to the bubble is effectively removed, at least transiently. These results reveal an unconventional mode for bubble growth that may dominate in materials with a high density of sinks.

Note
1. Under a typical ITER helium flux of 0.7-3 × 10 24 He/m 2 s, as reported in Ref. [35], and assuming that only a fraction f = 0.15 of the helium atoms reach the GB, for the simulation supercells used in this study with a GB area of 0.3-1.2 × 10 −16 m 2 , the helium insertion rate to the GB is equivalent to 1 He/(17-300 ns).