F-TRIDYN: A Binary Collision Approximation Code for Simulating Ion Interactions with Rough Surfaces

TL;DR

F-TRIDYN is an advanced simulation tool that models ion interactions with rough surfaces using fractal geometry, improving understanding of sputtering phenomena and enabling better plasma-material interaction predictions.

Contribution

It introduces a fractal surface algorithm into TRIDYN, enhancing simulation efficiency and accuracy in modeling surface roughness effects on ion irradiation.

Findings

Reproduces experimental sputtering data from rough surfaces

Shows roughness significantly affects sputtering yield and angles

Provides a flexible tool for multiscale plasma-material modeling

Abstract

Fractal TRIDYN (F-TRIDYN) is a modified version of the widely used Monte Carlo, Binary Collision Approximation code TRIDYN that includes an explicit model of surface roughness and additional output modes for coupling to plasma edge and material codes. Surface roughness plays an important role in ion irradiation processes such as sputtering; roughness can significantly increase the angle of maximum sputtering and change the maximum observed sputtering yield by a factor of 2 or more. The complete effect of surface roughness on sputtering and other ion irradiation phenomena is not completely understood. Many rough surfaces can be consistently and realistically modeled by fractals, using the fractal dimension and fractal length scale as the sole input parameters. F-TRIDYN includes a robust fractal surface algorithm that is more computationally efficient than those in previous fractal codes and which reproduces available experimental sputtering data from rough surfaces. Fractals provide a compelling path toward a complete and concise understanding of the effect that surface geometry plays on the behavior of plasma-facing materials. F-TRIDYN is a flexible code for simulating ion-solid interactions and coupling to plasma and material codes for multiscale modeling.

Full text

††journal: Journal of Nuclear Materials

1 Introduction

Plasma-material interactions (PMI) are crucial to the operation of all plasma devices, including experimental fusion reactors. In fusion experiments, interactions with the the wall are responsible for drastic changes in plasma performance and confinement, including radiation losses via impurities. Energetic ions are primarily responsible for many phenomena important to fusion reactor performance, including sputtering and tungsten fuzz growth [1]. In addition to semi-empirical formulas for calculating sputtering yields such as the Yamamura formula [2], the two main strategies for studying ion-solid interactions are Molecular Dynamics (MD) and Binary Collision Approximation (BCA) simulations [3].

While semi-empirical equations for sputtering yields can be used for simple systems, simulations allow one to study a significantly larger volume of the many-dimensional parameter space for ion-solid interactions. Simulations provide information, such as implantation depth profiles and Frenkel pair damage profiles, that cannot be captured in a single equation. MD simulations directly integrate the equations of motion for a number of interacting atoms. This direct integration, generally relying on ab initio or experimentally derived pair interaction potentials, is computationally difficult. MD simulation of timescales relevant to ion-solid interactions in the fusion community requires access to high-performance computers and, in some cases, prohibitively long computational times. For this reason, MD codes are not practical to use as directly coupled components of multiscale codes. BCA codes, on the other hand, treat the ion-solid interaction as a combination of elastic two-body ion-atom collisions and inelastic electronic stopping [4]. BCA codes are fast enough to run on-line with other codes and are statistically and physically accurate for a wide range of materials and ion energies [5].

While SRIM, based on the TRIM code, is arguably the most popular BCA code in use today, it has been shown to be inaccurate, especially for sputtered atom angular distributions [6]. TRIDYN, and its successor SDTrimSP, do not share the errors present in SRIM. They are upgraded versions of TRIM to include a more accurate interaction potential and features such as dynamic surface composition [5]. Fractal TRIDYN (F-TRIDYN) is a new, further upgraded version of TRIDYN to include explicit fractal surface roughness and additional output modes. Fractal surface roughness was first implemented in a BCA code in FTRIM [7],[8]. FTRIM, in addition to being based on the older TRIM code, utilized a fractal surface algorithm that was of $O(N^{2})$, where $N$ is the number of points that comprise the fractal surface, making it prohibitively expensive to use for large-scale contemporary simulations. The new surface algorithm implemented in F-TRIDYN, however, is of $O(N)$. In addition to the new fractal surface algorithm and the inherent advantages in accuracy of using TRIDYN as a base code, F-TRIDYN adds features that FTRIM does not have, including output lists of stopped projectile locations, stopped Primary Knock-on Atom (PKA) and Secondary Knock-on Atom (SKA) locations, Frenkel pair damage locations, and energy-angle distributions of the sputtered particles, all in three dimensions. These additional output features in particular make F-TRIDYN suitable for coupling to a wide variety of plasma and material codes as a wall boundary condition for energetic ions.

