Modelling the impact of argon atoms on a tungsten surface

Sputtering from plasma-facing surfaces upon particle impact is an important process in material science. It is especially relevant in the diverter region of fusion devices, which nearly always consist of tungsten. Besides the main plasma components, argon is used in fusion devices to improve energy confinement. As a consequence, hot Ar atoms interact with W surfaces and can cause sputtering and other material degrading events. Atomistic simulations of the plasma-wall interactions make it possible to carry out a detailed analysis of sputtering, reflection, and retention processes. We report the results of molecular dynamics simulations with neural network potential energy expressions modelling the bombardment of tungsten samples by argon atoms in the energy range from 100 to 800 eV. The obtained sputtering results are in good agreement with available literature data. Furthermore, our data provide additional insight into atomic details of the processes involved in sputtering. We also investigate the effect of surface temperature on sputtering and reflection probabilities, which significantly affects the irradiation process at higher impact energies.


Introduction
Erosion of plasma-facing materials (PFMs) occurs under the impact of energetic plasma particles. The thickness and stability of the PFM is reduced and, equally important, the plasma is contaminated with particles that lead to radiative losses [1][2][3][4]. Tungsten as a plasma-facing material has many useful properties, the low erosion yield, high melting point [5,6], and favourable activation behaviour [7] being the most important ones. In particular, the ITER divertor will consist of tungsten. It will be in contact with fluxes of hydrogen isotopes (H/D/T) and inert gases such as helium (He), argon (Ar), neon (Ne) and also nitrogen (N_2). The energy of ion bombardment in the divertor region is relatively low, usually ranging from a few eV to a few keV [4]. Seeding of inert gases is used to diminish the plasma energy by radiation and dissipation over a large surface area [8]. Consequently, the interaction of these gases with tungsten needs to be known, in particular the sputtering yield as a function of environmental parameters. In addition, inert gases affect the retention of hydrogen isotopes in PFMs.

Supplementary Information
The online version contains supplementary material available at https://doi. org/10.1140/epjd/s10053-022-00495-3. a e-mail: michael.probst@uibk.ac.at (corresponding author) Decrease or increase in retention due to impurities is also determined by other environmental parameters, such as tungsten sample temperature and surface structure [8]. The surface temperature of tungsten in ITER should be less than the recrystallization temperature of W (1300-1600 K [9]). However, this temperature could be exceeded in fusion operations due to various plasma instabilities [10]. Tungsten surfaces will also be exposed to fluxes from Be-main chamber erosion. This can cause the formation of various Be-W alloys, such as Be_22W, Be_12W, and Be_2W with melting temperatures of 1300, 1750 and 2250 K respectively.
The sputtering of tungsten by Ar has been measured in several experimental works [11][12][13][14]. Stuart and Wehner [11] studied the sputtering yields of polycrystalline targets at low impact energies. The lowest impact energy of Ar that causes sputtering is 35 eV with an extremely low (and therefore inaccurate) yield of 10-5 atoms/ion. In this work, besides Ar, bombardments by Ne, Kr, and Xe gases were measured. At lower impact energies, the differences between impact gases are insignificant. These very low yields in the 20-40 eV region increase quickly with increasing kinetic energy. In [12], sputtering yields of 28 elements were measured as a function of the incident energy by bombarding polycrystalline spherical targets with a diameter of about 6 mm with Ar ions. The bombardment of smooth tungsten surfaces and fuzzy ones with Ar was studied by Nishijima et al. [13] experimentally. The yields for the smooth surfaces were close to the theoretically predicted ones, while those for the fuzzy surfaces were significantly smaller. Golczewski et al. [14] measured the sputtering from polycrystalline nanolayers deposited onto cylindrical quartz crystals. Since different techniques and samples led to the same yields, this investigation seems particularly robust. Parts of these results are also presented in the discussions below. The influence of noble gases on D retention was investigated in other studies [8,15,16]. In Ref. [8], the effect of plasma impurities on retention of D was reported for two different temperatures of the tungsten samples. The results show that Ar impurities can create escape channels for hydrogen isotopes and can remove damaged surface layers. It was also found that at higher temperatures less D is retained, implying that knowledge of the temperature dependence of these processes is important.
Computational methods have taken an important place in materials science in general and for fusion material research in particular. To obtain atomic scale information about PFMs under irradiation, computational methods such as Monte Carlo (mostly based on the binary collision approximation, BCA) and molecular dynamics (MD) are used. BCA simulations very efficiently provide information at the microscopic level [17][18][19] but do not cover chemical effects [20]. Chemical bonding and reaction processes are especially important at lower impact energies [21]. MD simulations based on accurate interaction potentials, on the other hand, are restricted due to computational requirements but access a wide range of often complementary data. For instance, Marenkov et al. [22] studied tungsten sputtering yields under low energetic Ar atoms by MD simulations and reported angular and energy distribution of sputtered particles. The energy spectrum in this case, as in work [23], is well described by Thompson's formula. This formula enables an estimation of the surface binding energy (SBE), an important input quantity for BCA simulations. Ref. [23] showed that SBE can be extracted from MD simulations. The sublimation energy is usually used as SBE in BCA simulations. They may, however, differ considerably since the former one is a bulk property whereas the SBE is defined as the binding energy of a single-atom to the surface. Besides, MD can provide information about sputtering mechanisms [24]. MD results strongly depend on the interatomic potential energy function. Earlier studies used so-called embedded atom (EAM) potentials. Their parameters were obtained from density functional theory calculations (DFT) on tungsten structures [25]. Despite being quite efficient, EAM potentials do not reach DFT-level accuracy themselves, due to their restricted flexibility. Machine learning potentials can achieve such accuracy. Previously, Gaussian approximation potentials (GAPs) [26] and spectral neighbour analysis potentials [27] were used for this purpose. We use a similar technique, highdimensional neural network potentials (HDNNP) [28] in the present work.
We model the non-cumulative sputtering from tungsten surfaces by Ar projectiles, as well as retention events. Based on former works [29][30][31], we continue to develop neural network potential energy functions using the Behler-Parrinello approach [28] as implemented in the n 2 p 2 code [32,33]. With neural networks representing the potential energy surfaces, sputtering, reflection, and retention probabilities for various incident angles α (0°, 20°, 40°, and 60°) on bcc-W(110) are calculated for projectiles with 100-800 eV kinetic energy and surface temperatures of 300, 1500, and 2500 K.

