Simulation of how pressure influences the reactive sputtering process

UPTEC F 15035

Examensarbete 30 hp

Juni 2015

Simulation of how pressure influences

the reactive sputtering process

Erik Strandberg

Teknisk- naturvetenskaplig fakultet UTH-enheten

Besöksadress:

Ångströmlaboratoriet Lägerhyddsvägen 1 Hus 4, Plan 0

Postadress:

Box 536 751 21 Uppsala

Telefon:

018 – 471 30 03

Telefax:

018 – 471 30 00

Hemsida:

http://www.teknat.uu.se/student

Abstract

Simulation of how pressure influences the reactive

sputtering process

Erik Strandberg

Sputtering is a physical vapor deposition (PVD) process used to create thin films, i.e very thin layers of material. To form compounds, such as oxides and nitrides, it may be beneficial to add a reactive gas to the process which is known as reactive sputtering. This thesis focuses on the simulation of the reactive sputtering process and, more specifically, the effect of the process pressure.

Two models have been developed. A Monte Carlo model simulates the distribution of sputtered material throughout the chamber. It is based on the binary collision model with initial conditions acquired from simulations in TRIM. The hard-sphere potential is used as interaction potential in the scattering calculations. The effect of the process pressure is studied for two different elements, sulfur and tungsten. It is found that the distribution of material is heavily influenced by the pressure. A high pressure gives a more diffusion-like distribution compared to a low pressure. As the pressure is increased the deposited material’s energy distribution is found to be shifted towards lower energies until it reaches the energy of thermalized atoms.

The second model developed is an extended Berg model that incorporates the effect of redeposition of sputtered material on the target, implantation of reactive ions in the target and preferential sputtering. Using simulations the effect of these extensions is discussed. It is found that an increased pressure may eliminate the hysteresis region which has been observed experimentally. Finally an outline is presented on how the two models can be unified into a Berg-model that takes the non-uniform distribution of sputtered material into account.

Ämnesgranskare: Sören Berg Handledare: Erik Särhammar

Populärvetenskaplig sammanfattning

Tunnfilmer, väldigt tunna lager av material, har en stor spridning i vår moderna värld. Exempel  på användningsområden inkluderar metallskiktet i CD­skivor, transistorer i elektriska kretsar,  solceller och hårda beläggningar på borrar. För att skapa väldigt tunna lager finns det flera  tekniker varav den som studerats i detta examensarbete heter sputtering. Vid sputtering 

placeras det objekt som skall beläggas (substratet) tillsammans med materialet som skall utgöra  beläggningen (targeten) i en vakuumkammare. Genom att beskjuta targeten med energirika  partiklar så kommer det att “slås ut” atomer från targeten. Dessa kommer att färdas genom  kammaren för att fastna på substratet och därigenom skapa tunnfilmen. För vissa beläggningar,  exempelvis oxider, är det även önskvärt att tillsätta en reaktiv gas, såsom syre, till kammaren i  en process som kallas reaktiv sputtering. 

För att modellera reaktiv sputtering finns det olika modeller varav en är Berg­modellen. 

Berg­modellen inkluderar dock inte processtrycket som i experiment visat sig påverka  processen. Det här examensarbetet har därför fokuserat på att utöka modelleringen för att  inkluderar processtrycket. Utökningen består av två separata modeller som utvecklats och  implementerats som simuleringsprogram. 

Den första modellen beräknar hur atomer som lämnar targeten sprider sig i kammaren. Ungefär  tio miljoner partiklar simuleras individuellt. Varje partikel följs från det att den lämnar targeten  tills dess att den träffar någon yta i kammaren. Partiklen kommer under sin färd att krocka med  gasmolekyler i kammaren vilket medför att partikelns riktning och energi förändras. När alla  partiklar har simulerats så aggregeras informationen och resultatet blir en distribution över hur  partiklarna sprider sig i kammaren. Simuleringar visar att processtrycket spelar en stor roll för  spridningen där ett högt tryck medför att partiklarna sprider sig mera än vid ett lågt tryck. Vidare  medför ett högt tryck att partiklarnas energi är lägre då de träffar kammarytan. 

Den andra modellen är en utökning av Berg­modellen med två fysikaliska effekter som är  beroende på processtrycket. Det är dels joniserad reaktiv gas som inplanteras under ytan på  targetet samt även atomer som, på grund av kollissioner med gasmolekyler, återdeponeras på  targeten. Utifrån simuleringar har effekten av dessa mekanismer studerats. Resultaten visar att  den utökade modellen kan förutsäga de trender som observerats experimentellt. Slutligen  presenteras hur de två modellerna kan kombineras till en Berg­modellen som även tar  spridningen av atomerna i kammaren i beaktande. 

Contents

1. Introduction​ ………1 

2. Sputtering process​……… 3 

3. Material distribution​……….. 6 

3.1 Initial conditions​……… 6 

3.2 Scattering throughout the gas phase​……… 11 

3.3 Aggregating the results​………18 

3.4 Interaction potentials​……… 19 

3.5 Simulations and analysis​……… 22 

4. Reactive sputtering simulation​……… 26 

4.1 Berg model​……… 26 

4.2 Dissociative sputtering​……… 29 

4.3 Ion implantation​……… 32 