2 Fractal Surface Roughness

Fractals are mathematical objects that possess either self-similarity or self-affinity. Self-similarity is the property that elements of an object at different scales are exactly or approximately the same. Self-affinity is a weaker quality that allows for an object to scale differently along different spatial dimensions. Fractals can be characterized by a so-called fractal dimension. Fractal dimensions describe how the size of the space fractals occupy relates to the size of the space within which they are embedded. In particular, the fractal dimension of a curve describes how its total measured length changes depending on the length scale of the measurement. For example, smooth curves will be measured to have the same total length regardless of the length scale considered. Therefore, the space smooth curves occupy scales according to a power law of order equal to their topological dimension. Measurements of the length of fractal curves, however, will scale according to a power law of order greater than their topological dimension, meaning that their total length is not well-defined and depends on the observed length scale. In particular, they will scale according to a power law of order equal to their fractal dimension [9]. Intuitively, fractal dimension can be seen as a measure of the complexity or roughness of a curve. Determination of the dimension of simple geometric objects, such as points, lines, and planes is trivial. For fractals, however, there are multiple methods of measuring their dimensions, which are not, in general, equal. Throughout this paper, the fractal dimension referred to will be the box-counting dimension [10] as described hereafter.

Calculation of the box-counting dimension is straightforward. For a given curve, separate the space it occupies into boxes of side length $L$. The number of boxes that are at all occupied by any portion of the fractal curve are counted as $N$. If one plots on a log-log scale the number of occupied boxes versus the box side length, the slope of any linear section of the resulting line will be equal to the negative fractal dimension in that region. Some curves may exhibit fractal behavior over a limited range of length scales. Outside of this range the box-counting plot will either have slope zero or be nonlinear [9]. The range of length scales for which a curve has a measurable, non-zero fractal dimension will be referred to as the fractal length scale. In Fig. 1, the Koch Snowflake and its box-counting plot are shown.

Many natural objects and physical phenomena are fractal at some length scale [10]. Particularly, fractals provide a consistent and realistic model for rough molecular surfaces [11], [12]. For surfaces such as these, the roughness can be characterized by just two parameters, the fractal dimension and the fractal length scale. This allows complex geometry representative of the rough surface to be generated on the fly for any arbitrary fractal dimension. Additionally, fractal dimension is practically measurable in-situ. For most PMI experiments, it is not possible to measure surface morphology directly without breaking vacuum. Using the methods of Fractal BET adsorption theory [13] [14], the fractal dimension of a rough surface can be determined via the amount of adsorption onto the surface of differently sized gas molecules. This offers a compelling methodology to study the effect of rough surfaces on ion-solid interactions in a single machine by measuring surface roughness, performing ion irradiation, and performing subsequent surface diagnostics.

3 Generation of Fractal Surfaces for F-TRIDYN

Fractal surfaces in F-TRIDYN are composed of two orthogonal cross-sections. Because a fractal surface will have an overall fractal dimension equal to the sum of the dimensions of each of its cross sections, two-dimensional fractal surfaces of dimension $2.0\leq D_{2D}\leq 3.0$ can be represented by a combination of a dimension 1.0 cross-section and a dimension $D_{2D}-1.0$ cross section. For this reason, fractal surfaces in F-TRIDYN are represented using a single fractal cross section, spanning the $Y-Z$ plane. As long as the fractal dimension of a surface is homogeneous across the target, this is a reasonable approximation for a Monte Carlo code such as F-TRIDYN. Although surface roughness in F-TRIDYN is two-dimensional, the addition of complete spatial particle tracking in F-TRIDYN over TRIDYN allow for the possibility of three dimensional surface implementation in the future.

Fractal cross sections are created for F-TRIDYN using the fractal generator method popularized by Mandelbrot [9]. This method iteratively combines an initial piecewise curve, known as the generator, to construct a fractal curve. A true fractal would require iteration ad infinitum, but for the description of atomic surfaces the generator need only be iterated until the distance between points is approximately equal to desired minimum interatomic distance. A recursive algorithm is used to construct the curve as input to the code. The initial curve is translated, rotated, and scaled to the local angle and segment length for each line segment of the previous generation.

4 Implementation of Fractal Surface Roughness in F-TRIDYN