NN-driven MD simulations 2.1 Behler-Parrinello method
In HDNNP [28], the total energy of a configuration is obtained as a sum of atomic energies Ei , which are output of an atom-centred and element-specific neural network. A cut-off radius of 7Å defines the atomic environment. The Cartesian coordinates of the neighbour atoms are converted into parameters of weighted radial and angular Behler-Parrinello type symmetry functions [34] which reproduce the local environment with physically correct invariances. They serve as input to the feedforward neural network with two hidden layers of 25 nodes each. The nodes consist of a soft-plus activation function with an offset. After the NNP is created and trained, the MD simulations are performed with a modified LAMMPS code [35].

Density functional calculations
The energies and forces from DFT calculations serve as training and test data to determine the NN parameters, both in the initial and refinement steps. One can view the NN as a 'bridge' to generate fast DFT-quality forces that drive the MD simulation.
Total energy calculations and Born-Oppenheimer ab initio molecular dynamics generate the initial and further training data sets [36,37]. They employ Kohn-Sham DFT on the generalized gradient approximation (GGA) level with the PW91 functional [38]. Vienna Ab Initio Simulation Package (VASP) [39,40] was used for this purpose. A rather large supercell, with a Gamma-centred k-point mesh of 3 × 3 × 3 was employed. The Kohn-Sham orbitals are expanded in a periodic plane wave basis set. Spin-unpolarized calculations were performed with a cut-off energy of 350 eV [29]. GGA-PAW and the other methods had been used before for studying the interaction of noble gas atoms with tungsten [41]. As described below in detail, the initial set of training structures, energies, and forces were obtained from sputtering trajectories with different incident angles of Ar on small supercells. In optimizing the geometries of these cells, the lowest layer of atoms in the slab was fixed while all other atoms were allowed to relax. In this step, the convergence threshold for the forces was 10 −4 eV/Å.

Generation of training data and network training
In order to generate the training data, sputtering simulations were performed using direct density functional molecular dynamics simulations first. The initial training data consisted of 2712 configurations of Ar deposition on a W(110) surface with 72 and 144 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 [42,43]. The projectile Ar atom was initially placed at z = 5Å above the target surface with uniformly chosen x and y random coordinates and supplied with the initial velocity corresponding to the kinetic energy at impact.
We used an iterative refinement protocol as described in Ref. [29] to successively improve the training by supplying more data sets. The final set of reference data consists of 6230 configurations containing 72 and 144 W atoms and one Ar atom each. The final NNP is trained with 1,140,705 forces and 5655 total energies. The remaining 10% of structures containing 575 energies and 116 817 forces served as a test set to validate the interpolation capability of the neural network potential energy expression. With 50 training steps, the RMSE in the test set converged to 0.97 meV/atom for energies and 0.19 eV/Å for atomic forces. The corresponding values in the training set are 0.94 meV/atom and 0.17 eV/Å. There was no evidence of overfitting in the training process.
The correlation between NNP and DFT energies per atom and the corresponding ones for the atomic forces together with their distributions are shown in Fig. 1. It can be seen that the correlation of predicted vs. reference energies and forces is uniformly accurate for both within the large range spanned by them.

