Sputtering of the beryllium tungsten alloy Be2W by deuterium atoms: molecular dynamics simulations using machine learned forces

Material erosion and fuel retention will limit the life and the performance of thermonuclear fusion reactors. In this work, sputtering, reflection and retention processes are atomistically modeled by simulating the non-cumulative sputtering by deuterium projectiles on a beryllium–tungsten alloy surface. The forces for the molecular dynamics trajectories were machine learned from density functional theory with a neural network architecture. Our data confirms and supplements previous results for simulated sputtering rates. In the non-cumulative scenario we simulate, we did not observe reaction mechanisms leading to swift chemical sputtering. Thus, our sputtering rates at low impact energies are smaller than in comparable non-cumulative studies. The sputtering yields of the Be2W alloy are generally lower than those of pure beryllium. We found a strong dependence of the sputtering yield on the incident angle with an increase by about a factor of 3 for larger incident angles at 100 eV impact energy. In the pristine surface, a large majority of the impacting hydrogen projectiles at perpendicular impact remains in the surface.

S Supplementary material for this article is available online (Some figures may appear in colour only in the online journal) * Author to whom any correspondence should be addressed.
Original content from this work may be used under the terms of the Creative Commons Attribution 3.0 licence. Any further distribution of this work must maintain attribution to the author(s) and the title of the work, journal citation and DOI.

Introduction
Materials consisting of beryllium and tungsten have gained a lot of attention for their use in nuclear fusion reactors, especially since carbon-based materials had to be ruled out for the first wall and the divertor tiles due to their high tritium retention rates [1][2][3]. The first wall panels in ITER will use low-Z beryllium as armor material, while the divertor will be made of tungsten to withstand the high heat loads [4,5]. Beryllium that is sputtered from the first wall can therefore be transported [6] and deposited on tungsten. Indeed, beryllium-tungsten alloy formation on the surface has been found in experiments [7][8][9][10][11] and modelings [12][13][14]. We are interested in the particular niche of plasma-wall interactions that deal with sputtering and retention processes. Sputtered atoms or molecules erode the blankets and the divertor and also contaminate the plasma and lead to radiative losses. Hydrogen (H, D, T) projectiles can enter the wall material, loose their kinetic energy by successive collisions, diffuse in the material at higher temperatures or bind chemically to beryllium or tungsten at lower temperatures. These retention processes also influence the blanket stability, but the main concern is tritium activation of the wall materials which has to be limited [15,16].
The problem of D/T retention and the erosion of the wall materials by the impacts of energetic particles including H isotopes and Be has been intensively studied by experiments together with computational modeling [17][18][19][20][21][22][23][24][25][26][27][28][29][30]. Temmerman and Doerner performed experiments on D retention in codeposited W layers. Retention was enhanced with increasing incident D energy and suppressed by increasing substrate temperature or higher W deposition rates. An empirical equation was proposed to describe the influence of deposition conditions on the retention rate. Based on this equation, they suggested that for similar deposition conditions, D retention in W co-deposited layers is about two orders of magnitude lower than in Be layers [19]. In a recent work, Hakola et al studied D retention in laboratory-made Be containing samples with different composition including Be-D, Be-O-D, Be-N-D and Be-C-O-D. They reported that the inclusion of carbon by 10-15 at% can significantly increase retention (a factor of 2-10). Results from Be-D samples indicate that the fuel retention would be about 1-2 at% in co-deposits on ITER plasma facing components [28]. Björkas et al employed bond-order potential (BOP) based molecular dynamics (MD) simulations to study the sputtering process from a Be and Be 2 W surface by low-energetic (7-200 eV) D bombardments. Interestingly, BeD molecules were found as sputtering products on both surfaces, which was explained by the swift chemical sputtering mechanism [22]. Lasa et al performed MD simulations of W surfaces irradiated by Be and D projectiles with energies from 10 to 150 eV. W was sputtered by Be, but the yield was suppressed by Be deposition. The sputtering products included dimers (Be 2 and BeW) and Be-D compounds [23]. Later, Lasa et al compared the binary-collision approximation (BCA) method with MD to study the D reflection and Be erosion of Be-W mixed materials under D irradiation [29]. Sputtering from Be/W systems by H and He was also studied in reference [30]. The authors investigated a much higher energy range of mostly 500-2000 eV than in the present work with a binary collision method.
In this work, we model non-cumulative sputtering of a pristine beryllium-tungsten alloy surface by hydrogen (D) projectiles as well as retention processes. Using machine learning (ML) to develop an atomistic 'mathematic potential' has attracted more and more interest in recent years [31][32][33]. Contrary to more traditional (empirical or fitted low-body) potentials, ML potentials do not have a direct physical meaning, but rely on flexible mathematic forms with many parameters. Trained with ab initio data, they can achieve quantum chemical accuracy. On the other hand, both the generation of many-body training data and the actual training process is computationally demanding.
Based on our work in reference [27], we continued to develop neural network potential energy functions using the Behler-Parrinello approach [31] as implemented in the n2p2 code [34,35]. With neural networks providing the forces in this three-element scenario, sputtering, reflection, and retention rates for impacts at various incident angles α (0 • , 20 • , 45 • and 60 • ) on Be 2 W(001) are calculated for projectiles with 10-100 eV kinetic energy. Here, the incident angle is the angle between the projectile initial velocity and surface normal.

