Sputtering from rough tungsten surfaces: Data-driven molecular dynamics simulations

Knowledge of the interaction of plasma particles with plasma-facing materials (PFMs) is crucial for the successful operation of fusion reactors. Tungsten (W) owing to its high melting point has been chosen as the primary armor metal for the divertor region of the ITER machine.^1,2 Particle bombardment causing sputtering, reflection, and retention dominates the interactions of fusion plasmas with tungsten. These processes can lead to surface and bulk modifications.³ Surface morphology changes of W caused by various ion bombardments have been studied extensively. They include the formation of a fibrous nanostructure commonly referred to as “fuzz,” which has been observed at relatively low incident energies.⁴ In particular, tungsten fuzzy structures were discovered in helium irradiation experiments,⁵ and an overview of available data is given by Reinhart et al.⁶ These needle-like structures have since become the subject of much experimental and theoretical research,⁷ for example, by scanning tunneling microscopy. The main driving force of the formation of such structures is the diffusion of surface atoms.⁸ However, once such structures are formed, the response of W surfaces with a fuzzy layer with respect to sputtering is poorly understood since experimental techniques are very sensitive to surface changes during measurements. Nishijima et al.⁹ performed argon bombardment experiments on smooth and fuzzy W surfaces. They found that the sputtering yield from fuzzy surfaces is an order of magnitude lower than the one from smooth samples.

Computational studies of ion bombardment on W-fuzzy surfaces is mainly done by BCA-based Monte Carlo methods.¹⁰ Single-crystalline, poly-crystalline, amorphous tungsten, and tungsten with fuzzy layers interacting with various energetic plasma particles were examined.^11–13 Higher-level micro-scale simulations of rough tungsten surfaces have been detailed studied by Eksaeva et al.¹⁴ These computational works cover many scenarios of interest to fusion research, such as sputtering yields, reflection and retention properties, and penetration depth of plasma particles. Despite the differences between these systems, they all agree that surface modifications result in reduced sputtering yields.

In this study, we aim to gain a better understanding of the interactions between energetic plasma particles and tungsten surfaces with fuzzy layers. We performed molecular dynamics (MD) simulations at the atomic scale to investigate interactions between argon atoms and imperfect tungsten surfaces. MD simulations explicitly model the interactions between individual atoms, considering forces such as van der Waals forces, electrostatic interactions, and bond breaking and formation. In contrast, BCA-based models allow for faster modeling but somewhat simplify these interactions by using empirical formulas and collision cascades. MD simulations require only the potential energy surface of the system. Our simulations use neural network potential energy functions, which have been proven to be a powerful tool for modeling complex system.^15–17 Specifically, we investigate sputtering, reflection, and retention mechanisms.

We employ non-cumulative molecular dynamics simulations in which the effect of a defect is determined by averaging over many trajectories, differing by the initial conditions of the Ar atoms and the microstate of the W sheet. In our protocol, we work with pairs of trajectories. The first simulation of each pair starts with the Ar atom moving toward a planar surface. The x/y coordinates of the Ar atom are chosen randomly, while z, the initial distance from the surface, is always 5 Å. The trajectory proceeds in the same way as in our previous work,¹⁷ leading to sputtering, reflection, or adsorption. For the second simulation of each pair, we identify where Ar encountered the surface in the first trajectory. This determines the center of the region where we construct a defect. The defect is constructed by removing surface atoms, and the second trajectory is constructed. We call these trajectories “mode II.” We believe that they are a good model for fuzzy tungsten where the whole surface is spiky and nonplanar. The first part of the simulations is shown in Fig. 1(a). The point of impact is represented by an orange dot in Fig. 2(b), where a hole or defect is introduced in the surface (blue filling). In Fig. 2(c), the second trajectory is constructed, and the Ar atom interacts with the surface defect. This approach enables us to simulate a surface and to accurately analyze the differences between planar non-planar geometries.

We also study other geometrical shapes, namely, cylindrical shapes with radius r and height h [Fig. 2(b)], and conical shapes defined by base radius r₁ and height h [Fig. 2(c)]. To perform the simulations, we use a neural network potential energy function that was trained with an extensive dataset of structures computed through the density functional theory (DFT). Neural network potentials enable accurate reproduction of potential energy surfaces in complex systems, achieving DFT-level precision in our molecular dynamics calculations. Detailed information regarding the training and validation procedures can be found in our previous research.¹⁷ To obtain a reasonable statistical quality, 3000 (1000 for incident energies 500 and 800 eV) pairs of non-cumulative simulations of Ar irradiation each are performed, all at normal impact. The impact energies range from 60 to 800 eV. The geometrical characteristics of shapes used in this work are as follows: sphere r = 6.4 Å, cylinder r = 5.32 Å, h = 10 Å, cone r₁ = 6.4 Å, h = 10 Å, and r₂ = 3.72 Å. For the hemispherical defect [Fig. 2(a)], we also investigate the effects of changing the radius r. Although these shapes are by necessity only models for real tungsten surfaces in fusion reactors, they can capture some of their characteristics.

In addition to the MD simulations performed according to the protocol described above, simulations where simply a defect is present in the center of the periodic slab were carried out. In them, the trajectory of the Ar atom starts again at random x/y coordinates and z = 5 Å. Ar will then normally not collide with the slab at the center of the defect (as is always the case in for the second trajectory of the protocol described above) but at its side or even at planar surface region. For brevity, we refer to these simulations with a static defect as “mode I.”