Simulation details
The initial simulation box was created from 630 tungsten atoms in a body-centred cubic lattice (bcc) with the (110) surface in the z direction. Its dimensions are x = 22Å, y = 22Å and z = 36Å. It was relaxed according to the initial temperature of the sample in the canonical (NVT) ensemble for 0.1 ps. The subsequent sputtering simulations were performed in the microcanonical (NVE) ensemble as described in Sect. 2.2.2. The timestep was dynamically adjusted between 0.001 and 1 fs depending on the initial velocity of the projectile. The simulation continues until one stop conditions is met. These conditions involve, for example, the exit of a sputtered particle or reaching a limit of timesteps. At each energy 2000 non-cumulative irradiations were simulated. The trajectories were visualized with the OVITO software [44].

Sputtering yields at perpendicular impact
The irradiation of Ar atoms on a tungsten surface with a surface temperature of 300 K is modelled by various incident ion energies and impact angles by MD simulations. Table I reports the sputtering yields, reflection (reflection of Ar atoms per total impacts) and retention probabilities (adsorption of Ar per total impacts). The statistical error of the sputtering yield was estimated as Δ = σ/ √ N , where σ is the standard deviation from a Bernoulli distribution. N is the number of MD trajectories over which we average.
In Fig. 2 our results from perpendicular impact are compared with available experimental data [11][12][13][14]. The yield is defined as the number of sputtered W  [20]. At high impact energies, computed and experimental values agree perfectly. At 100 eV, the calculated yield is higher than the experimental one. At the same time, at lower energies, the errors are larger due to poor statistics at low impact energies. We believe that the computed yields predict rather accurately the experimental ones if the sputtering yields are integrated over all incident angles. This will be discussed in the next section. Figure 3a visualizes the sputtering and reflection processes. Figure 3b (as well as Table 1) shows the derived dependence of W sputtering yield on various Ar incident angles and energies for all incident energies, together with the Eckstein fitting formula. Sputtering yields differ only slightly between the normal incident angle (0) and 20°. However, at 100 eV, increasing the angle leads to a decrease in the yield, and at 60°, no tungsten sputtering was detected in the simulations. At 150 eV, the yield is lowest at 60°and equal at other angles. In the range of 200 and 500 eV, the highest yield is observed at 40°, but not more than ≈20% with concerning the normal incident. At 60°, the yield is minimal. As can be seen, the angular effect is especially strong below 500 eV, and starting from 500 eV the angle effect is small. Based on this, it can be concluded that a reasonable range for modelling the angular effect of sputtering in our work lies below the keV region. Moreover, in the keV region, the timesteps of MD simulation should be reduced to very small values. This, in turn, significantly increases the simulation time. In this case, it is preferable of necessary to use other methods, such as BCA.

Angular distributions of sputtered particles
The distributions of the outgoing angles of sputtered W atoms are shown in Fig. 4. Outgoing angles plotted relative to the surface normal (in the z direction). The results are presented for all impact energies at 4 incoming angles. The shape of the distributions does not depend on the incident angle at 0 and 20°and very slightly for higher angles. The vast majority of the sputtered particles leave the surface at angles between 30°and 60°. The maximal yield shifts from 45°to 55°w ith an increasing angle of incidence. The distributions resemble a Chi-square probability distribution. We also analysed the angles between sputtered particles and unit vector along the x-and y-axes. Graphs are presented in the Supporting Information in Figs. S3 and S4. It should be noted that the sputtering angles along the x -axis have two maxima. While along the y-axis, Fig. S4, this effect is not so strong. This can be explained by the asymmetric surface structure of W(110). On such surfaces, the nearest neighbour atoms are located at distances a 0 and √ 3 2 a 0 (a 0 is the lattice constant) along the x-and y-axes, respectively. Thus, second maxima correspond to the collisions with distant surface atoms in each dimension.
The mechanisms governing sputtering seem to be different for the Ar case, compared to light atom irradiation [45]. This is not unexpected, due to the mass and size of Ar which leads to negligible penetration at low energies and low penetration at higher energies. In turn, the sputtering mechanism is simpler than in the case of light atoms, since sputtering occurs due to collisions of Ar atoms with surface atoms where it knocks out a surface atom and creates a defect on the surface. Some defects observed in this study can even be regarded as nanochannels with several W atoms involved. The results have been partially confirmed by an experiment [8].