4.4 Redeposition​………. 34 

4.5 Non­uniform distribution of sputtered material​ ………35 

5. Conclusion​………. 37 

Acknowledgements​………... 38 

References​………. 39   

1. Introduction

Thin film technology is an important part of our modern lifestyle. Without it there would be no  CDs or electronic devices such as computers. So what is a thin film? A thin film is a very thin  layer of material. A thin film may be applied to give a surface certain characteristics that the bulk  material does not exhibit. Examples include extremely hard coatings on drill bits, which enables  to drill into hard materials, or low friction coatings in valves, to achieve better performance in  engines. Or the product may be so small that thin films are part of the product itself. Examples  include the transistors found in modern electronics and thin film solar cells. There are many  different techniques for creating a thin film and most of them belong to one of two major groups,  physical vapor deposition and chemical deposition. [1] 

There are a wide range of different chemical processes with different advantages and  drawbacks. What they do have in common is that they are based on chemical reaction. 

Processes that uses a precursor gas that reacts with the surface of the object to coat are known  as chemical vapor deposition (CVD) processes. CVD processes generally have good uniformity  of the deposited film, i.e the thickness of the thin film is almost the same regardless of the  structure of the coated surface. Another advantage of the CVD processes is that they have a  high throughput which is critical in an industrial application. However CVD processes usually  involve adding energy, in the form of heat, in order to stimulate the chemical reaction. This  makes it hard to grow thin films on materials that are sensitive to heating, such as plastics. The  chemicals involved in a CVD process may be toxic, requiring extra precautions to be taken in  production, and they can also be very expensive. [2] 

The other category of thin films deposition processes are physical vapor deposition (PVD)  processes. In a PVD process atoms are removed from the source material by either vaporizing it  or bombarding it with energetic particles. The former is known as an evaporation process and  the latter as a sputtering process. Both processes requires that the system is operated in a  low­pressure environment. In contrast to the chemical processes the physical deposition  processes are characterised by the directionality of the deposition [2]. It is hence hard to  achieve a uniform film on a substrate with a complex geometry such as holes. However most  coatings can be produced by PVD processes, the deposition rate is good and films can be  grown on heat sensitive substrates [1]. Examples of thin films that are produced by physical  deposition processes are the metallic layer in a compact disc and the diffusion barrier on the  inside of plastic food containers. 

This paper focuses on the sputtering process and, more specifically, on reactive sputtering. 

Reactive sputtering is a process where a reactive gas, such as O​2​, is added to the sputtering  chamber. The gas will react with the deposited metal atoms and aid in the formation of  compounds such as oxides. To get a better insight and to be able to predict the behaviour of  these complex processes it is important to have a model that described them. A well used  model for reactive sputtering is Berg’s model [3]. It is a simple, yet powerful, model that 

described the process at steady state. One property that is omitted in Berg’s model is the partial 

pressure of the sputtering gas. There are however experiments that indicate that it influences  the sputtering process in a non­neglectable way. Hence the purpose of this paper is to extend  Berg’s model so that it incorporates the partial pressure of the sputtering gas and this is done in  two parts. 

In the first part of this paper a new model, based on the binary collision model, is presented [4]. 

The model describes how the sputtered material is distributed throughout the chamber. Using  simulations the effects of different process parameters, such as the pressure and choice of  material, are studied. The second part of the paper describes the enhancements made to Berg’s  model. The physical mechanisms that have been considered are the effect of dissociative  sputtering, redeposition of sputtered material on the target and ion implantation. Results from  simulations are presented as well as an analysis of how these mechanisms influence the  sputtering process. Finally an outline is presented on how the two models can be unified into a  Berg model that takes the non­uniform distribution of sputtered material into account. The  implementation as a simulation software was however considered to be outside the scope of  this project due to the limited amount of time available. Hence no simulation, or analysis, for the  unified model is presented in this paper. 

2. Sputtering process

There are many different variations of the sputtering process but they all share the same  principle. The material to deposit, the target, and the object that is to be deposited, the 

substrate, are put in a sputtering chamber. The system is operated at a low pressure, usually in  the range of a few pascal, and hence the chamber must be sealed and connected to a pumping  system. To grow the film the target is bombarded by energetic particles. The energetic particle  will start a cascade of collisions between the atoms in the target. During these collisions some  of the target atoms will be given such an kinetic energy that they are ejected from the target i.e  they are sputtered. The average number of sputtered particles per incident energetic particle is  defined as the sputtering yield. The sputtered atoms will then traverse the vacuum chamber and  stick to the substrate where the film is grown. [1] 

