University of Tartu Faculty of Science and Technology Physics Institute Tarvo Metspalu CLIC electrical breakdown studies: Cu self-sputtering MD simulations for 0.1-5 keV ions at elevated temperatures Master’s thesis Approved for defense: Supervisors ....................................... Vahur Zadin, PhD ....................................... Flyura Djurabekova, PhD (signature, date) Konstantin Avchaciov, MSc Tartu 2015 Contents Acronyms and glossary of symbols 3 1 Introduction 5 1.1 Importance of the subject . . . . . . . . . . . . . . . . . . . . . . . . . 5 1.2 The aim and outlines of this thesis . . . . . . . . . . . . . . . . . . . . 6 1.3 Author’s contributions . . . . . . . . . . . . . . . . . . . . . . . . . . . 7 2 Overview of the subject field 8 2.1 Compact Linear Collider (CLIC) . . . . . . . . . . . . . . . . . . . . . 8 2.2 Electrical breakdowns . . . . . . . . . . . . . . . . . . . . . . . . . . . . 10 2.2.1 RF cavity EBD experiments . . . . . . . . . . . . . . . . . . . . 10 2.2.2 CERN DC-spark testing . . . . . . . . . . . . . . . . . . . . . . 12 2.2.3 EBD formation . . . . . . . . . . . . . . . . . . . . . . . . . . . 12 3 Effects concerning particle bombardment 14 3.1 Sputtering . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 14 3.1.1 Collision cascades . . . . . . . . . . . . . . . . . . . . . . . . . . 14 3.1.2 Sputtering yield . . . . . . . . . . . . . . . . . . . . . . . . . . . 15 3.1.3 Dependence on temperature . . . . . . . . . . . . . . . . . . . . 18 3.2 Nuclear and electronic stopping power . . . . . . . . . . . . . . . . . . 19 3.3 ZBL potential . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 20 3.4 Channelling . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 21 4 The Molecular Dynamics method 24 4.1 Basics principles of MD . . . . . . . . . . . . . . . . . . . . . . . . . . . 24 4.2 Interatomic potential and the embedded-atom model . . . . . . . . . . 25 4.3 Advantages and disadvantages of MD . . . . . . . . . . . . . . . . . . . 26 5 The simulation setup 27 5.1 Boundary and initial conditions . . . . . . . . . . . . . . . . . . . . . . 27 5.2 Incorporating ESP into LAMMPS . . . . . . . . . . . . . . . . . . . . . 29 5.3 Bombardment and sputtering yield . . . . . . . . . . . . . . . . . . . . 30 5.3.1 Incident angle considerations . . . . . . . . . . . . . . . . . . . . 31 5.3.2 Crystallographic orientations and averaged SY . . . . . . . . . . 32 5.3.3 Rough surface considerations . . . . . . . . . . . . . . . . . . . 33 1 6 Molecular dynamics simulation results 34 6.1 Cu EAM potential validation . . . . . . . . . . . . . . . . . . . . . . . 34 6.2 Angular dependence . . . . . . . . . . . . . . . . . . . . . . . . . . . . 36 6.2.1 Normal incidence . . . . . . . . . . . . . . . . . . . . . . . . . . 36 6.2.2 Small angle deviations . . . . . . . . . . . . . . . . . . . . . . . 39 6.2.3 Full angle range . . . . . . . . . . . . . . . . . . . . . . . . . . . 40 6.3 Surface roughness dependence . . . . . . . . . . . . . . . . . . . . . . . 43 6.4 Comparison with a previous model . . . . . . . . . . . . . . . . . . . . 45 7 Discussion 46 Conclusion and outlooks 49 Kokkuvõte (Conclusion in Estonian) 52 Acknowledgements 55 References 56 Appendix A: Source codes 60 Appendix B: Incorporating ESP into LAMMPS 64 Appendix C: Additional figures 65 2 Acronyms and glossary of symbols BCA — Binary Collision Approximation BCC — Body-Centred Cubic (crystal structure) BDR — (electrical) BreakDown Rate/probability CERN — European Council for Nuclear Research CLIC — Compact LInear Collider CTF — CLIC Testing Facility DC — Direct Current EAM — Embedded Atom Model/Method EBD(s) — Electrical BreakDown(s) ESP — Electronic Stopping Power FCC — Face-Centred Cubic (crystal structure) FGS — Fixed Gap System HIP — Hot Isostatic Pressing LAMMPS — Large-scale Atomic/Molecular Massively Parallel Simulator LHC — Large Hadron Collider MD — Molecular Dynamics NPT/NpT — constant pressure p and temperature T regime (volume V may vary) NSP — Nuclear Stopping Power NVE — constant energy E regime NVT — constant temperature T regime (volume V and/or pressure p vary) OFHC — Oxygen-Free High thermal Conductivity OM — Order of Magnitude PBC — Periodic Boundary Conditions PD — Penetration Depth RF — Radio-Frequecny SBE — Surface Binding Energy SP — Stopping Power SRIM — the Stopping and Range of Ions in Matter SY — Sputtering Yield TDE — Threshold Displacement Energy UH — University of Helsinki UT — University of Tartu ZBL — Ziegler, Biersack and Littmark (inter-atomic potential) 3 If not stated otherwise, the following symbols are used in the context of their defini- tion. Some symbols, that have not been defined here, have been defined inside the text. Ecoh — cohesion energy Ed — threshold displacement energy Ekin or Ek — kinetic energy of the ion fi — force (vector) acting on the i-th atom G or Gα — embedding energy functional mi — mass of the i-th atom mion — mass of the ion N — number of atoms in the system ri — coordinate vector (location) of the i-th atom ṙi or vi — velocity vector of the i-th atom r̈i or ai — acceleration vector of the i-th atom rcut — cut-off distance for a pairwise potential rij — distance between the i-th and the j-th atom rN = (r1, r2, ...rN) — complete set of 3N atomic coordinates R — surface roughness t0 — current time(step) ∆t — timestep length T — ambient temperature of the target surface uαβ or uαβ(ri, rj) or uαβ(rij) — pairwise potential energy function uN(r N) — N-body interaction potential U or U(rN) — total interaction potential of the system Unb or Unb(rN) — total interaction potential for non-bonded particles vion — total velocity of the ion vx, vy, vz — x-, y- and z-components of the velocity of the ion Y and δY — sputtering yield and its uncertainty z0 — height from the surface from which the atoms are considered to be sputtered zion — height from the surface from which the projectile starts moving to the target θ — incident angle of the ion (to the normal of the surface) λ — lattice constant ραj (rij) — the contribution to the electron charge density from j-th atom of type α at the location of the i-th atom 4 1 Introduction 1.1 Importance of the subject Nowadays, physics experiments at the highest energy range are mostly carried out in the Large Hadron Collider (LHC) at Geneva, Switzerland [1]. Nevertheless, there are numerous issues with high background signal levels, particle multiplicities, energy losses, unknown initial states of the collided particles etc [2]. In order to increase the precision and also the energy of such experiments, numerous institutes across the globe have been working on a new type of e+e− collider called the Compact LInear Collider (CLIC) since the beginning of this millennium. With this new collider, planned to be fully operational by the year 2030 situating next to LHC, scientists hope to reveal significant laws of fundamental physics beyond the standard model that have not yet been witnessed experimentally or maybe even theoretically, including measurements of the mass of the Higgs boson. [3] In order to bring this elaborate idea to life, a lot of physical, economic and engineering-related issues need to be surpassed. Carrying out all of the necessary ex- periments and tests to get the machine running would be tedious, time-consuming and at some point, even experimentally almost impossible or unreasonable. Fortunately, thanks to the vast advancements in computer technology, scientists are able to calcu- late the evolution of physical quantities in space-time faster and faster as the time goes by. On the other hand, as our understanding of physical processes improve, so do the models we construct get more and more complex — this is the frontier where science and computational speed race each-other day by day. All in all, computational models concerning CLIC components can help scientists save time and money on this massive project, while also getting a deeper comprehension of the underlying fundamental laws that take place in the processes they investigate. Although, special care has to be taken that the basis of these simulation models would stay physically sound throughout, as scientists and engineers rely on these data in order to make this dream a reality. The collaborating research groups of CLIC can by and large be divided into two divisions: people working on the detectors and people working on the acceler- ators themselves. Scientists at the Intelligent Materials and Systems Laboratory of 5 the University of Tartu (UT) and the Accelerator Laboratory of the Helsinki Insti- tute of Physics (UH) are the ones who have been working with the materials of the accelerator’s Radio Frequency (RF) cavities. These components are used to ramp up the energy of the particles that travel in the CLIC. As high voltages and therefore high electrical field strengths are used in these components, special fastidiousness has to be practised in order to avoid Electrical BreakDowns (EBDs), which are observed even at extremely smooth surfaces and very high vacuums after many trials of different RF cavity concepts. [4] 1.2 The aim and outlines of this thesis The aim of the current master thesis was to test the hypothesis whether the temper- ature drastically increases the self-bombardment process of Cu, main material of the RF cavities, producing more non-bound atoms upon which an EBD can form. For this purpose, in-depth research and numerous MD simulations were carried out. In addition to temperature dependence, other parameters were focused on including the crystallographic orientation of the target, the incident angle of the projectile and small scale surface roughness produced with a Kinetic Monte Carlo (KMC) [5] code. This will provide important information to the plasma formation model that has been devel- oped in the UH and used widely for the CLIC research. If any significant dependence is to be found, it would offer unprecedented ideas on how to reduce these unwanted harmful events during the accelerator operation. The master thesis begins with the overview of the subject, focusing on the CLIC, its accelerating RF cavities, the electrical breakdown events that take place, how they may form and how are they studied. After the introductory sections into the subject, the effects taking place during Cu self-bombardment are explained, followed by the simulation method and setup (the program and the methods used). This is followed by the results obtained in graphs, which are analysed in the section after that. Finally, there will be a conclusion of the work carried out and an outlook on this topic at the end of this thesis. In addition to the main content, there are three appendices: scripts used for the simulations, tutorial on incorporating an original add-on used in the simulation program and supplementary figures. All the references that have been placed before the punctuation mark have been used for the specific sentence or concept, otherwise, if they are after the punctuation mark, the references have been used for more content, if not stated otherwise. 6 1.3 Author’s contributions The author contributed to this thesis by: • researching and reviewing articles on sputtering, including its temperature de- pendence and interplay with other parameters; • planning and setting up, including coding scripts for the simulations; • validating the simulation model to produce physically acceptable results; • upgrading the Molecular Dynamics (MD) program LAMMPS (see subsection 5.2 and Appendix B); • conducting the simulations and logging the simulation results; • analysing the results and offering explanations on noticeable tendencies; • planning future simulations. 7 2 Overview of the subject field 2.1 Compact Linear Collider (CLIC) The Compact LInear Collider (CLIC) study is an international collaboration of 28 countries (53 institutes, including UT and UH) working on a concept for a machine to collide electrons and positrons (rather than protons like in LHC) head-on at energies up to 3 TeV [6] (see figure 1). Figure 1: Compact Linear Collider layout for nominal 3 TeV version [2]. The concept of the linear accelerator is used to avoid the loss of particle energy due to the bremstrahlung [2], which is observed e.g. in synchrotrons and LHC. Therefore, the particles in CLIC must be accelerated to their final energies only in a single passage through the machine, making the complex truly ”compact”. CLIC is based on the state of the art two-beam accelerator technology, where the main beam is propelled to high energy by an additional high current electron beam (drive beam) [7]. The concept involves RF structures — multicell linear travelling wave accelerator components — which are going to produce accelerating electrical fields as high as 100 MV/m (see figure 2) [8]. These RF components are not active-type meaning they are not constantly on, rather the RF power is generated by decelerating the drive beam in special Power Extraction and Transfer Structures (PETS). The generated power is in turn transferred to the main beam. [9] 8 Figure 2: An assembled high electric field RF accelerating structure manufactured by joining the disks together by Hot Isostatic Pressing (HIP) diffusion bonding [10]. The high electric field strength (ca. 20 times bigger than in LHC) combined with the associated power flow of more than 100 MW define many of the design challenges of CLIC accelerating structures. In addition to fatigue damage in the RF components due to pulsed surface heating, scientists and engineers have been also investigating unwanted Electrical BreakDown (EBD) events, which represent a second serious limit to the CLIC project in general. [11] Selected RF parameters for CLIC accelerating structures, with losses considered, are summarized in table 1. Table 1: Selected CLIC accelerating structure parameters [12]. Parameter Value Operating frequency 12GHz Average accelerating electric field strength 100MV/m Total RF pulse train length 156ns Luminosity ∼ 1034 cm−2s−1 Repetition rate 50 pulse/s Estimated power consumption 589MW 9 Figure 3: Basic RF cavity disk (d = 9 cm) composed of Oxygen-Free High thermal Conduc- tivity (OFHC) Cu (purity of at least 99.999%). The disk is diamond milled and turned with micron precision. The radial lines are the damping waveguides, the middle part including the hole and its surrounding perimeter is the iris. [10] 2.2 Electrical breakdowns Electrical BreakDown (EBD) is a sudden substantial increase in current through an insulating medium (in CLIC’s case, vacuum) resulting from failure of the medium to withstand an applied electric field [13]. In other words, EBD emerges from a rapid reduction of the insulator’s resistance when the voltage across it exceeds the breakdown voltage. Generally, this phenomenon can be optically observed as an electrostatic discharge through the insulator. When this occurs upon a surface, defects like electrical treeing, spot-melting, cratering and/or erosion may be observed. [11] 2.2.1 RF cavity EBD experiments Defects produced on the RF multicell component disks by EBDs have been mostly observed near the irises (see figures 3 and 4). Besides rendering the surface and there- fore the devices useless in the long term, EBDs are also responsible for deflecting the trajectories of accelerated beam, which is definitely unfavourable as the main beam must stay as focused as possible throughout the accelerator structure. [2] As can be see from figure 4, different materials such as molybdenum (Mo), tung- sten (W), but also stainless steel have been tested to replace the region most affected by EBDs (the iris and its close region). Although these materials exhibit greater break- down voltage, maximum stress and surface electrical field strength limits, they are 10 t Figure 4: Damage on iris after runs of the 30-cell clamped structures tested in CTF II. First (a, b and c) and generic irises (d, e and f) of W, Mo and Cu structures respectively. [12] more susceptible to fatigue, meaning these components would fail much faster (with less amount of cycles) and would have to be replaced more often than irises made of high-purity OFHC Cu. [12] This, in addition to the economical reasons, is why Cu has been considered as the sole material for the RF disks. Besides optical identification, EBDs have been measured electrically in various test set-ups as rapid reductions of the transmitted power and reflection of the incident power, but also acoustically as a result of the shock waves that are caused by the rapid discharges [14]. These events are always characterised by an emission of intense bursts of fast electrons (EKin ∼ 100 keV) with a build up time of around 20 ns. Moreover, the fast rise of gas pressure, emission of visible and UV light with the pulse being longer than incident RF pulse (about few ms) and emission of positive ions (EKin ∼ 0-200 eV) with the pulse longer than incident RF pulse may also be witnessed. [15] In order to reduce the likelihood of EBDs happening, different conditioning or material treatment mechanisms have been proposed [16], which reduce the surface roughness down to around 10 nm, which in turn lessens the field enhancement [17]. To further minimize the BreakDown Rate (BDR), electrical restrictions have been set, 11 because the higher the electric field strength and the longer the pulses are, the greater the BDR gets [12]. Despite the precautions, EBDs are still always present as stated before, which is why engineers have to take them into account. Scientists working on CLIC have calculated the acceptable BDR to be around 10−6EBD/(pulse/metre), where the distance represents the distance travelled by the particles inside the accelerator, while an ultimate goal of 3 · 10−7EBD/(pulse/m) at 100 MV/m loaded electric field strength at 230 ns pulse length has been suggested [12]. 2.2.2 CERN DC-spark testing Since the Direct Current (DC) and RF breakdown mechanisms exhibit similar as- pects [18], studies of the former have been carried out by CERN (European Council for Nuclear Research) researchers in order to test candidate materials and surface preparations, and have a better understanding of the EBD mechanism, as DC reruns are easier to set up and also cheaper to test. In a simple DC set-up upon which this thesis is mostly based on, there were two circular Cu electrodes in ultra-high vacuum at room temperature (T ≈ 300K) that were positioned face to face with only a 20µm gap between them. The DC discharge conditions were adjusted to match the RF case, so as 5 to 6 kV was applied on the electrodes, the average electric field strength rose to around E0 = 100 − 150MV/m, while the electric field strength close to the surface rose to around Esurf = 200 − 300MV/m increasing the field emission (explained in 2.2.3) and EBDs were observed. It has been noted from these experiments, that the choice of the cathode material has the biggest impact on EBD voltage, which in turn directly affects the maximum electric field strength threshold. This experimental set-up is sometimes also referred to as CERN Fixed Gap System or FGS in short. [18] 2.2.3 EBD formation The EBDs observed in a CLIC-like environment can be classified as spontaneous vac- uum (plasma) arcs or continuous discharges of quasi-neutral gas of charged particles that last as long as there is enough energy and particles fuelling the plasma [19]. The detected vacuum arcs in both the RF and DC set-up are triggered by field emission — emission of electrons from the surface into vacuum activated by high elec- trostatic fields [19]. Typically this occurs in local fields of 1 GV/m and above, being strongly dependent on the metal’s work function (≈ 4.5 eV for polycrystalline Cu [20]). 12 Table 2: Selection of parameters of the EBD event [4]. Parameter Value Critical local electric field strength Eloc 10− 11GV/m Current of field emission 10−10 − 10−9A Emission area (DC) - calculated 10−20 − 10−16m2 Current of a fully developed arc 10− 100A Voltage of a fully developed arc ∼ 100V Damaged cathode area after one EBD 10−10 − 10−8m2 It has been experimentally observed in DC case that EBDs always occur around a crit- ical local field strength of Eloc = 10 − 11 GV/m, which might refer to presence of nano-scale needle-like protrusions (with an aspect ratio of 50-100) or other field emit- ters, which are believed to be responsible for this field enhancement [15]. In addition, the high discharge current can also heat up the surface area where the arc ends, enhancing the emission of electrons and neutral atoms (neutrals). As these particles go through many collisions in a stochastic manner, Cu ions are produced from the neutrals and at a critical concentration a plasma arc forms. [15] Once the plasma is established (see table 2), the vacuum arcs are good at self- maintaining themselves as with it a plasma sheath, an arranged layer of electrons and ions, forms close to surfaces. The distribution of flow velocities of these constituents becomes balanced through a difference in densities, producing a sheath potential that decelerates electrons and accelerates ions. From this point onward, the cathode is continuously bombarded by ions — this leads to the aforementioned defect (e.g. craters, around 50-300 per disk [21]) creation which helps sustaining the arc. When the external energy source is exhausted, the vacuum arc extinguishes with it. [15] Nonetheless, it is yet not quite clear from where does the critical amount of neutrals come from in the first stages of plasma formation, during the field emission, when the ionisation probability is rather low. However, it is known that some neutrals must be accumulating above the field emitter as otherwise the plasma would never initiate — the ions would quickly return to the cathode. Moreover, the local electric field of 10 GV/m observed is not enough to get a direct field evaporation of Cu. It has been proposed that the increase of the temperature of the surface may intensify the amount of neutral Cu atoms emerging from Cu self-bombardment, moreover, through a reciprocity of the field emission electron current heating up the field emitter and the electric field pulling on the upper crystal layers of Cu. [4] 13 3 Effects concerning particle bombardment 3.1 Sputtering Sputtering is a radiation process whereby atoms are ejected from a target material due to bombardment by an energetic projectile (i.e. an atom, an ion or a cluster of atoms). Sputtering takes place when the kinetic energy of the incoming projectile is higher than conventional thermal energies (>> 1 eV) of the target atoms, which allows for an atom to overcome the Surface Binding Energy (SBE) barrier. When the projectile and the target consist of the same type of particles, the process is called self-bombardment. The minimum kinetic energy (of the projectile) required to permanently displace an atom in solid from its lattice site to a defect position is known as the Threshold Dis- placement Energy (TDE) Ed. As this energy consists of many components, it is usually higher than the lattice energy, the SBE (commonly equalled to the heat of sublimation in calculations) or the Frenkel (interstitial-vacancy) pair creation energy. [22] Moreover, TDE depends on the crystallographic orientation of the lattice in which case the minimum and the average threshold energies can be distinguished. Experi- ments have approximated that for Cu the average TDEs are Ed(100) = 27 eV,Ed(110) = 18 eV and Ed(111) = 29 eV [23], while Yamamura et al. have calculated Ed = 29.65 eV for the polycrystalline Cu [24]. In addition, for the neutrals to be actually sputtered, their kinetic energy normal to the surface must stay above the SBE (3.52 eV for Cu [25]) when it crosses the plane of the surface. If it is any lower, it causes the atom to refract to the surface. This is why it is easier to remove a target atom from an extended position above the surface, rather than from a flat surface. [26] 3.1.1 Collision cascades A collision cascade is a sequence of collisions between atoms induced by an energetic projectile. Two common mechanisms can be distinguished: linear collision cascade (low energy) and non-linear heat spike (high energy) sputtering mechanisms. [27] In the former instance, when the mass of the projectile is low and the target material has a small density, the rather infrequent collisions occur as two-body colli- sions [28]. This mechanism can be approached by the Binary Collision Approximation (BCA) at its most simplified form, whereby the many-body effects have not been taken into account (see figure 5). BCA can be easily compared with familiar billiard 14 Figure 5: Schematic illustration of binary collisions by K.Nordlund. The red dot represents the projectile, blue dots the secondary (knock-on) cascade generation atoms. ball mechanics, where the momentum of energetic atoms is dissipated along the path via collisions with generations of cascade atoms (see figure 6, left). Linear cascades are always present in bombardment events, but these can easily develop into thermal spikes when the initial energy of the projectiles is high enough (see figure 6, right). In the heat spike mechanism, when the projectile heavier and more energetic, and the target material is dense, the collisions between the particles can occur really close to each other that these can not be regarded as independent. The temperature in the region of these dense collisions may rise up to tens of thousands of Kelvins for a transient time causing a liquid/gas-like zone because of the energy concentration in a small volume. It has to be noted, that this temperature does not correspond directly to the thermodynamic equilibrium temperature, as no Maxwell-Boltzmann distribution of kinetic energies will initially form — this takes about 3 lattice vibrations. [26] For Cu self-bombardment, recoil energies from around 5 keV are highly likely to produce heat spikes, but the energy needed might even be lower, as long as a liquid-like zone is developed [29]. When the thermal spike takes place near the surface of the tar- get, a crater may form due to the liquid flow of atoms. This does, however, differ from the macroscopic mechanism of crater formation as long as the size of the projectile is below one hundred thousand atoms. [28] 3.1.2 Sputtering yield Sputtering is best described by the the Sputtering Yield (SY) which is the number of atoms emitted from the target per incident particle [22]. In the case of figure 6 the sputtering yield in the left case is 1, whereas in the right case of a sputtered dimer and 15 Figure 6: Schematic illustrations of a linear collision cascade on the left and non-linear heat spike sputtering on the right. The thick horizontal line illustrates the position of the surface and thinner lines the ballistic paths of the atoms. Red, orange and green dots illustrate primary, secondary and tertiary recoils or cascade generations respectively. In the thermal spike case, the linear cascade develops into a dense collision volume, where the atoms collide and move in complex paths. a single particle, the sputtering yield is 3. The sputtering yield depends on the incident angle, the energy and the mass of the projectile (see figure 7), but also on the mass and the SBE of the target atoms [30]. In addition, for a crystalline target as Cu, the orientation of the crystal axes or the crystallographic orientation (see subsection 3.4) in regard to the target surface and the moving direction of the projectile is also relevant [31]. Figure 7: SY data for Cu self-bombardment. On the left, dependence on the kinetic energy of the projectile (exp. results + fitted function at θ = 0 deg); on the right, calculated incident angle dependency functions based on the fit formula by Eckstein & Preuss [26, 32]. 16 At the time, no specific fundamental equations exist to describe the sputter- ing yield’s dependence on all of the variables mentioned before — there are only approximations developed upon experiments carried out throughout last five to six decades [33, 32, 34, 35]. These referenced models mostly work best when the elements in question are not too heavy (e.g. Au) nor too light (e.g. H, He), the projectile’s energy is not too low and the target is free of any additives. The mathematical-empirical model for SY at a normal incidence (θ = 0 deg) at room temperature (T = 300K) devised by Yamamura, Matsunami et al. presented in figure 7 has been widely spread in the sputtering community and proven to approximate SY experimental data well enough: √ Q(Z ∗2)a (M2/M1) Sn(Ek) Eth Y (Ek, θ = 0 deg) = 0.042 (1− )s, (1) Us 1 + Γke0.3 Ek where Ek is the kinetic energy of the projectile [eV], 0.042 is a numerical factor in units of Å−2, Q(Z2) is a dimensionless fit-function dependent on the atomic number of the target, a∗(M2/M1) is a dimensionless energy-independent fit-function dependent on the mass of the target atom M2 and the incident ion M1, Sn(E) is the nuclear stopping cross section (this phenomenon is explained in subsection 3.2), Us is the surface binding energy of the target solid, Γ is a factor dependent on Z2 and M1, ke is the Lindhard electronic stopping coefficient,  is the reduced energy function, Eth us the sputtering threshold energy (related to Ed) and s is a factor dependent on the target material [33]. On the other hand, Eckstein and Preuss have suggested that the formula for angular dependence of SY by Yamamura et al. does not ”...represent the available data for all cases, especially for problems that occur at low incident energies and for selfbombardment.” For this reason, they proposed a revised (Yamamura et al.) SY formula modified with additional fit parameters based on computer simulations. [32] What is more, yet newer SY fit functions have been submitted, which have aimed to become more and more generic, while getting more complex and requiring extra specific fit constants [35]. Despite this and the fact that these formulae are still only estimations and not experimental data, not many SY fit functions are required to provide a great benchmark for the results (more in sections 5 and 6) and advanced SY formulae developed in the future. 17 3.1.3 Dependence on temperature One of the pitfalls of the semi-empirical formulae mentioned is that they lack the dependency on the target’s temperature. Even so, many experiments concerning the SY versus temperature have been carried out and when it comes to Cu (as well as Al) as a FCC metal, most of the research results claim that the SY either stays rather constant or experiences a small-scale deviation with the increase of the ambient temperature of the target surface (up to 20-30%) [36]. Furthermore, it is said that with the increase of the temperature, bombarding the (100) and (110) planes (see subsection 3.4) at a normal incidence, no change is observed, while bombarding the target surface at the (111) direction exhibits a minor decrease in the SY, which confirms the anisotropic effect of the temperature. This is to be contrasted with other kind of materials, e.g. Mo or W (both BCC metals), which have shown linear increases in yield up to 40 % in the temperature range from 350 to 1000 Kelvins. [37] It is known that with the increase of the target’s temperature, the mean Penetra- tion Depth (PD) for the projectile becomes shallower. In theory, this should indicate that more collisions are taking place near the surface of the target, possibly ejecting more atoms and therefore increasing the SY. Moreover, with higher temperatures, the thermal vibrational energies of the atoms should also increase and the effective binding energy of the surface decrease, so as the projectile hits these atoms, they should be kicked out of their lattice more easily. However, the experiments suggest that this seems to be only a partial explanation. [26] It has been hypothesized that with increasing temperature, the kinetic energy of the projectile is dissipated more, i.e. there are less focused few-body collisions and the SY drops [37]. In contrast, the increasing effect that the ambient temperature might have on heat spikes has been excluded as its Order of Magnitude (OM) compared with that of a thermal spike is negligible. Also, annealing of the defects as a SY decreasing effect has been excluded as well as small defects move around intensively at any temperature from room temperature and above [38]. On the other hand, computer simulations have connoted the increase (up to ∼10 times in comparison with room temperature) of SY with increasing tempera- ture [39, 40]. Additionally, an exponential increase of the SY due to additional evap- oration has been observed near melting point temperatures, but these processes have to be distinguished from one another [31]. It has to be kept in mind, that in these experiments, the OM of the set-up and often the flux of the projectiles has been many times bigger than in aforementioned simulations [26]. 18 3.2 Nuclear and electronic stopping power As discussed in subsection 3.1, an energetic particle deposits its energy through multi- ple interactions with the atoms of the target — this energy loss is generally measured in terms of Stopping Power (SP), which is the decelerating force acting on these (charged) particles [41]. The linear SP of a particle moving inside the (target) material is quan- titatively equal to the loss of energy E per unit path length x: dE S(E) = − , (2) dx which is conventionally indicated with eV/Å, MeV/mm or similar units. Figure 8: NSP and ESP for Al-in-Al versus particle energy per nucleon by H.Paul [42]. The total non-relativistic SP consists of the Electronic and the Nuclear Stopping Power — ESP and NSP (see figure 8). The former (ESP) refers to the deceleration of a projectile due to the inelastic collisions between the projectile and the bound electrons in the target. This component of the general SP is prevailing in the shallow part of the target surface when the energy of the projectile (and the knock-on atoms) is high enough to excite electrons. It has to be noted that ESP is not a form of force or energy, but a description of (kinetic) energy losses to electron excitations. [41] As the charge state of the projectile traversing the matter may change frequently, it is inconvenient to describe all possible interactions taking place in between, which is why ESP is often expressed as a simple averaged function for different charge states. At energies lower than 100 keV/nucleon, it is difficult to determine ESP theoretically, which is why empirical fit functions are preferred in low-energy situations. The most thorough fit curves for ESP for specific elements can be obtained from SRIM database, 19 which have been averaged upon numerous experimental results. [43] The other component of the total non-relativistic SP is the Nuclear Stopping Power (NSP), which concerns the elastic collisions between the projectile and the knock- on atoms in the sample. It has to be noted, that NSP has nothing to do with nuclear forces, rather it depicts the interaction of the projectile with the nuclei in the target. NSP can be calculated from the repulsive potential energy between two atoms if its form is known, growing with the mass of the projectile. This component of the general SP is prevailing in the deeper part of the target surface. [41] At extremely high energies (GeV and more), relativistic radiative SP also should be considered, but this remains out of the scope of the current thesis and is not regarded. All in all, the total non-relativistic SP which can be expressed linearly as Ftotal(E) = Felectronic(E) + Fnuclear(E), (3) consists of averaged fits to experiments, which are not perfectly accurate, because these do not only depend on the type and density of the material, but also on the substance’s microscopic structure and SP cross-section. Nevertheless, these fit curves give an acceptable approximation on the OM of the effects taken into account. [41] 3.3 ZBL potential The NSP has a substantial effect when the distance between atom nuclei is small enough (less than an Å) for them to interact. At minuscule distances like this, the repulsive interaction can be regarded as Coloumbic, but at bigger distances (i.e. more than half of the lattice constant λ), the atomic orbitals screen the nuclei from one another. This is why it is useful to consider the Coloumbic repulsion which has been multiplied with a screening function ϕ(r/a), where ϕ(r)→ 1 when r → 0: Z 21Z2e V (r) = ϕ(r/a), (4) 4πε0r where Z1 and Z2 are the charges of the interacting nuclei, r the distance between them and a the screening parameter. This equation is also known as the Ziegler, Biersack and Littmark (ZBL) repulsive potential, which has been formulated by fitting an universal screening function to theoretical potentials calculated for a huge diversity of atom pairs. This potential is usually fitted to regular interatomic potentials used in simulations with a cut-off range 20 of around few Å-s when more accurate close-range nuclei interactions are required — this is generally the case with non-equilibrium bombardment events. [44] The ZBL screening parameter and the ZBL screening function can be further expressed as follows: 0.8854a0 a = (5) Z0.23 0.231 + Z2 and ϕ(x) = 0.1818e−3.2x + 0.5099e−0.9423x + 0.2802e−0.4029x + 0.02817e−0.2016x, (6) where x = r/a and a0 is the Bohr atomic radius 0.529 Å. In the case of Cu self-bombardment, where Z1 = Z2 = 29, the ZBL repulsive potential and the screening parameter simplify to the following forms: Z2e2 VCu(r) = ϕ(r/a), (7) 4πε0r 0.8854a0 0.8854 · 0.529Å 0.4683766Å aCu = = = ≈ 0.108Å. (8) 2Z0.23 2 · 290.23 4.338918925 3.4 Channelling Channelling is the phenomenon in which a projectile experiences much lower stopping power in certain crystallographic directions than others (i.e. channelling is anisotropic). Crystallographic orientations and planes are conventionally described by Miller indices according to the substance’s unit cell (Face-Centered Cubic (FCC) for Cu) edges [45]. The most common FCC planes usually investigated are the low-index (100), (110) and (111) (see figure 10). All of the information about channelling presented in this subsection was gathered from [46]. It is important to note, that the ”crystallographic face” or the (100), (110) and (111) faces in the context of this work is used to emphasise that the crystallographic planes investigated are also the faces that are bombarded. Due to energy conservation laws the moving particle’s path inside the crystalline solid may be constrained by the channelling as the projectile is steered through a ”corridor” of atoms in the crystal lattice (see figure 9). This is described by a so-called ”steering potential” which constantly repels the projectile away from the lattice atoms nuclei instead of allowing the moving particle to collide with them, or at least minimizes the collision rate. 21 Figure 9: Adapted exaggerated schematic illustration of (a) channelling and (b) blocking as a projectile comes in from the left. The red spot indicates the collision area. [46] Channelling has been evidently witnessed during both experiments and MD sim- ulations of mean penetration depth of a projectile bombarding monocrystalline Cu [36]. The most important parameters affecting channelling are the incident angle and the energy of the projectile. For FCC metals (e.g. Cu), the (110) direction has, on av- erage, the strongest and the (111) direction the weakest channelling effect, especially at θ = 0 deg. This becomes conspicuous when the crystal lattices are viewed down parallel tu the surface normal (see figure 10 on page 23). When the projectile penetrates deeper into the solid, it gradually looses energy and at one point is susceptible to scattering. On average, this happens much deeper inside the target material than for a similar setup at a random projectile incident angle. In this case, there are many layers of atoms between the cascade and the surface, in other words, more atoms that can dissipate the kinetic energy of the cascade. This, furthermore, means that at the surface, atoms have lower kinetic energies and are less open to be ejected. 22 Figure 10: Schematic of the a) (100), b) (110) and c) (111) planes in a FCC crystal. Above: planes inside a FCC unit cell. Below: the same faces looked down in the direction normal to the plane, where A indicates atoms in the first layer, B and C the following stacking positions. [47] 23 4 The Molecular Dynamics method 4.1 Basics principles of MD In order to investigate Cu self-bombardment close-up and statistically, Molecular Dy- namics (MD) computation method as a N -body simulation was used. The following explanation on basics of MD is based on [48], [49] and the lectures on MD given at UH in 2014 by K.Nordlund and F.Djurabekova. This will exclude the discussion on algorithms that perform and optimize the MD calculations. Molecular dynamics, at its very core, has been built up on solving classical New- ton’s equations of motion for N atoms in a predefined simulation cell. These equations for the force fi [eV/Å] and the location ri [Å] of the i-th atom can be described in the following simplified form : ∂ fi(t0) = mir̈i(t0) = − U , (9) ∂ri(t0) 1 ri(t0 + ∆t) = ri(t0) + ṙi(t0)∆t+ r̈i(t0)[∆t] 2 + ... (10) 2 These time-dependent calculations ((t0) denoting the parameter’s value at the current time), where the acceleration r̈(t0) of the i-th atom is dependent on the potential U (generalized form), are solved step-wise for each atom ending up with finite-sized evolution steps. The amount of time that passes with each step is called a timestep ∆t and it can either vary or be fixed depending on the simulation in hand. When the system obeys the ergodic hypothesis, macroscopic thermodynamic properties can be determined over specific statistic thermodynamic ensembles. The most important ensembles used in MD are the following: • NVE — the total energy E is conserved as well as the number of atoms N and the volume V , while temperature T and pressure p are unregulated; • NVT — T is controlled by a thermostat, N is constant; V and/or p may vary; • NpT — p is controlled by a barostat, N is constant; V and/or T may vary. Usually, there are a few thousand up to a couple of million of atoms investigated in a classical MD simulation as it runs for around t ∼ 103−106 timesteps, corresponding up to tens of nanoseconds in experimental time scale (this means that ∆t ≈ 10−15s). 24 4.2 Interatomic potential and the embedded-atom model The interaction potential U used in equation 9 for non-bonded particles (an atomic system) is traditionally split into 1-body (an externally applied potential field), 2-body (pairwise potential), ... N -body terms as: ∑ ∑∑ ∑∑ ∑ U Nnb(r ) = u1(ri) + u2(ri, rj) + ...+ ... uN(ri, rj, ...rN). (11) i i j>i i j>i N>(N−1) Usually, in the simplified cases, the interaction potential energy of atoms is con- sidered to consist only of the two first terms with higher order interactions excluded. Although this makes the system more crude, it can still represent the essential physics and be computed faster. The most common pair-potential used in MD is the semi- empirical Leonard-Jones potential model which approximates the interaction between a pair of neutral atoms or molecules. However, the pair potential scheme alone leaves out the critical essence of the physics of metal bonding including elastic constants, vacancy formation energy, melt- ing points and cohesive energy, which makes it a non-complete tool for metallic ma- terial simulations. In order to take these coordination-dependent many-atom aspects into account, the embedded-atom model/method (EAM) has been introduced. In this approach, ”the energy of the metal is obtained by embedding an atom into the local electron density provided by the remaining atoms of the system” [50]. EAM interactions can be summarised with the following equation for potential energy of the i-th atom: ∑ ∑ 1 ∑ Ei = Gα( ρβ(rij)) + uαβ(rij), (12) 2 i i 6=j i,j(i=6 j) where Gα is the embedding energy functional indicating the energy required to place the i-th atom of type α into the electron cloud, ρβ(rij) is the contribution to the electron charge density from j-th atom of type β at the location of i-th atom, and uαβ(rij) the electrostatic, pairwise interaction [51]. The embedding energy in this case is considered as the interaction of the atom with the background electron gas, which is created by its neighbouring atoms [50]. The EAM potentials used by MD software are usually obtained through semi- empirical approximations with experimental data like the prementioned many-body 25 interaction energies, with the most thorough of them also including ab inito calculations based on well-established laws of nature [50, 52]. These potentials consist of three parts: the pairwise potential, the embedding function(al) and the electronic density function — all of them given in one file as a discrete lookup table [53]. 4.3 Advantages and disadvantages of MD There are many substantial reasons why the MD approach was chosen over the BCA approach to investigate Cu self-bombardment. To begin with, the former method follows the complete atomic-level dynamics of the system based on a representation of the interatomic interactions, while the latter approximates the full atomic dynamics of a material by a series of binary collisions [54]. Furthermore, although the most commonly used BCA code SRIM [44] can simu- late linear collision cascades and sputtering yields thereof for all ion in all disordered materials up to projectile energies of 1 GeV, it excludes many-body effects and has a schortcoming when the material becomes amorphous during ion bombardment [55]. Additionally, regular BCA does not usually account for different crystalline ori- entations. This automatically excludes the influence channelling, grain boundaries (interface between different crystalline orientations) and surface roughness may have on ion bombardment. Even if these effects have been considered in a BCA model, they are usually added up as another set of variables in the two-body collision equation which already tries to encompass multiple interaction variables. [56] On the other hand, MD simulations are computationally more expensive — they take a longer time to compute. Besides, MD simulations have a short simulation time (and therefore, experimental time) constraint because of the numerical methods used for solving the motion equations, meaning that longer simulations become ill- conditioned as even the most minuscule errors accumulate during numerical integration. Moreover, if the timestep is too small, general timespan and useful computational power is lost, or if the timestep is too large, one may end up with a non-physical too discretely evolving system. This is where careful calibration and adjustment of different system parameters has to be made depending on the physical properties of interest. Nevertheless, MD simulations can be easily remodelled and controlled to represent configurations of interest and provide a lot of data that can be analysed. This and other advantages mentioned in this subsection is why the MD method was chosen. 26 5 The simulation setup 5.1 Boundary and initial conditions The Cu self-bombardment MD computer simulations were performed using the Large- scale Atomic/Molecular Massively Parallel Simulator (LAMMPS) open-source software by Sandia National Laboratories, as it is widely used and takes use of the Message Pass- ing Interface (MPI) for parallel communication [57, 58]. The corresponding LAMMPS input scripts used for most of the simulations are provided in Appendix A. In the beginning, a simulation box was created by defining the lattice type to be FCC with a lattice constant of λ. By establishing the dimensions of the box in unit cells, atoms were placed to the lattice points and their coordinates were recorded. In order to make the results more comparable, the simulation boxes were created so that their dimensions (in Å-s) in all directions stayed approximately the same for any chosen FCC Cu crystallographic orientation. The conversion for the (110) face and the (111) face simulation boxes were based on the dimensions of the (100) face simulation box of (X, Y, Z)unit cells, using the ceiling function on irratio√nal results:x(100) = X x (110) = X x 8(111) = dZ/ e3  √ √ y(100) = Y y(110) = dY/ 2 e√z y(111) = dY/ 2 e (13) √ (100) = Z z(110) = dZ/ 2 e z(111) = dZ/ 3 e To make the simulation box represent a physical system, it was allowed to reach a thermal equilibrium via thermal expansion using a integration algorithm based on the NpT ensemble. This slowly gave the atoms a distribution of velocities to statistically emulate temperature — during this phase the simulation box could expand in all directions. The temperatures studied in the simulations were the following: • 300 K — room temperature, at which the complex works; • 500 K — heated accelerator material (during operation); • 1200 K — material spots heated up by the field emitters and/or vacuum arcs. Next, Periodic Boundary Conditions (PBC) were introduced, meaning that vir- tual copies of the simulation box were added on each side of the original simulation box. This was required in order for the system to represent bulk material properties 27 Figure 11: Schematic 2-dimensional illustration of the simulation setup (see text). at a lower computational cost. In other words, atoms that resided near the boundaries of the simulation box ”saw” the atoms on the opposite surface and when they exited from one side, they entered from the other side. Otherwise, the system would have floated in space and had surface forces acting on each side and not only on the side open for the ion bombardment. In the middle of the equilibration process, the simulation boundaries (not the box of atoms itself) were expanded in the positive z-direction (see figure 11) and the system was allowed to reach a more stable configuration. These two steps (i.e. the relaxation or the equilibration) during which the atoms reached to their minimum potential energy were calculated (i.e. run) for 8000-20000 timesteps, once per ambient target temperature, crystallographic orientation and simulation box size configurations. Thereafter, thermally controlled boundary conditions (1-2 unit cells thick) were applied to the sides of the simulation box, leaving only the topmost surface uncon- trolled. The atoms in this layer kept around the predefined ambient temperature of the target surface as due to the PBC in x- and y-direction (see figure 11), the system could otherwise heat itself through virtual copies whenever energetic motion of atoms (shock waves of the cascade) is involved. In order not to restrict the physics by do- ing so, different simulation box size configurations were tested out so that the kinetic energy of the atoms reaching the thermally controlled walls would be less than ∼ 1 eV. By default, the simulation boxes were created as big as (30, 30, 30) unit cells (or 28 ∼ (110, 110, 110) Å) according to the (100) face for projectile energies up to 1 keV and (35, 35, 40) unit cells (or ∼ (125, 125, 145) Å) for 2 and 5 keV in order to let the collision cascade dissipate well inside the simulation cell. Finally, as the dimensions of the simulation box became stable, the bottom layer was fixed by confining its movement, to avoid the movement of the simulation box after being hit with an energetic projectile. Right before the projectile (i.e. the ion) was introduced into the simulation, the target was randomly, yet uniformly, shifted between zero and the maximum dimension length of the simulation box in the respective direction (x&y) in order to avoid the ions hitting the same place on the surface in every simulation. 5.2 Incorporating ESP into LAMMPS The MD program LAMMPS had one important drawback when it came to energetic particle bombardment: it did not include the Electronic Stopping Power (ESP) by de- fault, as this effect is considered insignificant at lower energies (see figure 8 on page 19) at which most of the classical equilibrium MD simulations are carried out. Although, Cu self-bombardment as a more energetic process is definitely a non-equilibrium MD process, which made it necessary to incorporate ESP into the LAMMPS source code. The additional ESP in-line command was built up so that when an atom/ion exceeded a certain kinetic energy threshold (e.g. 10 eV by default) in whichever di- rection, a decelerating force (as the ESP is not a force by itself - see subsection 3.2) corresponding to its energy was subtracted from the force acting upon the particle via other interactions in that direction. Moreover, this was only calculated when the energetic projectile had more than 5 neighbours around it (the amount of atoms on an arbitrary face of a FCC structure) in the region of the target, i.e. it was penetrating the material. This excluded the possibility of ESP calculations in dimers, trimers and bigger sputtered clusters that might be sputtered out of the target. The force calculations depended only on the magnitude of the kinetic energy and the distance travelled by the projectile. New forces were calculated with the help of a fitted function of ESP values obtained from SRIM experimental fit databases [44]: f [i][n] = f [i][n]− sgn(v[i][n]) · 100.217148·log(Ekin[eV ])−1.1887√ (14) f [i][n] ≈ f [i][n]− sgn(v[i][n]) · 0.064759 · Ekin[eV ], (15) 29 where f [i][n] is the force acting on the i-th atom in the n-th direction, sgn(v[i][n]) is a signum function of the velocity of the i-th atom in n-th direction and the constants used are the fit parameters. The latest version of ESP add-on for LAMMPS uses the SRIM lookup tables directly for faster calculations. This means that any element-in-element ESP interac- tion approximations which are supported by SRIM can be used by the ESP LAMMPS add-on. Nevertheless, equation 14 was tested in many situations and was found to produce physically correct approximations for the ESP of Cu-in-Cu in the following simulations. The tutorial on incorporating the latest ESP add-on (written in C++) into LAMMPS and how it can be used can be found in Appendix B. 5.3 Bombardment and sputtering yield Whilst the ion has been mentioned as the experimental projectile, no charge was applied on the projectile atom in the MD setup, as the non-relativistic heavy (compared to electrons and protons) ion was assumed to neutralise almost immediately when coming in contact with the target surface atoms. The initial location of the projectile atom was 5−15 Å above the topmost target layer with a random, yet uniformly distributed, shift from the surface’s centerpoint in the range of half of a unit cell in both x- and y-directions. Thereafter, its total velocity was calculated according to the kinetic energy formula: m · v2ion E ionkin = , (16) √ √ 2 2 · Ekin[J ] 2 · 1.6 · 10−19 · Ekin[eV ] vion = = [m/s], (17) mion[kg] mion[kg] where the Ekin was chosen according to the FGS (see subsection 2.2.2) in the range of 100 to 5000 eV. Doubly charged ions were not considered, as according to the research, their amount was imperceptible. Although preliminary simulations were carried out for the full energy range, the focus was shifted to lower energies of 100-500 eV in the proceeding simulations, as most of the ions that are not yet in the plasma regime lie in the lower part of that energy span and are of practical interest [59]. Since bombardment MD simulations deal with fast moving particles with high kinetic energies, adaptive timestep had to be introduced in order to avoid unphysical behaviour due to too long discrete ”jumps”. Each particle was allowed to move a 30 maximum of 0.1 Å (OM of atomic vibrations [60]) during one timestep, so the duration of the timestep ∆t was reduced accordingly, with the range being ∆t = [0.01, 2] fs. Each bombardment simulation was allowed to run for around 2000-4000 timesteps with 60 up to 200 reruns proportional to the kinetic energy of the projectile. Longer simulation runs were necessary for equalibration of produced collision cascades. At the end of each simulation, neutrals at the distance of z0 = 5 Å above the surface were counted to be sputtered. See an exemplifying simulation in figure 28 in Appendix C. 5.3.1 Incident angle considerations Although the ions were mostly accelerated in normal direction to the surface in the FGS, different incident angles were of interest in this thesis in order to get a wider understanding of the temperature’s interplay with the incident angle and to emulate surface roughness (R ∼ 10nm). The total velocity was distributed between the x-, y- and z-direction components accordingly: vx = vion · sin(θ) · sin(Φ), (18) vy = vion · sin(θ) · cos(Φ), (19) vz = vion · cos(θ), (20) where θ is the incident angle of the ion with the sin(θ) and cos(θ) having non-uniform distributions and Φ is a randomised azimuthal angle in the xy-plane which produces an uniform distribution for its sine and cosine, respectively: Φ(0, 2π) = cos−1(1− 2 · Pu(0, 1)) [rad], (21) where Pu(0, 1) is an uniform distribution from 0 to 1. Equation 21 disproves a common misconception of choosing the random azimuthal angle as Φ = 2π · Pu(0, 1), as this would produce a non-uniform distribution for its sine and cosine [61]. The fixed angles and angle ranges studied with the MD simulations were moti- vated by the subsequent considerations: • 0 deg or (0;3) deg — replicating normal incidence sputtering experiments — as these were found to produce similar results, θ = 0 deg was used by default whenever normal incidence was studied; • 7 deg [46] or 20 deg [62] — reducing the effect of channelling; 31 • (17;23) deg or (27;33) deg — emulating inclined features in the OM of 10 nm; • (0;20) deg — emulating rough surfaces in the OM of 10 nm; • 30, 45, 55, 65, 75, 85 deg — investigating the angle dependence of the SY. For some simulation series, the Penetration Depth (PD) was also logged as the projectile’s final z-coordinate minus average z-coordinate of the topmost layer. When the projectile bounced back from the surface, it produced a negative PD, so these results (< −5 Å) were excluded from average PD calculations. In addition, for an- gular dependence series, Backscattering Ratio (BR) was found as the ratio of recoiled projectiles to the incident projectiles. 5.3.2 Crystallographic orientations and averaged SY The Cu crystallographic faces studied — (100), (110) and (111) — were oriented so that their surface was tangent to the z-direction or the normal bombardment direction (θ = 0 deg, see figure 10 on page 23). The polycrystalline or OFHC Cu was assumed to consist equally of third of each crystallographic orientation, so the average sputtering yield for each incident ion energy was found as a weighed average: Y100(EK) · g100 + Y110(EK) · g110 + Y111(EK) · g111 Y (EKin) = , (22) g100 + g110 + g111 where g is the orientation-dependent (depicted by the subscript) weight function: 1 g(E 2Kin) = ( ) , (23) δY (EKin) where δY (EKin) as the uncertainty or (statistical) standard error of the mean of the average SY was found for each orie√ntation as:∑N 2 δY (E ) = k · i=1(Yi(EKin)− Y (EKin))Kin , (24) N · (N − 1) where k = 1.96 for 95% confidence level, N was the number of simulation in series, Yi(EKin) was the SY of the i-th simulation and Y (EK) was the arithmetic average of the SY over the certain simulation configuration series. The uncertainties and averages for PD and BR values described in the previous subsection were found using the same formulae, replacing SY values for the corresponding PD or BR values. 32 5.3.3 Rough surface considerations Besides studying the SY for atomically smooth surfaces, dependence of rough sur- faces on the SY was of interest, which was motivated by an article by Makeev and Barabási [63]. It stated that ”... ion bombardment of solid targets induces a random (self-affine) morphology on the ion-eroded surfaces”. Moreover, through theoretical calculations it was found that ”... for certain ranges of the parameters variations, sur- face roughness leads to substantial enhancements in the yield, with magnitude of the effect being more than 100%, as compared to the flat surface value. Furthermore, we find that, depending on the interplay between these parameters, the surface roughness can both enhance and suppress the sputter yields.” Rough surfaces used for the Cu self-bombardment MD simulations were generated with a Kinetic Monte Carlo (KMC) simulation code (R ≈ 1nm) by V.Jansson from UH (see figure 12). In order to be certain of the physical authenticity of these simulation cells, they were equilibrated in the same way as the atomically smooth surfaces. Bigger surface roughness (R ≈ 10nm) was emulated by angle ranges as previously mentioned and was not created similarly to small scale roughness simulation boxes as this would have produced a drastic increase in the amount of atoms in the simulation box and therefore the computational time. As the experimental surface profiles in the FGS have not been fully studied, they could not be used in this thesis. Figure 12: Rough surfaces produced by a KMC code, where the different crystallographic faces are clearly visible. Here, z0 denotes the lowest point on the surface. 33 6 Molecular dynamics simulation results 6.1 Cu EAM potential validation First of all, a proper copper (Cu) EAM potential had to be chosen and validated for the following sputtering yield (SY) simulations. All of the Cu EAM potentials files were obtained from the Interatomic Potentials Repository Project [53] and will be referenced by the main author’s name in the succeeding text. The central EAM potential used throughout this thesis was produced by Mishin et al., which has been fitted to the cohesive energy, the bulk modulus, vacancy formation energy, and lattice parameter of Cu. It is considered to be trustworthy as it has been used in more than a thousand research projects according to its citations and is said to be perfect for sputtering simulations [52]. Nonetheless, the pairwise potential part of the EAM potential was splined to the ZBL potential (see subsection 3.3) at r0 = 1.5 Å for a more valid short-distance interaction (see figure 13). This did not alter any significant parameters the potential can produce (see table 3), rather it would become more apparent when the atoms underwent a highly energetic collision (see figure 14 on page 35). Figure 13: Close-up of the Cu-Cu pair potentials of the pristine Mishin potential (blue) and the Mishin potential splined to the ZBL close-interaction potential (green). Before proceeding with the chosen Mishin-ZBL Cu EAM potential (Mishin-ZBL, for short) in the current study, validation of the potential had to be conducted. Since 34 Table 3: Bulk properties of Cu at T = 0K above the Tm: target (experimental) and EAM calculation values [64]. Target values for Mishin were a bit different [52]; as calculated, Mishin-ZBL had the same presented parameters as pristine Mishin. Property Exp. Mishin(-ZBL) Mendelev Lattice constant a (FCC) [Å] 3.640 3.615 3.639 Cohesion energy Ecoh [eV/atom] -3.540 -3.540 -3.283 Melting temperature Tm [K] 1356 1326 1353 most of the experiments in sputtering have been carried out at the normal incidence (θ = 0 deg) of the incoming ions, the validation simulations were also run under the same conditions. Average sputtering yield at 300 Kelvins at normal incidence using different EAM potentials 10 1 100 1000 10000 Energy of the ion [eV] Mishin-ZBL Mendelev Mishin Yamamura semi-empirical function Figure 14: Averaged SY at 300K for normal incidence (θ = 0 deg) using variations of Cu EAM potentials in comparison with SY fit function Y (Ek, θ = 0 deg) by Yamamura et al. In figure 14, it can be seen that, on average, the SY results for Mishin-ZBL lie closer to the semi-empirical Yamamura fit-function for SY (see equation 1) at room temperature than for pristine Mishin or Mendelev potentials. Although Mishin-ZBL overestimates the SY at 100 and 200 eV, the qualitative OM of the results was more important in that energy range, which is why the potential was not exchanged. Furthermore, the evolution of the potential energy per atom was plotted out for the (100) simulation box (as (110) and (111) gave the same results) relaxation at 35 Sputtering yield different temperatures, which is presented in figure 15. It can be well observed, that an increase in the ambient temperature of the target surface reduced the potential energy per atom. For total energy per atom, the graphs were similar with a bit higher values. Potential energy per atom during (100) relaxation at different temperatures -3.33 30 x 30 x 30 simulation cell | Total of 108 000 atoms -3.36 -3.39 -3.42 -3.45 -3.48 -3.51 0 5 10 15 20 25 Elapsed time [ps] 300 Kelvins 500 Kelvins 800 Kelvins 1200 Kelvins Figure 15: Potential energy per atom during (100) relaxation at different temperatures. 6.2 Angular dependence 6.2.1 Normal incidence For normal incidence of the projectile (θ = 0 deg), the SY results for each crystallo- graphic face (including the weighed average) for the Mishin-ZBL validation simulation series introduced in the previous subsection are presented in figure 16 on page 37. In addition, comparison with the Yamamura SY fit function (see equation 1), Eckstein & Preuss SY fit function [32] and SRIM calculation [44] results are included. Firstly, it can be well observed that for normal incidence, the MD results are consistent with the Yamamura and SRIM SY results and visa versa. Per contra, fit function provided by Eckstein & Preuss clearly underestimates the average SY. Secondly, it can be well observed that for the (110) face, the SY had the lowest value at higher energies (2 keV and 5 keV) compared to other orientations, which is explained by the channelling phenomena (see subsection 3.4). On the contrary, the (111) face evinced the biggest SY at higher energies, while the (100) face displayed 36 Potential energy per atom [eV] Figure 16: SY at 300 K for θ = 0 deg. For comparison, results from SRIM and SY fit functions by Yamamura et al. and Eckstein & Preuss (EP) have been provided. Above: full energy range (0.1-5 keV), below: low energy range (0.1-0.5 keV). highest SY in the energy mid-range (0.3-1 keV). For projectile energies less than 1 keV, the (110) and (111) faces produced rather comparable results. To investigate the temperature dependence, 300, 500 and 1200 K simulation series for 100-5000 eV were carried out — these results are presented in figure 17 on page 38. From this figure, it can be seen that the SY varies very little with temperature and 37 some of the differences may even be considered coincidental. In addition, increasing the temperature had no strong effect on channelling in the (110) direction as was expected, although an decrease in the PD with increasing temperature can be witnessed. Figure 17: Above: averaged PD and PD for the (110) face, below: averaged SY and SY for the (110) face at 300/500/1200 K for θ = 0 deg. Note, that the averaged PD was more or less the same for all temperatures in the limits of uncertainty. Moreover, the (100) and (111) faces produced similar SY results as presented in figure 16 and similar PD results as the average PD for all temperatures. 38 6.2.2 Small angle deviations In order to reduce the observed channelling events, two fixed incident angles were considered: θ = 7 deg [46] and θ = 20 deg [62]. As the 300 K θ = 7 deg series displayed rather similar behaviour to figure 16 (although SY for the (110) face did indeed rise), 300 K θ = 20 deg graph is presented for contrast in figure 18. It can be seen, that the (111) face now produces, on average, the least amount of sputtered neutrals, while the (110) face displays no apparent channelling. Sputtering yield at 300 K at 20deg 14 12 10 8 6 4 2 0 100 1000 10000 Energy of the ion [eV] (100) (111) Yam. (0 deg) (110) Average SRIM (20 deg) Figure 18: SY at 300 K for θ = 20 deg for different crystallographic orientations with SY fit by Yamamura et al. (θ = 0 deg) and SRIM SY calculations (θ = 20 deg) for comparison. Note the increase in the average SY with increasing θ. Temperature dependence for full energy range at θ = 7 deg and θ = 20 deg is presented in figure 20 on page 40 for 300 and 1200 K. As it can be observed, the SY slightly increases with temperature for θ = 7 deg, but for θ = 20 deg, the difference is not clearly distinguishable. Moreover, for 500 and 1200 K, the biggest difference was only observed in the case of θ = 7 deg on the (110) face (see figure 19 on page 40). In addition, comparing the penetration depths for θ = 0 deg and θ = 20 deg in figure 21 on page 41, it can be seen that the PD also changes with increasing θ for different crystallographic orientations. On the other hand, on average, the PD seemed to be very insensitive to the temperature of the target and the incident angle of the ion (see figure 29 in Appendix C). 39 Sputtering yield Average sputtering yield for (110) at different temperatures 7 deg 9 8 7 6 5 4 3 2 1 0 100 1000 10000 Energy of the ion [eV] 300 K 500 K 1200 K Figure 19: SY at 300/500/1200 K for θ = 7 deg on the (110) face. For other crystallographic orientations for θ = 7 deg and θ = 20 deg, the results at different temperatures were more similar to one another, i.e. 300 K (see figure 18). Average sputtering yield at different temperatures 7 deg vs 20 deg 12 10 8 6 4 2 0 100 1000 10000 Energy of the ion [eV] 300 K 7 deg 1200 K 7 deg 300 K 20 deg 1200 K 20 deg Figure 20: Averaged SY at 300/1200 K for θ = 7 deg and θ = 20 deg. 6.2.3 Full angle range For a fuller understanding of the angular dependence of the SY, θ was varied for the full energy range at T = 300K — SY results are presented in figure 22 on page 41. 40 Sputtering yield Sputtering yield Average penetration depth at 300 K with 0/20 deg on the (100)/(110)/(111) face 100 10 1 100 1000 10000 Energy of the ion [eV] 300 K (100) 0 deg 300 K (110) 0 deg 300 K (111) 0 deg 300 K (100) 20 deg 300 K (110) 20 deg 300 K (111) 20 deg Figure 21: Penetration Depth (PD) [Å] at 300 K with θ = 0/20 deg on the (100)/(110)/(111) face. The results for the (100) face at lower energies have been excluded as for both incident angles, the PD results were similar to 300 K (111) θ = 20 deg results. Moreover, results for higher temperatures have been ruled out of the figure because no apparent temperature dependence was noticed besides 1200 K (110) in figure 16. Average SY at 300 K at different incident angles or angle ranges 12 10 8 6 4 2 0 100 1000 10000 Energy of the ion [eV] Normal incidence 20 deg 17-23 deg 7 deg 0-20 deg 27-33 deg Figure 22: Averaged SY at 300 K for different angles. Notice the growing trend. Due to practical interest, the dependence was further studied for energies 100- 500 eV in the proceeding simulations (see section 5.3). The incident angle of 100, 200 41 Penetration depth [ang] Sputtering yield and 500 eV ions were varied respectively for both 300 K and 1200 K bombardment simulations (θ = [ 0, 7, 20, 30, 45, 55, 65, 75, 85 ] deg). For these series, the simulation boxes were also increased in the x- and y-directions ((42x42x20) unit cells for the (100) face) in order to avoid damping of the collision cascade caused by energetic atoms reaching the walls of the simulation boxes due to increased θ. The results for SY have been summarised in figure 23. The angle dependent SY fit proposed by Eckstein & Preuss mentioned in subsection 3.1.2 and visualised in figure 7 on page 16 was excluded as it produced even lower SY results than the SY fit by Yamamura et al. or SRIM calculations. Additionally, average backscattering ratios were calculated — these are sum- marised in figure 24 on page 43. The SY and BR results for different crystallographic orientations evinced similar behaviour as the averaged results, so these will not be presented. MASTER SPUTTERING YIELD 5 4.5 4 3.5 3 2.5 2 1.5 1 0.5 0 0 10 20 30 40 50 60 70 80 90 Incident angle of the ion [deg] 300 K 100 eV 300 K 500 eV SRIM 500 eV 1200 K 100 eV 1200 K 500 eV YAM 100 eV 300 K 200 eV SRIM 100 eV YAM 200 eV 1200 K 200 eV SRIM 200 eV YAM 500 eV Figure 23: Averaged SY at 300/1200 K for the 100/200/500 eV ion with various angles. SRIM calculations and angule dependent SY fit by Yamamura et al. added for comparison. 42 Sputtering yield MASTER BACKsCATTER 1 0.9 0.8 0.7 0.6 0.5 0.4 0.3 0.2 0.1 0 -0.1 0 10 20 30 40 50 60 70 80 90 Incident angle of the ion [deg] 300 K 100 eV 300 K 200 eV 300 K 500 eV 1200 K 100 eV 1200 K 200 eV 1200 K 500 eV Figure 24: Averaged BR at 300/1200 K for the 100/200/500 eV ion with various angles. 6.3 Surface roughness dependence The morphology of the rough surfaces used in the following simulations for the crys- tallographic faces (100), (110) and (111) can be seen in figure 12 on page 33. The normal incident angle MD setup was again set up according to the FGS and as it can be seen in figure 25 on page 44, deviations between all of the 300 K results (besides θ = (27, 33) deg) are negligible, while hitting the 1200 K rough surface produced less sputtered neutrals than the atomically smooth surface at the same temperature in the energy range of 200-400 eV. The additional 300 K θ = (0, 20) deg line in this case represents an emulated fairly smooth surface with shallow inclined (up to 20 degrees) features and the 300 K θ = (27, 33) deg line represents an emulated surface consisting only of inclined features (also in the same range as θ). These two emulated R ≥ 10nm with no edges nor apexes. Investigating the effect temperature may have on rough surfaces, it can be seen in figure 26 on page 44, that for 300 K, the KMC generated rough surfaces decreased the SY for all projectile energies in comparison to 1200 K with no substantial deviations. For the orientation dependence for the 300 K θ = (27; 33) deg series, a decrease in the SY was mostly observed for the (110) and (111) faces, while the (100) face produced roughly the same SY for both the flat and the rough surface. 43 Backscattering ratio Average sputtering yield at different temperatures normal incidence 4 3.5 3 2.5 2 1.5 1 0.5 0 100 200 300 400 500 Energy of the ion [eV] 300 K 0 deg smooth 1200 K 0 deg rough 300 K 0 deg rough 300 K (0;20) deg smooth 1200 K 0 deg smooth 300 K (27;33) deg smooth Figure 25: Averaged SY at 300/1200 K: smooth versus rough surfaces. Average sputtering yield at different temperatures with (27:33) deg 4 3.5 3 2.5 2 1.5 1 0.5 100 200 300 400 500 Energy of the ion [eV] 300 K (27;33) deg smooth 1200 K (27;33) deg smooth 300 K (27;33) deg rough 1200 K (27;33) deg rough Figure 26: Averaged SY at 300/1200 K for θ = (27; 33) deg. 44 Sputtering yield Sputtering yield 6.4 Comparison with a previous model Finally, a simulation series for 300, 800, 1300 K, θ = 5 deg, Ekin = 50, 100, 150 eV was carried out in order to compare the results with the MD results presented in [39]. As it can be seen, for 100 and 150 eV , the SY results at 300 K are coinciding, while with increase in the temperature, the results differ drastically. Average SY at 300/800/1300 K with 5 deg 100 10 1 0.1 0.01 200 400 600 800 1000 1200 1400 Temperature of the target [K] 50 eV 150 eV Ins-100eV 100 eV Ins-50 eV Ins-150eV Figure 27: Averaged SY results produced with the MD setup described in this thesis compared with the SY results produced by Insepov et al [39] at 300/800/1300 K for the 50/100/150 eV ion with θ = 5 deg. Notice the log-scale on the y-axis (SY). 45 Sputtering yield 7 Discussion In the first place, the Mishin-ZBL Cu EAM potential was proven to produce physically very sound Sputtering Yield (SY) results for normal incidence (θ = 0 deg) in figure 13 on page 34 in comparison to the widely used SY fit function by Yamamura et al., which is why it was also considered fit for the following simulations. According to the figure, the ZBL potential which was splined to the pristine Mishin potential, gave atoms an additional repulsive force, i.e. more atoms were sputtered. Continuing on with the 300 K normal incidence MD simulation series in figure 16 on page 37, SRIM calculations for the average SY at θ were found to produce compa- rable results with the MD results and Yamamura fit function, while the SY fit function proposed by Eckstein & Preuss underestimated the SY and was excluded from the following simulations. Moreover, the channelling phenomenon was witnessed for the (110) face but only clearly for higher energies such as 2 and 5 keV. In comparison, the other Sputtering Yield (SY) results for projectile energies less than 1 keV were quite comparable with the (100) and (111) faces. This might indicate that there is a certain threshold energy for the channelling event, much higher than thermal energies, which is confirmed by the fact that channelling experiments have been usually carried out at higher energy ranges, in which case the phenomenon is more evident [46]. The bombardment simulations onto the (111) face also evinced a sudden change in the SY in the higher energy range with an almost exponential increase. These observed high-energy effects for the (110) and (111) surfaces are most probably caused by the relative packing density of each crystallographic face (see figure 10 on page 23). As the incident angle of the projectile was increased from θ = 0 deg to θ = 20 deg, the channelling effect previously witnessed for the (110) lost its effect as the projectile’s direction was not in line with the surface’s ”channels” anymore. In contrast, the (111) face began to produce the least amount of sputtered atoms, which might indicate a projectile angle range for which the (planar) channelling probability for this crystallographic face should increase. An article has suggested this angle to be around θ ≈ 31 deg [36], although this was not confirmed within this thesis. In addition, looking at the 300 K Penetration Depth (PD) graph in figure 21 on page 41 and placing it next to figure 16 on page 37 (300K, θ = 0 deg) and figure 18 on page 39 (300K, θ = 20 deg), it can be clearly observed that deeper PD tends to correspond to lower SY. This is why, for normal incidence, the (111) face indicating 46 the lowest and the (110) indicating the highest PD produce the highest and the lowest SY, respectively. Furthermore, focusing on the θ = 20 deg results, it can be seen that the PD decreased with the increase of the incident angle, which correlates to the increase in the average SY. Nonetheless, it has to be noted, that a difference of approximately 1 − 2 Å in the PD results does not have any meaningful essence (see figure 29 in Appendix C) in comparison with the lattice constant λ ≈ 3.615 Å. To begin with the sputtering yield’s temperature dependence, looking at the equilibration graph in figure 15 on page 36, a linear decrease in the potential energy per atom was observed with increasing temperature. This strengthened the speculation that with increasing temperature, it should be easier to kick atoms out of the surface and therefore the SY should increase as well. Accordingly, a linear increase for the SY with increasing temperature was to be expected. On the other hand, the increase in the SY was confirmed only by θ = 7 deg (110) results in figure 19 on page 40, as other configurations in the θ = 7, 20 deg series evinced no discernible change in the SY and figure 20 on page 40, where the increase in the average SY was confirmed for the θ = 7 deg series respectively by comparing the SY for 300 K and 1200 K (∆T = 900K). Moreover, it can be clearly seen in figure 22 on page 41 that, on average, by increasing the incident angle, the SY also tends to increase. The data extracted from these MD simulations is in line with the experimental data [65], which also suggests that this increase in the SY is bound to decrease from beyond some incident angle. The latter was well confirmed with the 100, 200, 500 eV ion bombardment series in figure 23 on page 42 at 300 and 1200 K. Additionally, this shows that for both, small and oblique incident angles, higher target temperatures may increase the SY, especially with the increase of the kinetic energy of the projectile, while for other angles in between (i.e. in figure 17), although the increasing tendency was witnessed, the relative difference stayed much smaller and within the limits of uncertainty (see figure 30 on page 66). The decrease observed for the oblique incident angles can be explained by the surface binding energy of the target, which makes the atom travelling close to the surface, refract towards it. This has also been confirmed theoretically and experimentally. [26] In figure 23 on page 42, it was also confirmed, that SRIM and the angle dependent SY fit function by Yamamura et al. cannot produce comparable results with the MD simulations, at lower projectile energies (Ekin < 500 eV ), as they have been built up on more assumptions and fit constants than the current simulations. 47 As mentioned, the influence of temperature on the SY was expected to be linear, whilst the figure 24 on page 43 designates it to be non-linear. Besides, in figure 19 on page 40, the results for the mid-range temperature of 500 K indicate it might be dependent on the interplay of the energy of the projectile and its incident angle, at least for the 1 keV energy of the projectile. For other energies, the difference between 300 K and 500 K results stays within the uncertainty limits. On other hand, this might also be considered a coincidental anomaly, as the SY results are even more varying with temperature in figures 17 on page 38, although for the latter, the decrease in average PD seems to correlate with increasing temperature. Coming to the surface roughness simulation series, noticeable difference in the SY was expected. This was based on the reasoning that whenever a projectile hits a sharp tip (an apex/a peak) on the surface, with the atoms in it being held in place by less neighbours than in the bulk, those atoms can be hit out of their place with less force (i.e. the SY increases). On the other hand, when a projectile hits a dell, the sputtered atoms (with less kinetic energy than the projectile) may hit inclined features nearby and get ”stuck” (i.e. the SY decreases). It should be clear that depending on the morphology of the surface, the incident angle and the kinetic energy of the projectile, either the former or the latter phenomenon dominates, i.e. an increase or a decrease can be observed in the SY. However, at best, small negative deviations could be observed by introducing roughness into the MD model (see figures 25 and 26 on page 44). Additionally, as the SY is dependent on a complex interplay between the parameters studied before and as the MD morphology model for the FGS system is missing, nothing certain can be said about the influence temperature has on roughness and it in turn on the SY in the experimental set-up. All in all, it has become clear that the SY is not significantly dependent the temperature, on the potential energy per atom and furthermore, that the temperature dependence is most likely not linear. Besides, it can be assumed that even if the smaller PD and/or potential energy per atom for higher temperatures should increase the SY, the thermal vibration energies that increase with temperature help to dissipate the kinetic energy of the projectile as much, annihilating the positive effect. This would explain why the deviations observed in the MD simulations were rather minuscule and at some configurations, even with a negative effect on the SY. This is to be contrasted with the model that was used for comparison in section 6.4. 48 Conclusion and outlooks With the international development of new high-energy particle accelerators, includ- ing the Compact Linear Collider (CLIC), scientists and engineers have had to face many economical, but also physical problems, which have emerged due to pushing the materials used inside the complex to their physical limits. One of the biggest issues that has been holding back the construction of the CLIC is that there are too many Electrical BreakDown (EBD) events (more than 10−6 EBD/pulse/m) taking place inside the accelerating Radio-Frequency (RF) cavities. These components, working at room temperature at electrical field strengths of around 100− 300MV/m, are made of extremely high-purity copper with a surface roughness of R ≈ 10nm. The BreakDown Rate (BDR) for this material should theoretically be much lower, which is why scientists have been looking for reasons, under what circumstances can the copper surface induce and maintain these electric discharges. The most widely spread hypothesis has stated that the EBDs are triggered by the field emission or emission of electrons from sharp tips on the surface by high electric fields, followed by the plasma arc formation. The latter is believed to be fuelled by atoms sputtered out of the surface by ion bombardment. As the field emission and the plasma arc heat the surface, there has been an interest on how the temperature of the surface influences the Sputtering Yield (SY) or the amount of atoms hit outside of the surface by the energetic ions accelerated by the plasma sheath potential. This Master’s thesis was motivated by the controversy between the previous ex- perimental and the computer simulation results - the former have suggested that the temperature has almost no or up to 15% increasing/decreasing effect on the SY, while the simulations have suggested more than a 200% increase. The aim of this thesis was to study the temperature dependence of the SY using the Molecular Dynamics (MD) program LAMMPS in order to offer useful data to the plasma formation model that has been developed in the UH. The following was done: • the electronic stopping power was incorporated into LAMMPS; • the Ziegler-Biersack-Littmark close interaction potential was splined to the widely spread Cu Embedded-Atom Method (EAM) potential by Mishin et al.; • the Mishin-ZBL potential was validated against the semi-empirical SY formula by Yakamura et al. and other Cu EAM potentials (pristine Mishin and Mendelev); 49 • around 70,000 MD simulations at different initial conditions were set up and carried out for statistics to investigate the sputtering yield’s dependence on: – the crystallographic orientation; – the incident angle of the (projectile) ion; – the surface roughness; – the temperature and its interplay with the previous mentioned parameters. The conclusions on the results are the following: • the average potential energy per atom is directly dependent on the temperature of the target: with an increase in the temperature, the average potential energy per atom drops; • the SY for each crystallographic orientation studied ((100), (110) and (111)) is different, especially when θ is small; • channelling phenomena was witnessed for the (110) face for higher energies: the effect was found to decrease with lower projectile energies (Ekin < 2 keV ) and with increasing incident angle θ; • for other crystallographic orientations, no obvious observations about the chan- nelling phenomena were made; • the average SY was found to increase with increasing θ with the maximum SY produced around θmax ≈ 45 deg for the 200 and 500 eV ions), after which the SY begins to decrease again; • the SY angular dependence simulations were found to produce more accurate results at lower projectile energies than the SY fit function by Yamamura et al. or SRIM calculations; • the Backscattering Ratio was found to increase with increasing θ and decrease with increasing energy of the projectile, being fairly dependent on the tempera- ture; • the mean Penetration Depth (PD) of the projectile was qualitatively found to be in negative correlation with the SY, but rather for higher energies; 50 • the temperature dependence of the SY is minuscule and most likely not linear: the increasing effects can be annihilated by e.g. thermal vibrations, which dissipate the kinetic energy of the projectile; • the surface morphology can either increase or decrease the SY; • the increase in temperature did not evince any apparent addition to the effect morphology had on the SY; All in all, no major anomalous deviation in the sputtering yield was observed with varying the crystallographic orientation, the incident angle of the ion, the surface roughness and most of all, the temperature of the target surface. Besides using the knowledge gained with this thesisto advance the CLIC research, it might prove itself useful in other high electric field applications as well, such as Free-Electron Lasers (FEL) or medical accelerators. In addition, it was found that some MD simulations have previously overesti- mated the effect temperature has on the SY by up to more than 1000%, confirmed by no experimental considerations. In the light of this, the author is planning to set up additional MD simulations to establish a wider understanding of the processes concern- ing the temperature dependence of the SY. The parameters to be studied will include the dimensions of the collision cascades and the tensile stress induced by the electric field, as this is present in the RF cavities and the FGS. Moreover, the author and the research groups working in UT and UH are interested in the influence grain bound- aries, copper oxide and other additives may have on the SY, which remain far out of the scope of this thesis and shows that there is yet a lot to be investigated. 51 CLIC’i elektriliste läbilöökide uuringud: vase enesepihustuse molekulaardünaamika simulatsioonid 0.1-5 keV ioonide korral kõrgendatud temperatuuridel Tarvo Metspalu Kokkuvõte (Conclusion in Estonian) Seoses kaasaegsete kõrgenergiliste osakeste kiirendite, k.a. Kompaktse Lineaarpõrguti (ingl. k. Compact Linear Collider ’i ehk CLIC ’i), kontseptsioonide rahvusvahelise arendustööga, peavad teadlased ja insenerid ületama majanduslike probleemide kõrval ka füüsikalisi komplikatsioone, kuna kõrgtehnoloogilised seadmed seavad kasutatavatele materjalidele kõrgeid nõudmisi. CLICi valmimisel on üheks suurimaks takistuseks olnud raadiosageduslikes õõsre- sonaatorites (ingl. k. radio-frequency cavity) asetleidvate elektriliste läbilöökide (ingl. k. electrical breakdown e. EBD) liiga kõrge sagedus (ingl. k. breakdown rate e. BDR) — rohkem kui 10−6 läbilööki ühe impulsi kohta kiirendi meetri kohta. Need komponendid on tehtud erakordselt puhtast oksiidivabast vasest (OFHC — puhtusega vähemalt 99,999%), mille pinnakaredus jääb suurusjärku R ≈ 10nm ning need peaksid töötama toatemperatuuril kõrges 100 − 300MV/m elektriväljas. Teoreetiliselt peaks OFHC läbilöögisagedus olema palju madalam, mistõttu otsitakse nii eksperimentaalselt kui ka teoreetiliselt, millistel tingimustel on vasepind võimeline taolisi elektriläbilööke tekitama ja hoidma. Kõige laiemalt tunnustatud hüpoteesi kohaselt saavad elektrilised läbilöögid al- guse kõrgetest elektriväljadest tingitud väljaemissioonist ehk elektronide emissioonist, mis võimendub pinnal asetsevates teravates tippudes. Arvatakse, et algusetapil hoiavad tekkivat kaarlahendust üleval ioonpommitamise tagajärjel pinnast välja löödud aato- mid, misjuhul pommitavad plasma kattepotentsiaalis kiirendatud vaseioonid vasepinda ehk leiab aset vase enesepihustus. Kuna väljaemissioon ning kaarlahendus kuumutavad vasepinda, on teadlaste seas tekkinud suur huvi, kuidas kõrgendatud temperatuur mõjutab aatomsaagist (ingl. k. sputtering yield e. SY ) ehk välja löödud aatomite hulka ühe pommitava iooni kohta. Käesolev magistritöö teema, ”CLIC’i elektriliste läbilöökide uuringud: vase ene- sepihustuse molekulaardünaamika simulatsioonid 0.1-5 keV ioonide korral kõrgendatud 52 temperatuuridel”, oli ajendatud vastuolulistest eksperimentaalsetest ning arvutisimu- latsioonide tulemustest. Esimesed neist leidsid, et temperatuuri tõusuga võib aatom- saagis jääda muutumatuks või muutuda kuni 15% (erinevates teadustöödes kas tõusvas või langevas joones), siis simulatsioonid on hinnanud kuni 200% tõusu. Antud magistritöö eesmärgiks oli uurida vask-vask ioonpommitamise ehk vase enesepihustuse aatomsaagise temperatuurisõltuvust kasutades molekulaardünaamika (MD) programmi LAMMPS, et pakkuda kasulikku teavet Helsingi ülikoolis välja tööta- tud teoreetilisele plasmatekke mudelile. Töö käigus tehti järgnevat: • LAMMPS-i jaoks kirjutati ning lisati elektroonse pidurdusjõu (ingl. k. electronic stopping power ehk ESP) täiendus; • laialdaselt levinud Mishin et al. poolt välja töötatud vase varjestatud aatomi meetodi (ingl. k. embedded atom method ’i e. EAM ’i) potentsiaalile (täpsemalt paaripotentsiaalile) sobitati Ziegler-Biersack-Littmark’i lähimõju potentsiaal; • Mishin-ZBL EAM potentsiaal valideeriti võrdluses Yakamura et al. poolt välja töötatud poolempiirilise aatomsaagise valemiga ning algelise Mishin et al. ja Mendelev et al. poolt välja töötatud EAM potentsiaaliga; • seati üles ning viidi läbi ligi 70000 molekulaardünaamika simulatsiooni erinevatel algtingimustel, et uurida aatomsaagise sõltuvust: – kristallograafilisest tasanditest (100), (110) ja (111); – iooni langemisnurgast θ; – pinnakaredusest; – temperatuurist ja selle koosmõjust erinevate eelmainitud teguritega. Vase ioonpommitamise molekulaardünaamika simulatsioonidest saadud statistika analüüsimisel tehti järgnevad järeldused: • keskmine potentsiaalne energia aatomi kohta on otseselt sõltuv sihtmärgi (vase pinna) temperatuurist: temperatuuri kasvades keskmine potentsiaalne energia aatomi kohta väheneb; • aatomsaagis on erinevate kristallograafiliste suundade jaoks erinev, eriti madalate θ väärtuste korral; 53 • (110) kristallograafilise pinna puhul täheldati kõrgemate iooni energiate korral (Ekin > 2 keV ) kanaldumise efekti, mis vähenes iooni langemisnurga kasvades ning energia vähendamisel; • teiste kristallograafiliste pindade puhul ((100) ja (111)) kanaldumise efekti ei täheldatud; • ioonpommitamise aatomsaagis kasvas iooni langemisnurga θ kasvades, saavu- tades maksimaalse väärtuse iooni energiast sõltuva nurga korral, misjuhul edasise langemisnurga kasvatamisel aatomsaagis taas kahaneb; • keskmine iooni läbitungimissügavus (ingl. k. penetration depth e. PD) kor- relleerub kvalitatiivsel tasemel negatiivselt aatomsaagisega; • ioonpommitamise aatomsaagise temperatuurisõltuvus on tühine ja tõenäoliselt mittelineaarne: aatomsaagist tõstvad efektid tasakaalustatakse näiteks soojus- võnkumiste poolt, mis hajutavad temperatuuri tõustes põrkekaskaadi suurema efektiivsusega; • pinnamorfoloogia võib ioonpommitamise aatomsaagist tõsta või väheneda, kuid see sõltub suuresti konkreetsest pinnast; • sihtmärgi temperatuuri tõstmine ei ilmutanud nähtavat lisa pinnamorfoloogia poolt ioonpommitamise aatomsaagisele osutatavale mõjule. Kokkuvõttes ei täheldatud kristallograafiliste orientatsioonide, iooni langemis- nurga, pinnakareduse ega ka sihtmärgi temperatuuri varieerides ioonpommitamise aatom- saagise puhul ühtegi suurt anomaalset kõrvalkallet, mis on kooskõlas varasemate eksper- imentaalsete tulemustega. Magistritöö raames leiti, et mõningad varasemad MD simulatsioonid on temper- atuuri mõju aatomsaagisele hinnanud antud tööst rohkem kui 1000 % suuremaks, toetu- mata sealjuures eksperimentaalsetele kaalutlustele. Seetõttu plaanib käesoleva töö au- tor üles seada täiendavad MD simulatsioonid, et saada ioonpommitamise aatomsaagise temperatuurisõltuvusest sügavam arusaam. Selleks uuritakse tulevikus läbiviidavates arvutisimulatsioonides elektrivälja poolt pinnale ostutatava rõhu, kristallograafiliste piirpindade, muude kohalike defektide, vaskoksiidi ja teiste lisandite mõju ioonpommi- tamise aatomsaagisele, mis jäid oma mahukuse tõttu antud magistritööst välja ning ilmestavad, et antud teema kohta on veel palju uurida. 54 Acknowledgements First of all, I would like to thank my inspiring supervisors. I would like to thank Vahur Zadin for accepting me into the international CLIC project by taking me under his supervision with my Master’s thesis and also for the recommendations he made on my work. I had the chance to visit our colleagues in the Accelerator Laboratory in University of Helsinki. I would like to thank Flyura Djurabekova for welcoming me in Finland, for helping me get started in the world of Molecular Dynamics and for giving me very specific tasks, always keeping a hand on the pulse of my progress. I would especially like to thank Konstantin Avkhachev for sharing me his Python scripts for splining and plotting ZBL, introducing me into MD, LAMMPS, Shell script- ing, C++, Git and of course for the unprecedented contribution he made in upgrading the newest ESP add-on for LAMMPS. In addition, I would like to thank the co-workers in UH and UT for their tips and cookies, particularly Ville Jansson for the rough surface KMC models; Mihkel Veske and Simon Vigonski for helping me set up my computations; Stefan Parvilainen and Kristian Kuppart for consultations on MD. Finally, I would like to thank Elen for helping me get through the hard times. 55 References [1] CERN, “Cern homepage: Basic information about the large hadron collider (lhc),” 2015. [2] A. Levy, “Clic detector and physics study overview.” 2015. [3] CERN, “Cern homepage: Basic information about the clic project: explanations, collaborators, goals,” 2015. [4] H. Timko, K. N. Sjobak, L. Mether, S. Calatroni, F. Djurabekova, K. Matyash, K. Nordlund, and R. Schneider, “From field emission to vacuum arc ignition : a new tool for simulating copper vacuum arcs,” Contributions to Plasma Physics, pp. 1–17, 2014. [5] A. F. Voter, “Introduction to the kinetic monte carlo method,” Radiation Effects, vol. 235, p. 1–23, 2005. [6] D. Brandt, “Cern accelerator school rf for accelerators,” tech. rep., 2006. [7] D. Hopkins, A. Sessler, and J. Wurtele, “The two-beam accelerator,” Applied Physics Letters, 1984. [8] J. Raguin, D. Schulte, I. Syratchev, I. Wilson, and W. Wuensch, “A new damped and tapered accelerating structure for clic,” Proceedings of LINAC2002, pp. 307–309, 2002. [9] I. Wilson, “The compact linear collider,” Tech. Rep. December, 2004. [10] W. Wuensch, “Introduction to rf acceleration: Basic rf concepts,” 2013. [11] H. Timko and K. Nordlund, “Arcing and plasma-wall interactions,” vol. 64, 2008. [12] F. Tecker, “Clic progress report at the 8th international accelerator school for linear colliders,” tech. rep., 2013. [13] D. Halliday, R. Resnick, and WalkerJearl, Fundamentals of Physics. John Wiley & Sons, 9th ed., 2013. [14] S. Matsumoto, M. Akemoto, S. Fukuda, T. Higo, N. Kudoh, H. Matsushita, H. Nakajima, T. Shi- dara, K. Yokoyama, and M. Yoshida, “The status of nextef: the x-band test facility in kek,” Proceedings of LINAC08, pp. 906–908, 2008. [15] H. Timko, Modelling vacuum arcs: from plasma initiation to surface interactions. PhD thesis, 2011. [16] M. Kildemo, S. Calatroni, and M. Taborelli, “Breakdown and field emission conditioning of cu, mo, and w,” Physical Review Special Topics - Accelerators and Beams, vol. 7, no. 9, pp. 82–88, 2004. [17] A. Moilanen, Creep effects in diffusion bonding of oxygen-free copper. PhD thesis, 2013. [18] N. C. Shipman, Experimental study of DC vacuum breakdown and application to high-gradient accelerating structures for CLIC. PhD thesis, 2014. [19] H. Timko, F. Djurabekova, K. Nordlund, L. Costelle, K. Matyash, R. Schneider, a. Toerklep, G. Arnau-Izquierdo, a. Descoeudres, S. Calatroni, M. Taborelli, and W. Wuensch, “Mechanism of surface modification in the plasma-surface interaction in electrical arcs,” Physical Review B - Condensed Matter and Materials Physics, vol. 81, pp. 1–8, 2010. [20] P. A. Anderson, “The work function of copper,” Physical Review, vol. 76, no. 3, pp. 388–390, 1949. [21] A. P. Fontenla, “Post-mortem analysis : Sem imaging review,” tech. rep., 2015. [22] P. Sigmund, “Recollections of fifty years with sputtering,” Thin Solid Films, vol. 520, pp. 6031– 6049, 2012. 56 [23] E. A. Kenik and T. E. Mitchell, “Orientation dependence of the threshold displacement energy in copper and vanadium,” Philosophical Magazine, vol. 32, no. 4, pp. 815–831, 1975. [24] Y. Yamamura, Y. Itikawa, and N. Itoh, “Angular dependence of sputtering yields of monatomic solids,” 1983. [25] B. Delley, D. Ellis, A. Freeman, E. Baerends, and D. Post, “Binding energy and electronic structure of small copper particles,” Physical Review B, vol. 27, no. 4, pp. 2132–2144, 1983. [26] R. Behrisch and W. Eckstein, “Sputtering by particle bombardment: Experiments and computer calculations from theshold to mev energies,” in Topics in Applied Physics (R. Behrisch and W. Eckstein, eds.), vol. 110, pp. 1–526, Springer, 2007. [27] R. S. Averback and T. D. de la Rubia, Displacement Damage in Irradiated Metals and Semicon- ductors, vol. 51. 1997. [28] E. Bringa, K. Nordlund, and J. Keinonen, “Cratering-energy regimes: From linear collision cascades to heat spikes to macroscopic impacts,” Physical Review B, vol. 64, no. 23, pp. 1–12, 2001. [29] T. Diaz de la Rubia, R. S. Averback, R. Benedek, and W. E. King, “Role of thermal spikes in energetic collision cascades,” Phys. Rev. Lett., vol. 59, no. 001, pp. 1930–1933, 1987. [30] E. Y. Zykova, a. S. Mosunov, J. S. Colligon, and V. E. Yurasova, “Mass dependence of sputtering characteristics for a number of metals with different properties,” Vacuum, vol. 85, no. 4, pp. 546– 549, 2010. [31] P. Sigmund, Fundamental Processes in Sputterin g of Atoms and Molecules. 1993. [32] W. Eckstein and R. Preuss, “New fit formulae for the sputtering yield,” Journal of Nuclear Materials, vol. 320, no. 3, pp. 209–213, 2003. [33] Y. Yamamura and H. Tawara, “Energy dependence of ion-induced sputtering yields from monatomic solids at normal incidence,” Atomic Data and Nuclear Data Tables, vol. 62, no. 0005, pp. 149–253, 1996. [34] C. Garćıa-Rosales, W. Eckstein, and J. Roth, “Revised formulae for sputtering data,” Journal of Nuclear Materials, vol. 218, no. 1, pp. 8–17, 1995. [35] K. I. Grais, A. A. Fakhry, D. E. Atta, R. M. Boutros, A. E.-a. M. S, and J. E. Ibrahim, “A new model describes the simultaneous effect of the incident ion energy and angle,” Journal of Applied Sciences Research, vol. 9, no. 6, pp. 3890–3904, 2013. [36] J. J. Ph. Elich, H. E. Roosendaal, and D. Onderdeltnden, “Copper single crystal sputtering in the temperature range from 50k to 600k,” Radiation Effects, vol. 10, no. 3, pp. 175–184, 1971. [37] C. E. Carlston, G. D. Magnuson, A. Comeaux, and P. Mahadevan, “Effect of elevated tempera- tures on sputtering yields,” Physical Review, vol. 138, no. 3A, pp. 759–763, 1965. [38] P. Sigmund and M. Szymonski, “Temperature-dependent sputtering of metals and insulators,” Applied Physics A Solids and Surfaces, vol. 34, no. 4, p. 247, 1984. [39] Z. Insepov, J. Norem, and S. Veitzer, “Atomistic self-sputtering mechanisms of rf breakdown in high-gradient linacs,” Nuclear Instruments and Methods in Physics Research, Section B: Beam Interactions with Materials and Atoms, vol. 268, no. 6, pp. 642–650, 2010. [40] C. Staudt, R. Heinrich, P. Mazarov, A. Wucher, V. I. Tugushev, and N. K. Dzhemilev, “On the temperature dependence of sputtered cluster yields,” Nuclear Instruments and Methods in Physics Research B, vol. 164, no. 165, pp. 715–719, 2000. [41] J. Schou, “Slowing-down processes , energy deposition , sputtering and desorption in ion and electron interactions with solids .,” no. 1, pp. 1–10. 57 [42] H. Paul, “Nuclear stopping power and its impact on the determination of electronic stopping power,” AIP Conference Proceedings, vol. 1525, pp. 309–313, 2013. [43] D. Duffy and A. Rutherford, “Including the effects of electronic stopping and electron-ion inter- actions in radiation damage simulations,” J. Phys.: Condens. Matter, vol. 19, pp. 1–11, 2007. [44] J. F. Ziegler, M. D. Ziegler, and J. P. Biersack, “Srim - the stopping and range of ions in matter (2010),” Nuclear Instruments and Methods in Physics Research, Section B: Beam Interactions with Materials and Atoms, vol. 268, no. 11-12, pp. 1818–1823, 2010. [45] M. Baker and H. Cohen, “An explanation of miller indices introduction : Miller indices,” tech. rep., 2004. [46] D. S. Gemmell, “Channelling and related effects in the motion of charged particles through crystals,” in Rev. Modern Physics, vol. 46, pp. 1–107, 1974. [47] K. E. Jensen, D. Pennachio, D. Recht, D. a. Weitz, and F. Spaepen, “Rapid growth of large, defect-free colloidal crystals,” Soft Matter, no. di, pp. 320–328, 2013. [48] W. Cai, J. Li, and S. Yip, “Molecular dynamics,” in Comprehensive Nuclear Materials (R. Stoller, ed.), ch. 128 Molecu, pp. 1–36, 2010. [49] D. Brandell, Understanding Ionic Conductivity in Crystalline Polymer Electrolytes. PhD thesis, 2005. [50] M. S. Daw, S. M. Foiles, and M. I. Baskes, “The embedded-atom method: a review of theory and applications,” Materials Science Reports, vol. 9, no. 7-8, pp. 251–310, 1993. [51] “Lammps documentation: Eam potential.” [52] Y. Mishin, M. Mehl, D. Papaconstantopoulos, a. Voter, and J. Kress, “Structural stability and lattice defects in copper: Ab initio, tight-binding, and embedded-atom calculations,” Physical Review B, vol. 63, pp. 1–16, 2001. [53] N. M. M. Laboratory, “Interatomic potentials repository project,” 2015. [54] S. M. Foiles, “Comparison of binary collision approximation and molecular dynamics for dis- placement cascades in gaas,” Tech. Rep. September, 2011. [55] L. Bukonte, F. Djurabekova, J. Samela, K. Nordlund, S. a. Norris, and M. J. Aziz, “Comparison of molecular dynamics and binary collision approximation simulations for atom displacement analysis,” Nuclear Instruments and Methods in Physics Research, Section B: Beam Interactions with Materials and Atoms, vol. 297, pp. 23–28, 2013. [56] M. T. Robinson, “The binary collision approximation: Background and introduction,” Radiation Effects and Defects in Solids, vol. null, no. 1, pp. 3–20, 1994. [57] S. N. Laboratories, “Lammps molecular dynamics simulator by sandia,” 2015. [58] S. Bringuier, “Molecular dynamics primer : Lammps examples,” 2014. [59] H. Timko, K. Ness Sjobak, L. Mether, S. Calatroni, F. Djurabekova, K. Matyash, K. Nordlund, R. Schneider, and W. Wuensch, “From field emission to vacuum arc ignition: A new tool for simulating copper vacuum arcs,” Contributions to Plasma Physics, vol. 55, no. 4, pp. 299–314, 2015. [60] L. Cartz, “Thermal vibrations of atoms in cubic crystals ii: The amplitude of atomic vibrations,” Proceedings of the Physical Society. Section B, vol. 68, no. 11, pp. 957–967, 2002. [61] E. W. Weisstein, “Sphere point picking.” 58 [62] L. S. Yu, J. M. E. Harper, J. J. Cuomo, and D. a. Smith, “Alignment of thin films by glancing angle ion bombardment during deposition,” Applied Physics Letters, vol. 47, no. 9, pp. 932–933, 1985. [63] M. a. Makeev and A. L. Barabási, “Effect of surface morphology on the sputtering yields. i. ion sputtering from self-affine surfaces,” Nuclear Instruments and Methods in Physics Research, Section B: Beam Interactions with Materials and Atoms, vol. 222, no. 3-4, pp. 316–334, 2004. [64] M. Mendelev, M. Kramer, C. Becker, and M. Asta, “Analysis of semi-empirical interatomic potentials appropriate for simulation of crystalline and liquid al and cu,” Philosophical Magazine, vol. 88, no. 12, pp. 1723–1750, 2008. [65] J. M. Fluit, P. K. Rol, and J. Kistemaker, “Angular-dependent sputtering of copper single crys- tals,” Journal of Applied Physics, vol. 34, no. 3, pp. 690–691, 1963. 59 Appendix A: Source codes The following simplified LAMMPS input scripts mostly used were initiated by a BASH code which is not presented here. It provided a more user-friendly interface, ran the simulations using OpenMPI, logged the results in a well-mannered hierarchy, plotted necessary graphs per each simulation series and varied a random seed number and counter for each bombardment simulation. The scripts were kept apart for conve- nience and clarity, although they could be merged into one script. For rough surface simulations, slightly different scripts were used, also not presented here. As internally, there is just one random number generator for all equal-style vari- ables and one for all atom-style variables, by using the random() or normal() math functions, the internal random number generators will only be initialized once, which means only one of the specified seeds will determine the sequence of generated ran- dom numbers. By varying the seed value, different results are achieved and moreover, knowing the seed value, it is possible to reproduce the specific simulation results. First of all, the equilibraion script is presented, which was initiated with the fol- lowing predefined variables: orientation - 100, 110 or 111, POT - name of potential used, surftemp - temperature of the target, relax - length of the relaxation phase in timesteps and (Lx, Ly, Lz) - dimensions of the simulation box in (100) unit cells. # --- INITIALIZATION --- units metal dimension 3 boundary p p p atom_style atomic # --- VARIABLES --- variable latparam equal 3.615 # Assumed lattice parameter # --- Simulation box size --- if "${orientation} == 100" then & "variable sizeX equal v_Lx" & "variable sizeY equal v_Ly" & "variable sizeZ equal v_Lz" & elif "${orientation} == 110" & "variable sizeX equal ’v_Lx ’" & "variable sizeY equal ’ceil(v_Ly/sqrt (2))’" & "variable sizeZ equal ’ceil(v_Lz/sqrt (2))’" & elif "${orientation} == 111" & "variable sizeX equal ’ceil(v_Lx/sqrt (8/3))’" & "variable sizeY equal ’ceil(v_Ly/sqrt (2))’" & "variable sizeZ equal ’ceil(v_Lz/sqrt (3))’" & else "print ’Unknown orientation !’" "quit" # Simulation box size dependent variables variable Z_ion equal v_sizeZ *5/4 # Room for the ion # --- ATOM DEFINITION & REGIONS --- # Orienting according to the crystallographic orientation 60 if "${orientation} == 100" then & "lattice fcc ${latparam} orient x 1 0 0 orient y 0 1 0 orient z 0 0 1" & elif "${orientation} == 110" & "lattice fcc ${latparam} orient x 0 0 -1 orient y -1 1 0 orient z 1 1 0" & elif "${orientation} == 111" & "lattice fcc ${latparam} orient x 1 1 -2 orient y -1 1 0 orient z 1 1 1" & else "print ’Unknown orientation !’" "quit" # Defining the region that is going to be bombarded region regioon block 0 ${sizeX} 0 ${sizeY} 0 ${sizeZ} create_box 1 regioon create_atoms 1 region regioon group pind region regioon # --- EAM POTENTIAL --- pair_style eam/alloy pair_coeff * * Potentials/Cu.eam.Mishin.zbl # --- EQUILIBRATION AND SETTINGS --- compute Pot all pe/atom compute Kin all ke/atom # Initial velocities/temperature velocity pind create ${surf_temp} 482748 mom yes rot no # Giving the surface a temperature fix 1 pind npt temp ${surf_temp} ${surf_temp} 1 iso 0 0 1 reset_timestep 0 timestep 0.002 # Timestep of 1 to 5 femtoseconds thermo 50 # Printing out information after every ... timesteps # Displayed parameters thermo_style custom step cpu vol lz press pxx pyy pzz pe etotal temp dt run ${relax} # Relax for ... timesteps # --- RELAXATION 2 --- change_box all z delta -1.5 ${Z_ion} run ${relax} # Relax for ... timesteps write_dump all custom Relax/${POT}/${surf_temp}Kelvins/RELAX.orient${ orientation }.${surf_temp}K.$(( v_relax *2))TS.DUMP id type x y z vx vy vz fx fy fz variable NRA equal count(pind) # --- END OF SCRIPT --- Secondly, the bombardment script was initiated with the following variables in addition to the ones used in the equilibration script: seed - random integer (uniform distribution), count - number of the simulation run in the specific simulation series, energy - energy of the ion in eV and (low, high) - incident angle range in degrees. The first part of the script used the same structure as the equilibration script from the initialization until the simulation box size dependent variables and continued as follows: # --- SIMULATION BOX SIZE DEPENDENT VARIABLES --- variable sizeX2 equal v_sizeX /2 # Middle of initial box size X variable sizeY2 equal v_sizeY /2 # Middle of initial box size Y variable sizeZ2 equal v_sizeZ /2 # Middle of initial box size Z variable X_1 equal v_sizeX -1 variable Y_1 equal v_sizeY -1 # Defining simulation box region regioon block 0 ${sizeX} 0 ${sizeY} 0 ${sizeZ} 61 create_box 1 regioon # --- ION SPEED --- # Ion energy in eV -s (divided by 100 because of the default A/ps unit) variable m_Cu equal (63.546*(1/12) *1.998467052*10^ -26) variable iSpeed equal (sqrt (2* v_energy *1.602176565*10^ -19/( v_m_Cu)))/100 # --- EAM POTENTIAL --- pair_style eam/alloy # Potential type pair_coeff * * Potentials/${POT} Cu # Potential file # --- READING IN PRE -RELAXED STATE --- read_dump Relax/${POT}/${surf_temp}Kelvins/RELAX.orient${orientation }.$ {surf_temp}K.$(( v_relax *2))TS.DUMP $(( v_relax *2)) x y z vx vy vz box yes add yes # --- RANDOM VARIABLES: COORDINATES , ANGLES , VELOCITY COMPONENTS --- # Random X value in the middle of the simulation box (+-1/2 lattice) variable RanX equal (random(v_sizeX2 -1/2, v_sizeX2 +1/2, v_seed)) variable RanY equal (random(v_sizeY2 -1/2, v_sizeY2 +1/2, v_seed)) variable RandomX equal ${RanX} # Fixing the random value variable RandomY equal ${RanY} # Random incident angle variable RanAlpha equal random(${low},${high},v_seed) variable RandomAlpha equal ${RanAlpha} variable RanPhi equal (acos (1-2* random(0,1,v_seed))*360/ PI) variable RandomPhi equal ${RanPhi} variable iSpeedZ equal v_iSpeed*cos(v_RandomAlpha*PI/180) variable iSpeedXY equal sqrt(( v_iSpeed*v_iSpeed) -(v_iSpeedZ*v_iSpeedZ)) variable iSpeedX equal sin(${RandomPhi })*${iSpeedXY} variable iSpeedY equal cos(${RandomPhi })*${iSpeedXY} # --- RANDOM BOX SHIFTING --- variable RanDX equal (random(0,v_sizeX ,v_seed)) variable RanDY equal (random(0,v_sizeY ,v_seed)) displace_atoms all move ${RanDX} ${RanDY} 0 # --- ION CREATION --- create_atoms 1 single ${RandomX} ${RandomY} ${Z_ion} # Creating one ion region iooniregioon sphere ${RandomX} ${RandomY} ${Z_ion} 0.1 group ion region iooniregioon # Defining a ion region # --- NEW REGIONS --- region atomregion block INF INF INF INF INF INF group atoms region atomregion region bottomregioon block 0 ${sizeX} 0 ${sizeY} INF 2 group bottom region bottomregioon region anttiwalls block 1 ${X_1} 1 ${Y_1} -1 INF group antiwalls region anttiwalls group walls subtract atoms antiwalls # Finding the group of atoms that is going to be integrated over time group integraalgrupp1 subtract atoms bottom group integral union integraalgrupp1 ion # Bombarded target atoms + ion - bottom - walls group bomb_surf subtract integral walls # Bottom is fixed , top is shrink -wrapped change_box all boundary p p fs # --- SETTING ION VELOCITY --- velocity ion set ${iSpeedX} ${iSpeedY} -${iSpeedZ} mom yes rot no units box # --- RUN SETTINGS --- compute IonKin ion ke compute IonZ ion reduce max z 62 # Displayed parameters thermo_style custom step cpu dt time ke c_IonKin c_IonZ reset_timestep 0 fix adaptive_timestep all dt/reset 10 1.0e-5 0.002 0.1 units box fix 3 bomb_surf nve fix thermowalls walls nvt temp ${surf_temp} ${surf_temp} 1 neigh_modify check yes every 1 delay 0 # Check for building the regular neighbor lists more frequently to avoid errors # --- ELECTRONIC STOPPING POWER --- fix STOP all elstop 10 region regioon # --- BALANCING PROCESSOR LOAD --- comm_style tiled fix BAL all balance 200 1.05 rcb # ---RUN --- thermo 25 # Show parameters every ... timesteps # Dump location of every atom after every 10 timesteps. dump 3 all custom 10 BOMBDUMP.dump id x y z if "${energy} == 100" then & "variable bombrun equal 2000" & elif "${energy} == 200" & "variable bombrun equal 2250" & elif "${energy} == 500" & "variable bombrun equal 2500" & elif "${energy} == 1000" & "variable bombrun equal 3000" & elif "${energy} == 2000" & "variable bombrun equal 3500" & elif "${energy} == 5000" & "variable bombrun equal 4000" & else "variable bombrun equal 2000" run ${bombrun} unfix adaptive_timestep unfix 3 unfix thermowalls unfix STOP unfix BAL undump 3 # --- PENETRATION DEPTH --- variable Z_final equal c_IonZ variable Z_surf equal ’v_sizeZ * v_latparamZ ’ variable Depth equal ’v_Z_surf - v_Z_final ’ # --- SPUTTERING YIELD --- region sputterbox block INF INF INF INF ${Z_ion} INF group sputtergroup region sputterbox variable NRA equal count(sputtergroup) variable NRAtot equal count(atoms) # Logging the results of the run print "${POT} ${orientation} ${count} ${seed} ${surf_temp} ${energy } ${RandomAlpha} ${RandomPhi} ${RandomX} ${RandomY} ${NRAtot} ${Depth} ${ NRA}" append Logs/Sputter.log screen no # --- END OF SCRIPT --- 63 Appendix B: Incorporating ESP into LAMMPS In order to incorporate the Electronic Stopping Power (ESP) add-on into LAMMPS: 1. download the latest ESP fix-style files from https://www.ims.ut.ee/mediawiki/upload/6/67/LAMMPS_elstop.zip; 2. unpack the contents into LAMMPS/src folder; 3. in the same folder, open style fix.h and add #include ”fix elstop.h” to the list; 4. recompile LAMMPS (check different possibilities in LAMMPS/src/MAKE folder). The ESP fix-style command has the following syntax: fix ID group-ID elstop cutoff LUT style args • ID = user-assigned name for the fix • group-ID = ID of the group of atoms to perform the computation on • elstop = style name of this fix command • cutoff = cutoff energy from which the ESP is applied in eV-s • LUT = Look-Up Table file, from where the ESP values are extracted • style (optional) = region (region args = region-ID) The LUT, which can be obtained from the õ program [44], must have the following structure (each column separated with a tab; with exemplifying values): Energy of the atom in eV-s M1-in-M1 in eV/Å M2-in-M1 in eV/Å ... 9.99999 7.364E-01 2.794E-01 ... 10.9999 7.723E-01 2.930E-01 ... ... ... ... ... Here M1, M2, ... denote different materials or atom types. M1-in-M2 — atoms of the first element moving inside the target composed of the second element atoms — these atom types must be predefined inside the script. The command can be applied multiple times. Finally, here is an example of using it iside a LAMMPS script: fix STOP all elstop 10 ../scripts/elstop.PtCoGa region surface2 64 Appendix C: Additional figures Figure 28: An example of the MD setup cross-section with visualized kinetic energies using OVITO (Open VIsualization TOol): Ek = 1 keV projectile, θ = 3 deg on the (100) face (SY = 3). 65 Average penetration depth at 300/500/1200 K with 0/20 deg 10 1 100 1000 Energy of the ion [eV] 300 K 0 deg 500 K 0 deg 1200 K 0 deg 300 K 20 deg 500 K 20 deg 1200 K 20 deg Figure 29: Average PD for 300/500/1200 K for θ = 0 deg and θ = 20 deg. Notice the linear fashion. Average sputtering yield at different temperatures 0 deg vs (27;33) deg 4 3.5 3 2.5 2 1.5 1 0.5 0 100 200 300 400 500 Energy of the ion [eV] 300 K 0 deg 1200 K 0 deg 300 K (27:33) deg 1200 K (27;33) deg Figure 30: Averaged SY at 300/1200 K for θ = 0 deg and θ = (27; 33) deg. 66 Penetration depth [ang] Sputtering yield Non-exclusive licence to reproduce and publicize thesis I, Tarvo Metspalu (15.03.1990), 1. herewith grant the University of Tartu a free permit (non-exclusive licence) to: 1.1. reproduce, for the purpose of preservation and making available to the pub- lic, including for addition to the DSpace digital archives until expiry of the term of validity of the copyright, and 1.2. make available to the public via the university’s web environment, including via the DSpace digital archives, as of 28.05.2018, until expiry of the term of validity of the copyright, “CLIC electrical breakdown studies: Cu self-sputtering MD simulations for 0.1-5 keV ions at elevated temperatures”, supervised by Vahur Zadin, Flyura Djurabekova and Konstantin Avchaciov. 2. I am aware of the fact that the author retains these rights. 3. This is to certify that granting the non-exclusive licence does not infringe the intel- lectual property rights or rights arising from the Personal Data Protection Act. Tartu, 28.05.2015