Computational methods
The sputtering process was simulated by classical MD trajectories using machine-learned forces that were trained on density functional theory (DFT) data. In all sputtering simulations, the projectile D atom was placed 5 Å above the target surface with random, uniformly chosen x and y-coordinates. We start with the description of the underlying DFT method.

Density functional theory
Kohn-Sham DFT was used for total energy calculations and Born-Oppenheimer ab initio MD simulations to generate and to refine training data for the neural network training [36,37]. The Perdew-Burke-Ernzerhof [38] functional was used to describe the electronic interactions. All DFT calculations were performed with the Vienna ab initio simulation package (VASP) [39,40]. Projector augmented wave (PAW) [41] potentials for W (with six valence electrons), Be (with two valence electrons) from the VASP library replaced the inner electrons. A Gamma-centered k-point mesh of 3 × 3 × 3 was employed. The Kohn-Sham orbitals were expanded in a plane wave basis set with periodic boundary conditions. The cut-off energy was selected to be 550 eV [42]. As described below in detail, the initial set of training structures, energies and forces were obtained from ab initio MD sputtering trajectories on rather small supercells. In these trajectories, the lowest layer of atoms in the slab was fixed while all other atoms were allowed to relax. A convergence criterion of 10 −4 eV Å −1 was used for the forces.
The steep repulsive part of the potential energy is well described by the density functional with PAW even for small interatomic distances. Figure S1 (http://stacks.iop.org / NF/61/016031/mmedia) shows a comparison between DFT with PAW and an all-electron calculation using the CCSD(T) method with the cc-pVqZ basis set.

Neural network potential
The neural network potentials (NNPs) were developed according to the Behler-Parrinello model using the n2p2 code recently developed by Singraber et al [34,35]. An asset of n2p2 for constructing the NNP is that both total energy and atomic forces of reference structures are used for optimizing the weights and biases in the fitting process. The fitting algorithm is a multi-stream Kalman filter [35]. After the NNP is created, the MD simulations were achieved by an interface to the LAMMPS code [43].
In the Behler-Parrinello approach [31] the total energy E pot of one configuration is obtained by summing over 'artificial' atomic energies E i , which are the output of an atom-centered element-specific neural network. Atomic coordinates are converted to symmetry-invariant atom-centered symmetry functions which describe the local chemical environment of each atom. These symmetry functions are used as input nodes in the neural network. We have used a feedforward neural network topography with three hidden layers of 20 or 30 nodes each. The soft-plus activation function and a cut-off radius of 7 Å were employed [27]. The input nodes consist of radial and angular Behler-type symmetry functions [44], which are specified for the elements H, Be and W in tables S1 and S2 in the supplementary information.

Training data generation and refinement
By training our NN potential with configurations from direct DFT sputtering simulations on small systems, we have the advantage that sputtering within the relevant energy ranges and from relevant configurations is 'by definition' included, provided that the training data is sufficiently large. The initial training data was generated by performing ab initio MD simulations on the slab A of a Be 2 W(001) surface [45] with 48 atoms. The slab geometry was first optimized, and then equilibrated for 2 ps at 300 K within the canonical ensemble using the Nosé-Hoover algorithm [46,47]. A further equilibration time of 0.02 ps was subsequently used to generate different initial surfaces for the sputtering processes. Then, impacts of single D atoms with an energy of 50 and 100 eV at incident angles of 0 • and 45 • were simulated. We used a much smaller integration step of 0.1 fs than in reference [27] in accordance with the mass ratio between D and 9 Be of about 4.5 and the corresponding higher velocity of D atoms with the same kinetic energy.
The initial training data consisted of 20 210 configurations extracted from 160 ab initio MD trajectories. We have implemented the iterative refinement protocol as shown in figure 1 of reference [27] to successively improve the training data. At the beginning, two preliminary NNPs were fitted to the initial training data generated from ab initio MD as described above. These two NNPs have the same topography but differ in the starting parameters of the weights and biases. One of the fitted NNPs is then employed in short MD simulations (20 fs) of the sputtering on a Be 2 W surface with 48 atoms by D projectiles at energies from 10 eV to 100 eV with random incident angles between 0 and 45 • . Next, the energies and forces of these new configurations extracted from MD trajectories are predicted by the NNP not used for the simulation. The configurations with large differences of energies or forces from different NNPs were added back to the training set after running static DFT calculations on them.
After five iterative refinement steps, the final reference data set consists of 27 502 configurations containing 32 Be atoms, 16 W atoms and one D atom each. The energy distribution of the configurations in the reference data is given in figure  S2. 24 727 total energies and 3634 869 atomic forces are used to train the final NNP and the remaining data (∼10%) is part of the test set which is used to validate the potential and to avoid overfitting. The computational effort for training this NNP is detailed in table S3. With 50 iterative training epochs, the RMSE in the test set converged to 3.3 meV/atom for energies and 0.23 eV Å −1 for atomic forces, slightly larger than the corresponding values in the training set (2.7 meV/atom for energies and 0.18 eV Å −1 for atomic forces). Overfitting was not observed during the training process.

Non-cumulative sputtering simulations of a Be 2 W surface by D impacts
In reference [27] it was shown that an NNP fitted to relatively small systems could be applied for reactive sputtering simulations of a much larger system. It was also found that the sputtering yields from Be impacts at 100 eV on a Be(0001) surface consisting of 490 atoms were very close to the ones of a Be(0001) surface with 2000 atoms [27]. Here, for saving CPU time, the target surface super cell was created by relaxing a Be 2 W(001) surface, consisting of eight Be and eight W atomic layers (600 atoms) with a size of 22.3 × 22.3 × 27.7 Å 3 . After relaxation, the structure was equilibrated for 2 ps at 300 K in the NVT ensemble using a Nosé-Hoover thermostat [46,47]. Then the pristine Be 2 W(001) surface was irradiated by D impacts with kinetic energies of 10, 20, 35, 50, 75 and 100 eV. For each energy, four sets of simulations with incident angles of 0 • , 20 • , 45 • and 60 • were performed. The second polar angle was chosen randomly. 2000 or 5000 separate MD trajectories for 120 fs were generated with an integration time step of 0.1 fs. In those intermediate cases, where it could not be decided within 120 fs if sputtering has been taking place, further 120 fs integration time were added to the trajectories.
At all incident energies, the projectile has a high initial velocity. A D atom with a kinetic energy of 100 eV has a velocity of ∼1.0 Å fs −1 . It can happen that some projectiles pass straight through the atomic layers in the super cell, thus adding to the D retention budget.

Results and discussion
The   DFT both for energies per atom and for atomic forces are very close except for a few configurations, where (partial) melting occurred after impact in the smaller supercells used for training data generation. Even though equilibrium structures were not explicitly included in the training data, table 1 shows that the equilibrium lattice structure is well described by the NNP potential energy function. Figure S3 shows the Be-H dimer potential energy curve obtained from DFT and NNP. The agreement is very good considering that we have not explicitly included any dimer structures in the training data.
The sputtering yield (sputtering of Be atoms per number of total impacts), reflection rate (reflection of D atoms per total impacts) and retention rate (retention of D atoms per total impacts) depending on incident angles and energies are summarized in table 2. The beryllium sputtering yields are also shown in figure 2. The statistical error of the sputtering yield was estimated as Δ = σ/ √ N, where σ and N are the standard deviation from assuming a Bernoulli distribution and the number of MD trajectories, respectively.
For normal impacts, the sputtering yields are compared with MD simulations on a BOP force field from Lasa et al [29]. At 50 and 100 eV, the NNP based sputtering yield is about a factor of 2 lower than the BOP value. At very low incident energy (10 eV), the difference in sputtering yield is about a factor of 10. In contrast to Lasa et al [29], we did not encounter BeD chemical sputtering events at 10 eV. However, the low probability of BeD formation does not allow for any conclusions at this point.
We can compare the sputtering yield of the alloy surface with a pristine Be(0001) surface for perpendicular D impacts with the same incident energy. At incident energies from 50 to 100 eV, the sputtering yields on the Be(0001) surface reported by Björkas et al [49] are about 0.02 to 0.03, while for the Be 2 W(001) surface the yield is lower and slightly increases from ∼0.01 to ∼0.02, indicating that the tungsten admixture enhances the stability of the Be surface. From 35 eV to 50 eV the sputtering yields increase rapidly. At the lowest incident energy (10 eV) a near-zero sputtering yield (0.0004) is obtained for this surface. Its sputtering threshold energy E th defined here by no sputtering yield occurring in 5000 trajectories can be assumed to be slightly below 10 eV. Figure 2 also reveals a strong angular dependence of the sputtering yield for this surface. At high incident energies of 75 and 100 eV, the sputtering yield increases by more than threefold with increasing incident angle from 0 • to 60 • . An increase by a factor of two is found for 35 and 50 eV. Interestingly for 20 eV, the sputtering yield for all incident angles is very close to 0.001 except for 45 • with 0.003. We note that the error bars are already quite large for these rare sputtering events at 20 eV. It is worth noting that no sputtering events have been observed at an incident energy of 10 eV for non-perpendicular impacts. The Nordlund group studied the incident irradiation angle and surface orientation dependent sputtering of W by argon and helium projectiles with MD simulations [50]. Their observed angular dependence is not as high as in our case. Sputtering increases from normal incident to a maximum for an incident angle of ∼30-40 • , then decreases to zero for grazing angles. For example for 100 eV Ar, the sputtering yield of a W(001) surface is 0.18 for 0 • and 0.32 for 30 • [50].
The D reflection rates as function of the incident energies are given for the various incoming angles in figure 3. Again, the results of MD simulations from Lasa et al [29] are included for comparison. For normal impacts at 10 eV, the NNP based D reflection rate is 44% and lower than the BOP value of ∼57% [29]. At higher incident energies (20-100 eV), the NNP based reflection rate decreases linearly from 51% to 23%. The BOP data has only two points for 50 and 100 eV where the rate is about ∼30% for both [29]. A similar trend occurs for grazing angles. From low to high incident energies the reflection rate increases first, approaches a maximum at around 20 eV, and then goes down until 100 eV. The decrease in reflection yield on the high energy side can be attributed to the increase in kinetic energy that is available to deform the structure and make space for depositing the projectile on or below  [29] for a Be 2 W surface by non-cumulative D impacts at normal incidence are included for comparison (asterisk).
the surface. In MD simulations of sputtering from a berylliumoxide surface by D impact [51], the reflection probability even decreased by one order of magnitude between 10 and 200 eV. On the lower energy side, this decrease in reflection yield may be attributed to a chemical sticking effect, where the D-Be adatom binding energy starts to play a role. As can be expected, larger incident angles lead to higher reflection rates. At 100 eV, the reflection rate for 60 • is nearly three times that for 0 • .  We can also analyze the distribution of the angles of the sputtered Be atoms in dependence of incident angles and energies. In figure 4 the respective histograms are shown for 100 eV kinetic energy of D impact. One sees that Be is sputtered in directions between 0 and 60 • with a preference for small angles and with only little correlation to the incoming angle ( figure S4).
In the same way, we looked at the distribution of the angles of the reflected D atoms. In figures 5 and S5 the respective histograms are shown for four different incident angles at 20 eV and 100 eV. The distributions are broad, and D is reflected in all directions with little correlation to the incoming angle ( figure S4). One might expect that large incident angles cause large angles of the reflected D but interestingly this picture of mostly elastic reflections does not hold as, somewhat surprisingly, larger reflection angles play a minor role in the whole energy range between 20 eV and 100 eV.
One must keep in mind that in figures 4 and 5, for geometric reasons, an angle near zero (perpendicular exit channel) has a lower probability than a larger one. This stems from the unequal length of the latitude circles in a polar plot. One can renormalize the distributions so that a random angle gives a uniform distribution by scaling them by 1/sin(θ). The result is shown in figures S6 and S7 of the supplemental information which correspond to figures 4 and 5. It can be seen clearly that a perpendicular exit (∼0 • ) of the sputtered Be atom is not unlikely due to physical reasons as figures 4 and 5 seem to indicate. Indeed, for impacts with 0 • and 20 • they are preferred. On the other side, grazing angles (∼70-90 • ) are found much less than would be expected by random scattering.
We know of no angle-dependent modeling of our system, however, in a recent work [50], the dependence of the angle of incoming Ar atoms on the sputtering yield of W was studied. Two of the energies employed there, 85 eV and 100 eV are within the range of the present study. As expected, like in our alloy, at these energies no sputtering of W but only reflection occurs. The observed reflection yield was largely independent of the angle of the incoming Ar atom. At higher energies, 200 eV and 300 eV, sputtering of W is believed to set in [50]. Depending on the W surface, their sputtering yields differ but are not unlike ours with a maximum at 30-60 • and low sputtering yields at ∼90 • [50].

Conclusions
In this work, an NNP was developed for a surface of the Be-W-H system. This NNP was applied to non-cumulative sputtering simulations of the Be 2 W(001) surface by D impacts with kinetic energy ranging from 10 to 100 eV at various incident angles. At incident energies from 50 to 100 eV, the sputtering yields for the Be 2 W(001) surface slightly increase from ∼0.01 to ∼0.02. The yield of the Be 2 W(001) surface is overall lower than that of a Be(0001) surface, which means that the tungsten admixture enhances the stability of the Be surface. The sputtering threshold of a Be 2 W(001) surface irradiated by D impacts lies slightly below 10 eV. Both the sputtering yield and D reflection probability of the Be 2 W(001) surface reveals a strong dependence on incident angles. For 100 eV impact energy, there is about a factor of 3 stronger sputtering for 60 • impact angle than for perpendicular impact. Last, the dependence of the angle distribution of the sputtered Be/reflected D atoms on incident energy and angles were discussed. We found that Be is sputtered in directions between 0 and 60 • and D is reflected in all directions, and both have only little correlation to the incident angle of the projectile.