To generate the energetic particles the sputtering gas is ionized and the ions are accelerated by  an electric field. To achieve this there are different techniques, direct current (DC) sputtering  and radio frequency (RF) sputtering. In the case of DC sputtering a DC power supply is  connected to the chamber in such a way that the target is held at a negative voltage, in the  order of 1 kV, and the chamber is grounded. The substrate may also be grounded or it can be  under floating potential. Due to the difference in potential between the target and the chamber  there is an electric field. The free charge carriers (electrons and ions) will be accelerated by it  and, if accelerated for a long enough distance without colliding, gain a kinetic energy that is  greater than the ionization energy of the sputtering gas. The ions, due to their much larger size,  have a higher probability of colliding and they will thus reach a significantly lower kinetic energy  compared to the electrons. Hence the ionization of the sputtering gas is mainly accounted for by  the electrons. When ionizing the sputtering gas an electron is broken free from its atom and the  number of free electrons will increase. Since the electrons are negatively charged they will  always be accelerated towards the chamber wall. So if an electron is formed close to the target  number of electrons are increased as it is accelerated towards the chamber wall but the newly  created electrons will be formed closer and closer to the chamber wall. To be able to sustain the  ionization process additional electrons must be created close to the target. This is achieved by  the ions, accelerated towards the target, since they will have a great chance of producing new  free electrons upon the impact with the target. Due to this effect the ionization process can be  self­sustained. [1] 

RF sputtering works in a similar way to DC sputtering but due to the very fast oscillation of the  power supply, typically in the range of 10​^7​ Hz, the electric field in the chamber will quickly  alternate. The ions, due to their heavy mass, are not affected by this to any greater extent but  the electrons are. They will travel back and forth ionizing sputtering gas atoms in the process. 

This reduces the need of the free electrons created when the ions collide with the target’s  surface. As a consequence it is possible to sputtering insulating materials with RF sputtering  which is not possible with DC sputtering. This however comes at the cost of a decreased  deposition rate. [1]  

In conventional DC and RF Sputtering the formation rate of ions is rather low resulting in a low  deposition rate. To be able to sustain the ionization process a rather high voltage must be  applied and the process pressure must be kept high as well. In order to mitigate these effects  magnets can be placed behind the target and this is known as magnetron sputtering. The  magnets will result in a strong parallel magnetic field close to the target’s surface. Due to the  Lorentz force the electrons emitted from the target’s surface will be “stuck” in the magnetic field  contributing to the ionization of sputtering gas atoms close to the target’s surface. This results in  a higher density plasma that is closer to the target. A drawback of magnetron sputtering is that  the sputtering will not be uniformly across the target’s surface but there will be areas with a  higher erosion rate. This results in a non­flat surface and since the flux of sputtered particles are  highest at the angle perpendicular to the target’s surface a larger fraction of the sputtered  material will not be deposited on the substrate but on the chamber walls. This results in a less  efficient use of target material. [1] 

When the film to grow consist of a compound (oxides, nitrides etc) a reactive gas, such as O​2​ or  N​2​, is usually added to the chamber. The gas will react with the deposited metal atoms and form  compound. Depending on the partial pressure of the reactive gas more or less compound is  formed. For a very high partial pressure the film will most likely be stoichiometric and for a low  pressure it will not. However the dependence of the composition of the deposited film on the  partial pressure of the reactive gas is not trivial. And the partial pressure of the reactive gas  does not only affect the composition of the deposited thin film. If the partial pressure is high  enough then the target will be mostly compounded resulting in a significantly lower deposition  rate. To make matter worse it is hard to control the partial pressure of the reactive gas in the  chamber in practice. A more convenient process parameter to adjust is the flow of reactive gas  that is let into the chamber. However the dependence between the flow of the reactive gas and  the composition of the thin film is even more complicated than that between the composition  and the partial pressure and it may even included hysteresis effect. [1] 

Figure 1; Hysteresis behaviour with marked transitions (dashed) in a simulated deposition rate 

Hysteresis effects are not unheard of and occurs in several different fields of physics such as  ferroelectric materials [5]. The deposition rate versus the flow of reactive gas for a typical 

sputtering process is shown in Figure 1. Assume that the process is operating stably at point M​1  but it is sought to operate at point P. If the flow of reactive gas is increased, even so slightly, the  process will however not follow the curve to P but instead go directly to C​1​. If the process is to  return to M​1​ it is not sufficient to simply decrease the reactive gas flow to the previous setting. 

Instead the flow must be decreased so the process passes C​0​ at which point it will go directly to  M​o​. It is hence impossible to operate in stable mode in the part of the process curve that lies  between M​1​ and C​0​. This is known as the hysteresis effect in the context of reactive sputtering. 

[1] 

3. Material distribution

To gain further insight into how the particles involved in the sputtering process are distributed  across the chamber a Monte­Carlo simulation model was developed. A Monte­Carlo process is  a process that includes randomness, i.e two consecutive executions of the process will not  generate the exact same result. However the result of the two executions should only differ by  an ‘statistical error’ [6]. The model developed in this paper supports simulating sputtered 

particles, reflected sputtering gas and negative ions formed close to the target. The simulation is  carried out by simulating a large number of particles independently and then aggregating the  results. Each particle is followed until it hits a chamber wall and when it does the particle’s  energy, incident angle and position is saved. Hence each particle undergoes three major steps  during the simulation; obtaining its initial condition, scattering through the gas phase and finally  hitting the chamber wall. 

3.1 Initial conditions

For each sputtered particle its initial position, direction and energy must be defined. The position  is assumed to be uniformly distributed across the target’s racetrack. This is a simplification since  a typical racetrack profile, as seen in Figure 2, shows that there are more particles being 

sputtered from the center of this circular track than from the edges of it. However the effect of  this simplification is minor since the width of the racetrack is small in relation to the dimensions  of the chamber. Supporting a non­uniform distribution would require a more complex function  and hence more processing power resulting in less particles being simulated per unit of time.  