Rough surfaces are included in F-TRIDYN via a parallelized implementation of the even-odd algorithm developed for solving the point-in-polygon problem in computer graphics [15], [16]. To determine whether a given particle is inside or outside the surface, a ray is cast upwards from the particle’s position and the number of intersections it makes with the surface is determined. A simple fact, provable from the Jordan Curve Theorem [16], is that if the number of intersections is odd, the particle’s position lies inside the surface, and if it is even, it lies outside the surface. This algorithm is O(N), where N is the number of points in the piecewise curve. F-TRIDYN can alternatively use any arbitrary piecewise curve for the rough surface. Surfaces for this work were created exclusively using fractal rough surfaces created via the fractal generation algorithm.

Solid lattices in TRIDYN are treated as amorphous, and the locations of collision partners are chosen randomly from the number density of the local compositional layer [5]. If a collision partner is created above the surface in TRIDYN, determined by a simple inequality, the impact parameter for the collision is set to an arbitrarily high value such that the effect of the collision is negligible. The same method is used in F-TRIDYN, but the determination of whether the particle is inside or outside of the surface is performed using the methodology outlined above.

5 Results for Sputtering Yields

Figures 3 (a)-(c) compare the angular sputtering yield for two surface roughnesses of three ion-target systems: Helium on Beryllium, Deuterium on Beryllium, and Argon on Tungsten respectively. As observed in other simulations and experiments, the peak of maximum sputtering yield shifts to the right [17], [18], [19]. In Fig. 4, F-TRIDYN results are compared to an alternate simulation method by Küstner et al., namely using STM (Scanning Tunneling Microscope) images to reconstruct local angle of incidence distributions for the surface and using those distributions as inputs for TRIM, and an experiment. Sputtering yields are difficult to measure and prone to significant experimental error. Nevertheless, F-TRIDYN compares quantitatively favorably to Küstner et al.’s simulation technique and qualitatively recreates the trends seen in the experiment. The fractal dimension for this comparison was chosen heuristically based on the standard deviation of heights, as the fractal dimension in that study was not measured, only the deviation in surface heights. This leads to the possibility of further increases in accuracy if the fractal dimension and length scale are measured directly. In contrast to Küstner et al.’s simulation, F-TRIDYN allows one to characterize a surface with two parameters, and allows for one to simulate rough surfaces even when their exact morphology is unknown.

Fractal dimension has a nonlinear relationship to other measures of roughness, such as the RMS roughness. In particular, the standard deviation of heights of the points of a fractal surface created using the fractal generator method is proportional to the square root of the fractal dimension. The relationship between the standard deviation of heights squared and the fractal dimension is plotted for three different fractal generators in Fig. 6. The constant $a$ depends only on the shape of the fractal generator and the fractal length scale. With the limit of $1.0\leq D\leq 2.0$ for a fractal curve, this means that even small fractal dimensions can play a significant role in ion-solid interactions. For the same three systems as before, the effect of fractal dimension on sputtering yield for two angles of incidence is plotted in  5. With a fractal length scale of 40 Angstroms, the result for high angles of incidence is overall a slight reduction in yield.

Figure 7 (a)-(c) compare the sputtering yield versus energy for three surface roughnesses at normal incidence for the same systems as presented above. For most energies and systems, surface roughness will decrease the sputtering yield, as there is a greater chance that atoms emitted from the surface will be recaptured at a different location. For Helium and Deuterium on Beryllium, the sputtering yield increased with roughness for a large range of energies. Argon on Tungsten sputtering yields remained relatively constant with changes in surface roughness. Due to the speed of the fractal surface algorithm, reasonable simulation statistics to capture sputtering yields near threshold, an important energy regime for future fusion reactors, can be accomplished with F-TRIDYN.

6 Conclusions

F-TRIDYN is a flexible BCA code for the simulation of ion-solid interactions with a rough surface. It performs within experimental error for systems where there is available sputtering yield data that includes an explicit measurement of surface roughness. In order to further validate the fractal surface roughness model, however, significantly more experimental data is needed. Surface roughness is predicted to play a significant role in the performance of future fusion reactors as energetic ions drive significant morphological evolution in Plasma Facing Components (PFCs). F-TRIDYN is designed for simple coupling to other codes, adding 3D spatial tracking and output of projectiles, distributions of sputtered atoms, PKAs, and SKAs in addition to Frenkel Pair damage. As a component in a multiscale code, F-TRIDYN will be useful for simulation of the atomic timescale effects of energetic ions on a material, such as erosion, implantation, and backscattering, and will act as the boundary between plasma edge and material bulk codes.

Acknowledgments

This work was supported by the PSI-SciDAC Project funded by the U.S. Department of Energy through contract DOE-DE-SC0008658.