Figure 3(a) shows sputtering yields, defined as the number of sputtered W atoms per single Ar ion impact, together with experimental data. Our approach involves dynamic changes in the position of defects according to Ar initial positions, making them the opposite of a predefined hole. Thereby, mode I refers to the simulation with a permanent defect, while the smooth line corresponds to the planar surface and mode II to the second stage of the protocol, as described in Sec. II. The dashed line represents the sputtering theory or Eckstein sputtering formula,¹⁸ which acts as a qualitative standard. Since non-cumulative simulations generate poor statistics, the statistical error of the computed values was estimated through bootstrapping. As can be seen, the sputtering yields partially agree with experimental data measured by Nishijima et al. In Fig. 3(a) (left side), a cone-shaped defect is applied. Despite large error bars, the curve from mode II is lower than the others at all energies. Little effect is seen for mode I simulations, compared to the flat surface. The experimental data from Nishijima are only available for low energies. Their results for fuzzy surfaces agree well with our mode II simulations. Their curve for smooth surfaces reflects lower yields than we found. The experimental difference between fussy and smooth surfaces is as clearly present as in our simulations.

Figure 3(b) compares the sputtering yields from the three differently shaped defects. Significant differences are mostly observed at lower energies but are always smaller than the respective error bars. The conical shape yield is the lowest, and the yield from a spherical surface defect is the highest. We can conclude that a nonplanar surface always shows a decreased sputtering yield, compared to a planar one and that the exact shape of the defect matter less, at least within the three geometries investigated. It might seem obvious that sputtering into the “free space” is inhibited if occurring inside holes, compared to a planar surface, but the changes in surface binding energy and other factors would make quantitative predictions without numerical modeling very difficult.

Enlarging a defect increases the sputtering yield in all cases. We used spherical defects with varying radii to illustrate this pattern. In Fig. 3(c), this can be seen from the sputtering yields of surfaces with spherical defects of three different radii. At low energies, the differences are pronounced, although the statistical uncertainty is highest. In short, the yields from surfaces with larger defects resemble those of smoother surfaces.

In the same way as with sputtering, the trajectories can be analyzed with respect to reflection events. In Fig. 4(a) (left side), the reflection probabilities of argon atoms upon collision with the different surface samples are presented. One sees that the mode I scenario (impact at the center of the defect) has also the lowest probability for Ar to be reflected. The difference between planar and mode II (static defect) on one side and mode I is even larger than for sputtering events. Obviously, the defects can trap Ar atoms, which release their energy in the bulk, leading to a less reflection. Figure 4(b) (right side) compares the reflection probabilities between the defect geometries. The differences are somewhat larger than for sputtering. At low energies, the defect geometry has more influence than in the case of sputtering. At low impact energies, the hemispherical defect has the lowest reflection probability and the conical one the highest (the least difference from a planar surface). At high energies, difference between the defect geometries becomes smaller (∼5% compared to ∼20% at low energies). The cylindrical shape has the lowest reflection probability, and cone and hemisphere become equal.

The shapes of the distributions vary between the planar surface, the one with defects and the other within the latter. On a smooth surface, W atoms are primarily sputtered at an angle of slightly above 45° in a rather narrow distribution. Our distributions for defected surfaces differ from the ones previously analyzed for smooth surfaces.¹⁷ These distributions are close to previous observation made by Nishijima et al.⁹ and are also supported by other simulations.⁶ It was noted that at low impact energies, on a fuzzy surface, most of the W atoms are sputtered at angles from 0° to 45°, in contrast to the distribution between 0° and 60° observed on a smooth surface.

In summary, surface defects induce a shift of the peak of angular distributions toward smaller angles, with the most significant shift occurring on conical defects. This shift is associated with a narrower distribution. Cylindrical defects, on the other hand, result in a wider distribution with the peak centered around 35°. The widths of the distributions from spherical defects and smooth surfaces exhibit only minimal differences.

In this study, molecular dynamics simulations were employed to investigate the effect of surface roughness on the erosion yields of tungsten surfaces. Our models used conical, cylindrical, and spherical defects. The calculated sputtering yields could be compared with experimentally obtained yields from fuzzy and smooth samples. The simulation shows that the surface shape plays a significantly different role at lower energies but less so at higher ones. The reflection probabilities are also lower for non-smooth surfaces, and the shape of the defect significantly affects the them. Finally, angular distributions of outgoing W atoms and of reflected Ar atoms vary between the surface structures in a complicated way.

Our simulation results point toward both an advantage and a disadvantage of fuzzy structures in terms of sputtering and reflection probabilities. The advantage is the stability against sputtering, and the disadvantage the increase in Ar uptake due to less reflection. Although we use idealized defects and, mostly, sputtering and reflections caused by impacts at the center of them, we believe that our approach can be deemed valid if surface defects are the rule rather the exception like it is the case in fuzzy materials.

The work has partially been carried out within the framework of the EUROfusion Consortium. MP received funding from the Euratom research and training programme by Grant Agreement No. 101052200-EUROfusion. SS has received funding from the European Union's Horizon 2020 research and innovation programme under the Marie Skłodowska-Curie Grant Agreement No. 847476. The views and opinions expressed herein do not necessarily reflect those of the European Commission. The computational results have been obtained using the HPC infrastructure LEO of the University of Innsbruck and the Vienna Scientific Cluster VSC4.

The authors have no conflicts to disclose.

Shokirbek Shermukhamedov: Data curation (equal); Formal analysis (lead); Investigation (equal); Methodology (equal); Software (lead); Validation (lead); Visualization (lead); Writing – original draft (lead). Michael Probst: Conceptualization (lead); Funding acquisition (lead); Methodology (supporting); Project administration (lead); Resources (lead); Supervision (lead); Writing – review & editing (lead).

The data that support the findings of this study are available from the corresponding author upon reasonable request.