Figure 2 (top); A typical racetrack shown in  profile​^. 

Figure 3 (left); Model of racetrack with inner  radius   and outer radius  .r₀ r₁  

Let the racetrack be modelled as a circle of radius   with a circular hole of radius  . This isr₁ r₀   illustrated in Figure 3. To calculate the initial position of a sputtered particle the following  equation is used 

x = r cos β  

y = r sin β (3.1) 

where  is the azimuthal angle and   the radial distance from the center of the target. Since theβ r   azimuthal angle is uniformly distributed on the interval [ 20, π[ it can be calculated using the  equation 

  πR

β = 2 ₀₁ (3.2) 

where R₀₁ is uniformly distributed random number on the interval [ 10, [. The radial distance   isr   however not uniformly distributed on the interval [r , r ]_{0   1} . To see this consider the area of the  segment bounded by two circles centered at ( 00, ) and with radiuses ρ − 0.5dρ and ρ + 0.5dρ  where dρ is infinitesimal. The area dA of this segment is given by 

A (ρ .5dρ) (ρ .5dρ) πρdρ

d = π + 0 ²− π − 0 ²= 2 (3.3) 

The total area of the racetrack is given as   

r r

A = π ₁²− π ₀²    (3.4) 

Define the cumulative area fraction function CAF F (r )′  as the fraction of the total area of the  racetrack that is bounded by a circle of radius  . By using (3.3) and (3.4) r′ CAF F (r )′  can be  determined using the following function 

AF F (r ) πρ dρ

C ′ =_{πr −πr}¹

12 02

∫

r₀2 =_{r −r}^{r −r′2 02}

12 02 (3.5) 

To randomly pick a value of   on the interval r [r , ]₀ r₁  such that all points on the racetrack have  an equal probability of being picked is equivalent to solving the equation 

AF F (r)

C = R^′₀₁ (3.6) 

where R^′₀₁ is a random number on the interval [  0, 1]. This since points randomly spread across  the racetrack will have a uniformly distributed amount of area on its inside. Substituting (3.5) into  (3.6) yields 

r −r² ₀²

r −r₁² ₀² = R^′₀₁    

which can be rewritten as   

r =

√

To calculate the initial position of the sputtered particle (3.7) and (3.2) are now substituted into  (3.1) which gives 

x =

√

√

where the values of R′₀₁ and R₀₁ are randomized once respectively and then used in both  equations. 

In addition to an initial position each particle is also given an initial direction and energy. These  depends on the type of particle that is being simulated. Negative ions are assumed to form  close to the target surface. They are then immediately accelerated by the electric field and  hence obtains an energy corresponding to the target potential and a direction that is 

perpendicular to the target’s surface. The particle will now be a small distance away from the  target but in order to simplify the simulation model the particle is assumed to be infinitely close  to the target’s surface. It is also assumed that the particle will not undergo any collisions while  being accelerated. The effect of this simplification should be pretty minor since the chance of a  particle colliding is proportional to the distance that it travels during the acceleration which is  small. 

For sputtered and reflected particles the sputtering angle and initial energy are obtained from  TRIM simulations. TRIM (Transport of Ions in Matter) is a software used to simulate the effect of  ion bombardment [7]. It is based on the binary collision model where a flat surface of a user  defined composition is bombarded with ions. The ions have a predefined energy and they hit the  surface at a specified incident angle. The output of the simulation is, among other, a list of the  particles that are sputtered from the surface along with their energy and direction. Since the  sputtering angle and initial energy are correlated they must be saved together and hence the  result of the TRIM simulation is a list of pairs. To obtain the initial sputtering angle and energy of  a particle a pair from the TRIM simulation is picked randomly. The initial direction and energy,  as reported by TRIM, are those of the particle when it is at an infinitesimal distance away from  the target’s surface. However to leave the surface the particle will lose energy equal to the  surface binding energyE_sb. So the energy of the sputtered particle is calculated using the  equation 

E_sput= E_trim− E_sb   

where E_sput is the energy of the particle being sputtered and E_trim is the energy as reported by  TRIM. The energy loss will also affect the sputtering angle. Lets start by noting that the velocity 

 of the particle is parallel to its direction   and if   is normalized it follows that

vˉ Dˉ Dˉ  

vˉ = | ˉv|D (3.8) 

Combining equation (3.8) with the equation for kinetic energy   

E =^mv₂²    (3.9) 

yields the equation   

vˉ =

√

In TRIM the coordinate system is defined as the z­axis being perpendicular to the target surface  and the xy­plane coincides with the target surface. Using (3.10) the velocity, before leaving the  surface, along the z­axis is given by 

D  

v_trim,z=

√

Using (3.9) and (3.11) the kinetic energy, in the direction of the z­axis, before leaving the  surface is given by 

D

E_trim,z=^mvtrim,z₂ ² = Etrim trim,z2 (3.12) 

Since the surface binding energy is lost in the direction that is perpendicular to the surface the  kinetic energy of the particle ​after​ leaving the surface, E_sput,z, is calculated as 

D

E_sput,z= E_trim,z− E_sb= Etrim trim,z2− E_sb    (3.13) 

where (3.12) has been used. Using (3.9) and (3.13) the velocity of the particle in the direction of  the z­axis ​after​ leaving the surface, v_sput,z, is given by 

  v_sput,z=

√

√

sb (3.14) 

The velocity of the particle after leaving the surface, v_sput is given by   

v_sput=

√

√

Let the sputtering angle   be defined as the angle between the particle’s velocity vector, afterα   leaving the surface, and the the z­axis. It can be calculated using the equation 

cos α = _v

sput

v_sput,z

(3.16)   

By substituting (3.14) and (3.15) into (3.16) the equation    

  cos α =  

√

sb    (3.17) 

is obtained. Figure 4 shows the initial energy distribution for tungsten and sulfur along with the  theoretical distribution for a target potential of 400 V [4]. It can be seen that the values from  TRIM match the theoretical distribution very well. The sputtering angle distribution is shown in  Figure 5 for the same materials and target potential. This distribution also match the theoretical  value very well. However it should be noted that these calculations are made under the 

assumption that the target is a perfectly flat surface. In practice this is not the case since the  erosions is not uniform and hence the formation of a racetrack.  

  Figure 4; The energy distribution for sulfur  and tungsten. 

  Figure 5; The sputtering angle distribution for  sulfur and tungsten. 

To get the particle’s initial direction the sputtering angle   must be combined with an azimuthalα   angle  . The azimuthal angle is uniformly distributed on the interval β [  0, 2π[ and given by the  function 

  πR

β = 2 ₀₁ (3.18) 

where R₀₁ is a random number on the interval [  0, 1[. If the target’s surface is placed at the  xy­plane, i.e its normal vector is ( 0 10, , ), then the definition of spherical coordinates 

x = r sin α cos β   y = r sin α sin β   z = r cos α    

together with (3.17) and (3.18) yields the following function to calculate the normalized initial  direction of the particle 

πR   x = sin

(

√

sb

)

πR   y = sin

(

√

sb

)

z = 

√

sb 

3.2 Scattering throughout the gas phase

Once the initial conditions of the sputtered particle are determined the simulation enters its next  step which is to simulate how the particle is scattered throughout the gas phase. This is done by  repeating the following steps for each particle: 

1. Calculate the distance that the particle  will travel before its next potential  collision with a gas particle. The  distance is based on the particle’s  mean free path and a random number. 

2. Check if particle will hit a surface  before reaching the point of its next  collision. If so, calculate the incident  angle and save this together with the  energy of the particle. Move on to the  next particle. 

3. Determine the amount that the particle  is scattered and the energy that it will  lose due to the collision. Using this  information the particle’s new direction  and energy is calculated. 

A schematic overview of the particle’s  traversal is shown in Figure 6. 

  Figure 6; The traversal of a sputtered  particle (red) through the chamber where it  collides with gas particles (blue). 

Distance between collisions

The distance that a particle will travel between two successive collisions is distributed  non­uniformly.  To obtain a value   from this distribution the equationλ  

−

λ = λ₀ln R₀₁ (3.19) 

is used where λ₀ is the mean free path and R₀₁is a random number on the interval ] 10, ][4]. If  the gas is thermalized and the velocity of the sputtered particle is much larger than that of the  gas molecules then the mean free path is given by 

λ₀=_nσ¹   (3.20) 

where   is the concentration of gas molecules and   is the scattering cross­section of a gasn σ   molecule. If on the other hand the velocity of the sputtered particle and the gas molecules are  roughly the same, i.e the sputtered particle has been thermalized, then the mean free path is  given by 

λ₀=_√2¹_nσ   

The scattering cross­section can be seen as the area that a particle must be inside to scatter  against another particle. Its value depends upon what type of potential that is used to simulate  the scattering. The simulation uses the hard­sphere potential so the scattering cross­section is  given by 

(r )

σ = π _gas+ r_particle ² (3.21) 

where r_gas and r_particle are the radiuses of the gas molecule and sputtered particle respectively  [4].  

Using the ideal gas law   

V k T

P = N _B    

the concentration of gas molecules   can be calculated asn    

n = ^N_V =_{k T}^P

B (3.22) 

where  is the gas pressure, P k_B is Boltzmann’s constant and   is the temperature of the gas.T  

−

λ = _{cπP(r +r}^{k TB} ₎

gas particle2ln R₀₁ (3.23) 

where c = √2 when the energy of the sputtered particle is low, i.e the particle is thermalized, and   when the energy is high. The molecule’s energy in a gas is distributed according to the

c = 1  

Maxwell­Boltzmann distribution. The mean energy E_avg for a gas molecule is thus given by [8] 

k T

Eavg=³_{2 B} (3.24) 

From (3.24) it follows that the mean energy of a molecule in a thermalized gas with a 

temperature of 300K is about 0.26 eV. Based on this the sputtered particle is said to have a high  energy, thus using c = 1, if its energy is above 1 eV. 

Collision with wall

Let P be the point where the prior collision occurred or, if the particle has not yet collided, the  position on the target where it was sputtered. Let P’ be the position where the next collision will  occur. Let Fˉ(t)= (F (t), F (t), F (t))_{x   y   z} define the particle’s trajectory from P to P’.   is given by theFˉ   functions 

(t)   (P )t

F_x = P_x′ +   _x− P_x′  

(t)   (P )t

F_y = P_y′ +   _y− P_y′  

(t)   (P )t

Fz = Pz′ +   z− Pz′    

where   is a parameter on the interval t [ 1  0, ].  

The chamber is defined by a number of surfaces that are either bounded or unbounded. An  unbounded surface is a surface that extends infinitely and is only defined by its equation. A  bounded surface on the other hand is defined by its equation but also by at least one constraint  that defines the boundary. Consider one of the surfaces that make up chamber. Let the surface  be described by the equation S. If there exists a solution to the equation 

(t)  S

Fˉ =   (3.25) 

such that   is on the interval t [ 10, ] it means that the particle’s trajectory intersects the surface. 

Solutions that yield a value of t outside the interval [0, 1] are rejected since the intersection is  outside the particle’s trajectory i.e they start before P or end after P’. Using the value of   thet   position of the impact is calculated. If the surface is unbounded the particle will collide 

regardless of the point of intersection. If the surface on other hand is bounded it must be  checked whether the point lies on inside or outside of the surface’s boundary. If the point lies  within the boundary then it will collide, otherwise it will not. It is possible that the particle’s 

trajectory will intersect with more than one surface that defines the chamber. If that’s the case  then the solution with the lowest value of   is the surface that the particle will hit since it is thet   first surface that the particle will intersect with on its path from P to P’. 

If the particle does collide with a surface then the particle’s incident angle is calculated. The  incident angle   is defined as the positive angle between the surface’s normal and the particle’sβ   trajectory, i.e an angle of 0 indicates that the particle will hit the surface perpendicular to the  surface. Let the surface normal be Nˉ = (N ,_z N_y,N_z) and the normalized vector describing the  direction of the particle by Dˉ = (D ,_x D_y,D_z). Using the dot product definition [9] 

|A|| ||B||  (α) Aˉ ∙ Bˉ = | cos    

the incident angle can be calculated as   

  β = arccos

(

)

(

N D +N D +N Dx x y y z z

)

Using the fact that both   and   are normalized, i.e have a length of 1, equation (3.26) can beNˉ Dˉ   rewritten as 

β = π − arccos (N Dx x+ Ny yD + Nz zD )   

which is used to calculate the incident angle. The calculated incident angle is then saved  together with the particle’s energy and the point of the impact. 

Scattering against a gas particle

Since the energy of the gas molecule is distributed according to the Maxwell­Boltzmann  distribution it is given an initial energy of  

k T

E =³_{2 B} (3.27) 

where k_B is Boltzmann’s constant and   is the temperature of the gas. The reasons to notT   assume that the gas particle is at rest is that the energy of the scattered particle would then  approach zero for systems where it would undergo multiple scattering collisions. This since a  scattering collision will always decrease the energy of the scattered particle since some of it will  be transferred to the gas atom. In practice the sputtered particle will, after undergoing multiple  collisions, be thermalized and have an energy according to the Maxwell­Boltzmann distribution. 

Using the equation for kinetic energy and (3.27) the velocity of the gas molecule v_gas is  calculated as 

vgas=

√

The direction dˉ = (_gas dgas,x, d  _gas,y, d  _gas,z) of the gas molecule is calculated using   

dgax,x= sin α cos β  

dgas,y= sin α sin β (3.29) 

d_gas,z= cos α    

where   is randomly picked on the uniform interval β [  0, 2π] and   is randomly picked on theα   uniform interval [  0, π]. Scattering calculation is performed in a xy­plane where the gas molecule  is at rest at origo prior to the collision. However since the gas molecule is moving lets introduce  a new cartesian coordinate system that is fixed in the gas molecule, i.e it moves with the velocity 

 that is calculated as

v_coordˉ  

v v_coordˉ = dˉ_{gas gas}   

The velocity of the sputtered particle in the new coordinate system is calculated using (3.28)  and (3.29) together with the equation 

v

v_partˉ = v′_partˉ − v_coordˉ = v_′partˉ − dˉgas _gas (3.30)   

where v′_partˉ  is the particle's velocity in the initial frame of reference, v_partˉ  is its velocity in the new  coordinate system.  Since there exists infinitely many planes that contains v′_partˉ  there exists an  infinite amount of possible scattering planes. To uniquely fix the scattering plane in the initial  frame of reference an azimuthal angle   is randomly picked on the uniform interval ω [  0, 2π]. It is  now possible to convert the velocity of the particle, and thus its energy and direction, between  the scattering plane and the initial frame of reference.  

To calculate the scattering angle   the sputtered particle is initially placed at ϕ ( L  − , p) where p   is the impact parameter and   is sufficiently large so that there exists no significant interactionL   between the sputtered particle and the gas molecule at this distance. The impact parameter is  uniformly distributed on the interval [  0, p_max]  where p_max is the maximum distance at which  there is any significant scattering. This is typically defined as the cutoff where the scattering is  less than a certain threshold such as 0.01π. The reason for this definition is that there are  interaction potentials that interact at an infinitely long distance and hence there is no value of 

such that the scattering equals zero. An impact parameter of   is thus a head on collision

p_max 0  

and an impact parameter of p_max represents a collision where the particles barely touch each  other. To choose an impact parameter the following equation is used 

R p = p_{max 01} 

where R₀₁ is a random number on the interval [  0, 1]. Lets define   asA    

A =_m^mgas

part (3.31) 

where mgas is the mass of the gas molecule and m_part is the mass of the sputtered particle. 

When calculating the scattering angle a centre­of­mass coordinate system is usually introduced  to simplify the calculation. To transform the particle’s energy into this new centre­of­mass  system the following equation is used 

E

Ecom=_{1+A part}^A     (3.32) 

where E_part is the energy of the particle [4]. By using the equation for kinetic energy (3.32) can  be rewritten as 

v Ecom=^m_1+A^gas _part²   

where v_part is calculated using equation (3.30). If the collision is assumed to be inelastic then  the scattering angle ϕcom in the centre­of­mass system is given by the equation 

ϕcom= π − 2^∞

∫

R 1

r²

√

where V(r) is the interaction potential as discussed in Section 3.4 [4]. For the hard sphere  potential (3.33) reduces to 

ϕcom= π − 2 arcsin_p^p

max    (3.34) 

To transform ϕcom into the scattering angle  the following equation is usedϕ    

ϕ = arctan

(

)

where it is assumed, as before, that the collision is inelastic [4]. Using (3.35), (3.34) and (3.31)   can be calculated as function of the impact parameter  . In Figure 7 a plot of this is shown

ϕ p  

for the collision between an argon gas atom and a sulfur, a tungsten and an argon atom. The  elements are chosen such that the mass of the scattered particle is greater than, less than and  equal to that of the gas particle. From Figure 7 it can be seen that element that have a mass  greater than that of the gas particle will always scatter less than  . More generally it can be^π₂   seen that the scattering angle decreases with increased mass where a heavy elements such as 

atom the scattering angle ranges from   to  . In other words it is possible that a particle lighter0 π   than the gas particle will completely backscatter which is not possible for elements heavier than  the gas particle. 

Figure 7; The scattering angle for sulfur (red),  tungsten (blue) and argon (green) when  scattered against an argon atom. 

  Figure 8; The energy loss for sulfur (red),  tungsten (blue) and argon (green) when  scattered against an argon atom. 

The energy of the sputtered particle after it has been scattered is calculated as   

  E_after= 1

(

)

where E_after is the energy of the particle after it’s been scattered [4]. A plot of the energy loss of  a scattered particle as function of the impact parameter   is found in Figure 8. The mostp   important conclusion that can be drawn from this plot is that elements that are heavier than the  gas atom will never lose all of it’s energy in the scattering whereas elements that are lighter can. 

Also if one correlates the plots in Figure 7 and Figure 8 it can be seen that tungsten will not lose  its maximal amount of energy when it scattered the most. Instead the highest scattering angle  corresponds to a relatively low amount of energy being lost. 

3.3 Aggregating the results

The result of the simulation so far is a list of  impact points along with the incident angle  and the particle’s energy at the time of  impact. To be able to interpret this as a  meaningful result the list must be 

aggregated. First each surface that makes  up a part of the chamber is assigned at least  one unique mesh of user­defined size. For  areas of interest, such as the substrate, the  mesh can be very fine and for areas that are  not of great interest it can be more sparse. 

An example of a chamber with the mesh  defined is found in Figure 9 where the red  part indicates the target and each gray area  is a mesh element. Once the mesh is  defined the list of particles is iterated and  each particle is assigned to the mesh that it  contains its impact point. 

  Figure 9; A sputtering chamber with a 

user­defined mesh (different gray areas) and  the target (red). 

For each mesh it is now possible to calculate a distribution of the energy and incident angle that  the incoming particle will have. The number of particles that must be simulated in order to  achieve a statistically significant distribution is heavily dependent upon the mesh resolution but  usually at least 10​^6​ particles are needed. The number of particles that end up on the target is  counted in order to obtain the redeposition factor  , i.e the fraction of the particles that aref   sputtered but ends up at the target, which is calculated as 

  f = _N

total

N_target

(3.36)   

where N_total is the number of particles simulated and N_target the number of particles that are  deposited at the target. The fraction of the sputtered particles   that are deposited on mesh g_i i   is given by 

      g_i= _(N ^N_−Nⁱ ₎

total target (3.37) 

where N_i is the number of particles that are deposited on mesh  . Let the flux   be defined asi F   the number of sputtered particles per unit area and time. Using (3.37) the flux F_i for mesh i   can be calculated as 

R F_i= g_iA^Ai

total  

where   is the number of sputtered particles per unit time,   is the area of segment   andR A_i i    is the total area of all mesh segments excluding the target.

A_total  

3.4 Potentials

The simulated scattering of a particle in the gas phase will depend on which inter­atomic  potential V(r) that is used. There are a number of difference potentials to choose from but in  this paper the following has been given some extra attention; the hard­sphere potential, the  screened Coulomb potential and the Lennard­Jones potential. The potential chosen in the  simulation model is the hard­sphere potential. Since the scattering involves two species where  at least one of them is a neutral the hard­sphere potential should be good enough 

approximation. Using one of the other potentials would greatly decrease the performance of the  simulation since it is computationally much easier to implement and evaluate the hard­sphere  potential than the other potentials.  

Hard-sphere potential

The hard­sphere potential models each particle as a hard sphere. The potential is defined as    

(r)     r  r

V = 0 >   _gas+ r_part  (r)    r  r

V = ∞ <   _gas+ r_part   

where r_gas is the radius of the gas particle and r_part is the radius of the sputtered particle [4]. 

There exists multiple different definitions of the radius of a particle such as the van der Waal  radius and the covalent radius but they are all in the order of 0.1 nm and increases with the  number of nucleons in the atom. Using the hard­sphere potential it is possible to solve the  scattering equation (3.33) analytically which simplifies the computational implementation. 

Another unique property of the hard­sphere potential is that the scattering angle will be  independent upon the energy of the particle. 

Screened Coulomb potential

The non­screened Coulomb potential is defined as    

(r)

V = ^{Z Zgas part}_r ^q2   

where Zgas is the number of protons in the gas particle, Z_part is the number of protons in the  sputtered particle and   is the elementary charge [4]. However this yields a too large interactionq   at large inter­atom distances due to neglecting the effect of electrons. To cope with this the  potential is usually screened by a screening function Φ( )^r_a  where   is the screening length [4].a   There exists a wide variety of screening functions but one of the more commonly used ones is  the ​Molière​ screening function defined as 

( ) .35e .55e .1e

Φ ^r_a = 0 ^−0.3ra+ 0 ^−1.2ra+ 0 ^−6ra   

where the screening length   is defined asa    

.885a (Z )

a = 0 _{0 gas}^0.5+ Z_part^{0.5 −32}   

with a₀ being the Bohr radius [10]. In contrast to the hard­sphere potential the scattering angle   depends on the energy of the particle when the screened Coulomb potential is used. It is not

ϕ  

possible to analytically solve the scattering equation (3.33) for this potential but there exists a  number of numerical methods that can be used [4]. 

Figure 10; The scattering angle   for a tungsten atom scattered against an argon atom with aϕ   screened Coulomb potential used as the interaction potential. 

In Figure 10 the scattering angle   as function of the particle’s energy and the impactϕ  

parameter   is plotted for a collision between a tungsten atom and an argon atom. Figure 11p   shows the same plot but for a collision between a sulfur atom and an argon atom. Note that the  maximal amount of scattering is much less for the tungsten collision and that a head­on collision  (corresponding to a very small value of the impact parameter) will backscatter the sulfur but not  the tungsten.  

Figure 11; The scattering angle   for a tungsten atom scattered against an argon atom with aϕ   screened Coulomb potential used as the interaction potential. 

Lennard-Jones potential

The Lennard­Jones potential, also known as the 6­12 potential, is defined as   

(r) ε  

V = 4

( ⁽

)

(

)

)

where   is the depth of the attractive well and   is the distance at which the repulsive andε σ   attractive parts are equal, i.e the potential is zero [11]. The constants   and   are dependentε σ   upon the material and have been determined experimentally. The potential is loosely based on  quantum mechanics perturbation theory where the positive term represents the strong repulsive  force and the negative term the weak attractive force [11]. In contrast with the hard­sphere  potential and the screened Coulomb potential the Lennard­Jones potential can be either  repulsive or attractive depending on the scattering. This manifest itself in that the scattering  angle   may be negative for certain combinations of the impact parameter and the particle’sϕ   energy. However since the scattering plane, as described in the previous section, is symmetrical  and fixed in the initial frame of reference by a azimuthal angle   the negative scattering angle isω   equivalent to a positive scattering angle of the same size in the scattering plane with azimuthal  angle ω + π. A plot of the scattering angle is found in Figure 12 where it is clearly visible that the  scattering angle may be negative. If the Lennard­Jones potential, as seen in Figure 12 is 

compared to the screened Coulomb potential, as seen in Figure 11, it can be seen that they  resemble each other but there are some differences. The scattering angle approaches   faster0  

for the Lennard­Jones potential when the impact parameter decreases. It can also be seen that  the impact parameter must be slightly larger to achieve the same amount of scattering when the  Lennard­Jones potential is used. 

Figure 12; The scattering angle   for a tungsten atom scattered against an argon atom withϕ   the Lennard­Jones potential used as interaction potential. 

3.5 Simulations and analysis

A series of simulations have been performed to investigate how the distribution of the particles  across the chamber is affected by different process parameters. The chamber used is cylindrical  in shape with a diameter of 0.32 m and a height of 0.25 m. The substrate is 11.6 cm in diameter  and is placed on a substrate holder at the bottom of the chamber. The target is specified to have  a diameter of 4 cm and is placed at the center in the top of the chamber. It consist of WS​2​ with  the radius of tungsten set to 193 pm and the radius of sulfur to 88 pm [12]. The masses used in  the simulation are 183 u for tungsten and 32 u for sulfur. The target potential is 400 V, the  sputtering gas is argon and at least 10​^6​ particles have been simulated in each simulation. 

Difference depending on sputtered element

Different elements will scatter differently due to the difference in mass as discussed in Section  3.3. The flux of sputtered tungsten and sulfur is shown In Figure 13 and Figure 14 respectively. 

The gas pressure is 4 mTorr in both simulations. It can be seen that the flux of tungsten at the  substrate is slightly higher than the flux of sulfur. At the top of the chamber however the flux of  sulfur is much larger than the flux of tungsten. This since the mass of the tungsten will make it