Ar reflection and retention probabilities
In a way similar to sputtering, we also analysed argon atoms reflected after a collision. In Fig. 5a their reflection probabilities are shown in the same way as before as a function of the incident energy for different incident angles. The 'flat' incident angle of 60°causes by far the highest probability of reflection, which decreases regularly towards the angle approaching the surface normal. One also sees that the maximal reflection probability occurs at lower energies. The lowest probability is around 0.47 at higher energies. The decrease in reflection probability on the high-energy side can be attributed to the increased probability of retention of  the projectile on the surface. Figure 5b compares these probabilities. Reflection and absorption compete with each other. At higher energies, this quantity tends to have a constant value equal to ≈0.5. Together with adsorption, increase retention in the bulk or penetration ratio of Ar projectiles. At 800 eV maximal retention is about 0.1.

Distributions of reflected Ar
Equivalent to the sputtering case, we also analysed the angular distributions of reflected Ar atoms. The respective polar plots are shown in Fig. 6. Each panel shows angular distributions for all simulated energies at the same incoming angle. We see that the distributions are relatively wide for small incidence angles. The maxima   of reflection distributions correlate to incoming angles at 40°and 60°. With increasing energy, the contribution from other reflection angles increases. The maxima of the Ar reflection probabilities are correlated with the incoming angle, in contrast to what has been found in the sputtering case, just with some blurring reminiscent of diffuse reflection. As mentioned above (Fig. 5a), by far the highest reflection probability is found for the 60º impact angle. As can be seen from Fig. 6d, at this angle the probability of elastic ('specular') reflection is the highest. In contrast to sputtering, the reflection angles along the x-and y-axes are qualitatively the same.

Temperature effect
The effect of the PFM surface temperature on sputtering and related events warrants examination due to its practical importance. Depending on the temperature, both the properties of the material and the behaviour of the plasma particles can change. It is inaccurate to argue that the surface temperature would play a minor role because the temperature formally corresponding to the kinetic energy of the impinging particle is orders of magnitude higher. For instance, it was found that retention properties of a D/He mixture increased by a factor 10-100 depending on the tungsten sample temperature [8]. This shows the significant effect of the sample temperature. Furthermore, diffusion and penetration properties also depend on temperature [46]. In addition to the simulations at 300 K mentioned above, we studied the effect of sample temperature by also performing non-cumulative simulations at a perpendicular incident at 1500 and 2500 K. The trajectories were analysed in the same way as has been done at 300 K. Figure 7a shows the calculated sputtering yields. As can be seen, in the lower energy range, the temperature effect on our system is insignificant. Starting from 500 K an increase in temperature leads to an increase in the yield. We would expect sufficient differences in the keV range. The reflection probabilities show that at low particle energies a smaller number of atoms are reflected at a high surface temperature, which changes to the opposite for particle energies above˜600 eV (Fig. 7b). Figure 8 shows the angular distribution of sputtered W atoms at 800 eV and reflected Ar atoms at 100 eV for all surface temperatures. These energies were chosen as the most statistically reliable points. The shape of the distributions is similar to the ones at 300 K even though the sputtering yield at the two elevated temperatures is higher than at 300 K and a broader range of angles is found. The polar plot of the sputtering and angle of Ar reflections shows the trend of decreasing reflection probabilities at higher T observed above (Fig. 7). An increase in temperature increases mobility of surface atoms, which in turn makes W_surface-Ar collisions less elastic. It also causes Ar atoms to acquire additional velocity components in x and y directions, widening the distributions.

Summary
In this work, an accurate HDNNP was trained for the W-Ar system. This potential energy function was applied to non-cumulative sputtering simulations of the W(110) surface by Ar atoms with kinetic energies ranging from 100 to 800 eV. Various incident angles as well as three sample surface temperatures were investigated. The calculated sputtering coefficients at normal incidents are in good agreement with both experimental and predicted theoretical data. Reflection coefficients were also calculated. The angle of incidence becomes more and more insignificant with increasing impact energy but the effect on the reflection probability remains. The highest elastic reflection probability was found at 60°. The distribution of the outgoing angles of the sputtered particles has a maximum of around 45°. Reflected atoms still retain a slight correlation with their incoming angle. Increasing the surface temperature increases the yield at high incident energies and expands the angular distributions of both sputtered and reflected particles. All raw data are available from the authors upon request.

Author contribution
SS and MP contributed to the design and implementation of the research, to the analysis of the results and to the writing of the manuscript.
Funding Open access funding provided by University of Innsbruck and Medical University of Innsbruck.

Data Availability Statement
This manuscript has associated data in a data repository. [Authors' comment: The authors declare that all data supporting the findings of this study are available within the article and its supplementary information files].
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecomm ons.org/licenses/by/4.0/.