Study of particle transport in high pressure sputter deposition
process
Takeo Nakano
December 2001

Preface
The objective of this study is to propose the new methods to simulate the transport process of
atoms and molecules in high pressure (5–50 Pa) sputter deposition process. More speciﬁcally, we
extend the conventional Monte Carlo (MC) simulation from two standpoints. One is to introduce
the eﬀect of the thermal motion of ambient gas atoms/molecules, which has conventionallybeen ignored, into the MC calculation of collision and scattering process between the sputteredparticles and the gases. The other is to formulate the random diﬀusion process of the sputteredparticles, which occurs after they decelerated, by Poisson’s equation sand to solve it using theBoundary Element Method (BEM).
Sputtered particles have large kinetic energy when they are ejected from the target, but their
speed is decelerated by the collision with gas atoms so that the gas motion becomes graduallyimportant. This is especially signiﬁcant in high pressure environment, since the collision occurs
more frequently. The ﬁrst extension is to treat this gas motion properly based on the kinetic
theory of gases.
When the pressure (hence the collision event) increases more, the calculational complexity
proportional to the square of it is required to trace the sputtered particle by MC method until
they deposit onto the chamber wall. The second extension aims to reduce this. It treats thetransport of the sputtered particles by the diﬀusion equation, after they lose their large initialenergy (and momentum).
Simulation program has been implemented based on these ideas, and applied to the following
problems to check the validity of the model as well as to understand the physics behind these
phenomena.
• Deceleration process of high energy atoms in gas environment. It’s a typiﬁcation of thetransport of sputtered atoms from the target. Especially, the eﬀect of the thermal gas
motion on this transport behaviour is speculated by comparing the cases with/without
the gas motion.
• Thickness proﬁle of sputtered copper ﬁlms. Firstly, the eﬀect of the distance from thechamber wall on the thickness proﬁle is studied. At higher pressures, ﬁlms also deposit onchamber walls from which the target cannot be seen directly. The pressure dependence of
this diversion phenomenon is investigated, and its relationship with the particle transport
process is discussed.
• The deviation of the ﬁlm stoichiometry from that of the composite target material, andits dependence on the gas pressure. Lanthanum hexaboride (LaB6 ) is used as the modeli

PREFACE

iicase. The reason of this is discussed considering the diﬀerence in the transport process
between lanthanum and boron atoms, which reﬂects the mass diﬀerence of these.
• Evaluation of spatial density of sputtered atoms. Using the extended method introducedin this study, we can trace the motion of sputtered atoms precisely. We applied this feature
to calculate the atomic density of these sputtered atoms. When you observe the opticalemission from Cu sputter plasma, the emission including the ground state, as the ﬁnal stateof its transition, is re-absorbed by the ambient Cu atoms in ground state. Therefore, thehigher the atomic density is, the lower the intensity of these emissions are. The calculated
atomic density using our simulation model is shown to be in accordance with the one
expected from the decrement of this kind of emission.
The contents of this thesis are as follows.
Chapter 1 is the introduction. Researches on sputter deposition process are reviewed mainly
from the theoretical viewpoint, and the modeling/simulation of the sputter deposition is mentioned. Subsequently, studies treating the MC simulation of particle transport process are
summarized. By showing the problems in these previous studies, the objective of this study isclariﬁed.
In chapter 2, details of the particle transport simulation used in this study are described.
By reconsidering the physical backgrounds of them, their validity and limitation are considered.
Next, the two methods introduced in this study are explained; incorporation of the gas motion
into the collision of sputtered particles with gases, and the treatment of thermal diﬀusion processof sputtered particles via diﬀusion equation. The latter will be compared with the result of Monte
Carlo method, and its validity is discussed.
In chapter 3, the simulation program is applied to the problems listed above. We will check
how well the simulation reproduce the experimental results, and discuss by what reason thosephenomena occur.
Finally, the study is summarized in chapter 4.
The proposed model in this study is implemented into the FORTRAN/C program working on
Linux operating system. In appendix A, this program is described; availability, usage, algorisms,
structure of the program, etc.
In appendix B, atomic parameters used in this study are given in a table.

Contents
Prefacei

1 Introduction
1.1 Magnetron plasma . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
1.1.1 Magnetic conﬁguration near the target . . . . . . . . . . . . . . . . . . . .
1.1.2 Electron trajectory and ion production . . . . . . . . . . . . . . . . . . . .
1.1.3 Sheath formation in RF plasma . . . . . . . . . . . . . . . . . . . . . . . .
1.1.4 Gas rarefaction eﬀect . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
1.2 Processes at the target . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
1.2.1 Sputtering phenomena . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
1.2.2 Target roughening during sputtering . . . . . . . . . . . . . . . . . . . . .
1.3 Particle transport process . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
1.3.1 Thermalization of particles . . . . . . . . . . . . . . . . . . . . . . . . . .
1.3.2 Analytical calculation . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
1.3.3 DSMC method . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
1.3.4 Directional sputtering and Ionization sputtering . . . . . . . . . . . . . .
1.4 Processes on ﬁlm surfaces . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
1.4.1 Experimental studies . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
1.4.2 Simulation studies . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
1.5 Stoichiometric deviation in the sputtering of compounds . . . . . . . . . . . . . .
1.6 Monte Carlo simulation of particle transport processes . . . . . . . . . . . . . . .
1.6.1 Early developments of MC simulation . . . . . . . . . . . . . . . . . . . .
1.6.2 Eﬀects of the thermal motion of gases . . . . . . . . . . . . . . . . . . . .
1.6.3 Applications of MC simulation . . . . . . . . . . . . . . . . . . . . . . . .
1.6.4 Extensions of MC simulation . . . . . . . . . . . . . . . . . . . . . . . . .
1.7 The goal of this study . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

1
3
3
4
5
5
5
6
6
7
8
8
9
10
11
11
13
14
15
15
17
17
18
20

2 Theory
2.1 Concept of the Monte Carlo method . . . . . . . . . . . . . . . . . . . . . . . . .
2.2 Particle ejection from the target . . . . . . . . . . . . . . . . . . . . . . . . . . .
2.2.1 Ejection position . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
2.2.2 Energy distribution . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
2.2.3 Angular distribution . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

22
22
23
23
24
26

iii

CONTENTS
2.3

iv

Collision and scattering with gases . . . . . . . . . . . . . . . . . . . . . . . . . .
2.3.1 Mean free path . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
2.3.2 Potential scattering problem . . . . . . . . . . . . . . . . . . . . . . . . .
2.3.3 Scattering potential . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
2.3.4 Collision parameter . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
Trapping on the chamber boundary . . . . . . . . . . . . . . . . . . . . . . . . . .
MC calculation with the thermal motion of gases . . . . . . . . . . . . . . . . . .
2.5.1 Collision frequency and mean free path . . . . . . . . . . . . . . . . . . .
2.5.2 Determination of the velocity of colliding gas . . . . . . . . . . . . . . . .
2.5.3 Particle velocity after the scattering . . . . . . . . . . . . . . . . . . . . .
Particle transport after the thermalization . . . . . . . . . . . . . . . . . . . . . .
2.6.1 Diﬀusion equation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
2.6.2 Relationship with MC simulation . . . . . . . . . . . . . . . . . . . . . . .
2.6.3 Boundary integral equation . . . . . . . . . . . . . . . . . . . . . . . . . .
2.6.4 Axisymmetric case . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
2.6.5 Boundary element method . . . . . . . . . . . . . . . . . . . . . . . . . . .
2.6.6 Comparison with MC method . . . . . . . . . . . . . . . . . . . . . . . . .
2.6.7 Thermalization threshold . . . . . . . . . . . . . . . . . . . . . . . . . . .

27
27
28
30
30
32
35
35
39
40
42
43
43
45
46
47
49
51

3 Results and Discussions
3.1 Experimental Apparatus and Simulation policies . . . . . . . . . . . . . . . . . .
3.1.1 Sputtering chamber . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
3.1.2 Summary of simulation policies . . . . . . . . . . . . . . . . . . . . . . . .
3.2 Deceleration and diﬀusion of high energy atoms in gases . . . . . . . . . . . . . .
3.3 Thickness proﬁle of sputtered Cu ﬁlms . . . . . . . . . . . . . . . . . . . . . . . .
3.3.1 Experiment and simulation . . . . . . . . . . . . . . . . . . . . . . . . . .
3.3.2 Thickness proﬁle . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
3.3.3 Diversion of sputtered atoms onto the back face . . . . . . . . . . . . . . .
3.3.4 Summary . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
3.4 Composition of LaB6 sputtered ﬁlms . . . . . . . . . . . . . . . . . . . . . . . . .
3.4.1 Experiments . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
3.4.2 Simulation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
3.4.3 Mass diﬀerence eﬀect on energetic particle transport . . . . . . . . . . . .
3.4.4 Discussion and summary . . . . . . . . . . . . . . . . . . . . . . . . . . . .
3.5 Density of sputtered atoms in the plasma . . . . . . . . . . . . . . . . . . . . . .
3.5.1 Density calculation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
3.5.2 Optical emission measurements . . . . . . . . . . . . . . . . . . . . . . . .
3.5.3 Self absorption eﬀect of optical emission . . . . . . . . . . . . . . . . . . .
3.5.4 Summary . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

54
55
55
58
60
65
65
66
69
71
71
73
75
77
80
82
83
88
93
96

4 Conclusion
4.1 Future works . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

97
99

2.4
2.5

2.6

CONTENTSv

A Implementation
100
A.1 Algorisms and libraries . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 100
A.1.1 Integration . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 100
A.1.2 Random number generation . . . . . . . . . . . . . . . . . . . . . . . . . . 101
A.1.3 Spline function . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 101
A.1.4 Matrix calculation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 102
A.1.5 Special functions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 103
A.2 Simulation Program . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 103
A.2.1 Availability and requirements . . . . . . . . . . . . . . . . . . . . . . . . . 103
A.2.2 Directory structure of the distribution archive . . . . . . . . . . . . . . . . 104
A.2.3 Source ﬁle and routines . . . . . . . . . . . . . . . . . . . . . . . . . . . . 106
B Atomic Parameters

112

Acknowledgement (in Japanese)

115

Bibliography

116

List of Figures
1.1

Schematics of the sputter deposition process. . . . . . . . . . . . . . . . . . . . .

2

2.1
2.2
2.3
2.4
2.5
2.6
2.7

Potential scattering in the center-of-mass system. . . . . . . . . . . . . . . . . . .
Born Mayer potential between Cu and Ar, and its potential scattering . . . . . .
vp , vg , vR and θ. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
Dependence of the particle mean free path on its speed. . . . . . . . . . . . . . .
Dependence of vR distribution on vp . . . . . . . . . . . . . . . . . . . . . . . . . .
Comparison between the BEM simulation with MC. . . . . . . . . . . . . . . . .
Eﬀect of the thermalization threshold on the particle transport simulation. . . .

29
31
36
38
40
50
52

Sputter chamber for ﬁlm deposition experiments. . . . . . . . . . . . . . . . . . .
Radial depth proﬁle of erosion tracks formed on Cu and LaB6 targets. . . . . . .
Axisymmetrical chamber boundary used with MC simulation. . . . . . . . . . . .
Time evolution of the velocity space distribution of high energy Cu atoms injected
into Ar gas environment. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
3.5 Time evolution of the energy distribution function of high energy Cu atoms injected into Ar gas. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
3.6 Time evolution of the real space distribution of high energy Cu atoms injected
into Ar gas environment. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
3.7 Time evolution of the sputtered Cu atoms in deposition system. . . . . . . . . . .
3.8 Thickness proﬁle of sputter-deposited Cu ﬁlms. . . . . . . . . . . . . . . . . . . .
3.9 Ratio of deposition rates at front/back face of the substrate holder. . . . . . . . .
3.10 Distribution of thermalized positions of sputtered particles. . . . . . . . . . . . .
3.11 Result of ICP-Mass measurements of LaB6 sputtered ﬁlms. . . . . . . . . . . . .
3.12 Composition of sputtered LaB6 ﬁlms. . . . . . . . . . . . . . . . . . . . . . . . . .
3.13 Scattering angle of B and La atoms when they collide with an stational Ar atom.
3.14 Evolution of the velocity distribution of B and La atoms injected into Ar gas. . .
3.15 Propagation of B and La atoms in Ar gas. . . . . . . . . . . . . . . . . . . . . . .
3.16 Chamber model used in the MC simulation for atomic density evaluation. . . . .
3.17 Pressure dependence of deposition rate at the center of the substrate holder. . . .
3.18 Density proﬁle of Cu atoms. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
3.19 OES observation system. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
3.20 Optical emission spectra from Cu/Ar DC sputter plasma. . . . . . . . . . . . . .

56
57
59

3.1
3.2
3.3
3.4

vi

61
63
64
67
68
70
72
74
76
78
79
81
84
86
87
90
91

LIST OF FIGURES
3.21 Pressure dependence of optical emission intensities from Cu/Ar magnetron sputter plasma. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
3.22 Discharge power dependence of optical emission intensities from Cu/Ar magnetron sputter plasma. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

vii

92
94

A.1 Call tree of HISPUT subroutines. . . . . . . . . . . . . . . . . . . . . . . . . . . . 111

Chapter 1

Introduction
Film preparation using sputtering phenomena is classiﬁed in the Physical Vapor Deposition(PVD) process, and has been widely used in these days [105,161]. For industrial purposes, sputter deposition frequently accompanies the low pressure discharges. By applying a high-voltage
to noble gas of 0.1–10 Pa, you can generate a glow discharge plasma. Then a sputtering target
is introduced into the plasma, and negative bias is applied to it, which cause the bombardmentof the target by positive ions and the ejection of target materials via sputtering phenomenon.
These ejected particles travel through the gas environment, and some of them will deposit onto
the substrate and form the ﬁlm. The deposition system using this method can be constructed
rather easily; evacuation of the chamber, introduction of the discharge gas and the applicationof the high voltage. Though the sputter deposition had some disadvantages in the early stages,
they have been conquered through developments of new techniques: e.g., magnetron sputtering
has enabled the faster deposition, and RF sputtering has made this applicable to the insulatingtargets. And now sputtering is used with almost all the materials.
The merit of the sputtering method, compared to the ordinary vacuum deposition, is thepresence of high energy particles that hit the substrate. It has been reported that the denseness,
the adhesion, and/or the strain of the ﬁlm are improved or controlled [79, 116, 118, 164, 165] by
the bombardment with moderated energy and ﬂuence.
On the other hand, however, it is not easy to determine the optimum deposition condition in
sputter deposition. It is because that relationships are not straightforward between the experimental parameters and the internal parameters of the process (former include chamber shapes,
pressure of the gas, discharge power / voltage, etc., while the latter include electron temperature / density of the plasma, ion generation frequency, sputtering yield of the target, etc.).
They are summarized in ﬁgure 1.1. It is generally diﬃcult to predict the ﬁlm property or the
thickness proﬁle from experimental parameters, and cut-and-try method with many repetitionsof experiments is necessary to optimize such conditions.
The topic of this study is the theoretical / experimental research on the “Particle Transfer”
in the ﬁg. 1.1, where the particles ejected from the target are transported onto the substrate. In
this ﬁeld, the usefulness of the Monte Carlo (MC) method is widely accepted today. However,the conventional MC simulation suﬀer from the discrepancy between its results and experiments
at high pressure region, where much more collision of sputtered particles occur.

1

CHAPTER 1. INTRODUCTION

2

EXP. Conditions
Gas (element, pressure)
Magnetic Field
Discharge PowerParticle
Transport
Collision & Scatteringwith gas atoms

PLASMAIon Production
Electric / Magnetic
fieldsSUBSTRATE
Resputtering
MixingDefect formation
Surface Diffusion

Diffusion(post Ionization)

Chamber Wall:
Particle Flux

TARGET
Formation of theerosion track
Reflection ofincident ions

Incidence angle
Incidence Energy

Film Properties
Thickness Profile

Partielc ejection (θ ,E)

Figure 1.1: Schematics of the sputter deposition process.

CHAPTER 1. INTRODUCTION

3

In this study, we try to extend the MC framework and enable more reasonable and realistic
simulation at high pressures. We also proposed the method to reduce the calculational complexity, which increases at the square of the gas pressure in MC. By comparing the experimental
results with simulation we developed, we discuss the physical processes of sputtered particles
and their behavior in high pressure environment.
In this chapter, we review the studies on the sputter deposition. In the ﬁrst 4 sections,overall sputtering process is roughly categorized based on ﬁg. 1.1; plasma generation, target
processes, particle transport, and atomic processes at the substrate. Studies on those topics
will be summarized mainly from theoretical viewpoint, to clarify the position of the MC studyof the transport process. The next section refers some experimental results on the sputtering
deposition of compounds. It has been reported that the ﬁlm composition deviates from thetarget, which will be discussed in this study in conjunction with the particle transport process.
The next section reviews historically the Monte Carlo method applied to the transport process
in detail. Finally, the purpose of this study is presented in the last section.

1.1

Magnetron plasma

Fundamental processes occurring in low pressure discharges used for ﬁlm deposition have been
attracted many interests [69].
For sputter deposition process, development of the magnetron plasma has the signiﬁcant
importance. It has been proposed in 1960’s by Kay et al. to generate and maintain the plasma
at lower gas pressures [158,161]. It uses magnetic ﬁeld to conﬁne electrons near the target, which
enhances the ion production, enables the plasma ignition at low gas pressure, and increase greatly
the deposition rate of ﬁlms.
On the other hand, it has some disadvantages caused by its spatial inhomogeneity. For
example, non-uniform of positive ion production generates the “erosion track” on the target,
which reduces the eﬃciency of the target usage. The inhomogeneity of the plasma also makes
diﬃcult the theoretical treatment of it.
From the viewpoint of the transport process, the energy and ﬂuence proﬁles of incident ions
onto the target aﬀect the ejection of the sputtered particles. In this section, we mention several
reports on this topic.

1.1.1

Magnetic conﬁguration near the target

Since its early years, design and optimization of the magnetic conﬁguration of the magnetron
has vigorously been studied. It is because that it aﬀects the eﬃciency of the target consumption,
stabilization of the discharge, and the reduction of substrate heating via the conﬁnement of theplasma (for example, the topics in 1980’s has been summarized by Almeida [2]). In 1991, Raoet al. have studied the ring type magnetron with various combinations of magnets and yokes,
and pointed out [111]:• Stronger magnetic ﬁeld maintains the discharge at lower pressure, or generates larger
current at constant voltage.

CHAPTER 1. INTRODUCTION

4• Spatial distribution of the magnetic ﬁeld determines the uniformity of the erosion and/or
the trapping eﬃciency of electrons. If the latter becomes worse, the substrate is heated
up by those escaped from the conﬁned plasma region.
• If electrons escape, cathode current decreases in high voltages.
On the contrary, “unbalanced magnetron” conﬁguration is proposed in 1986, where the yoke
is saturated and the magnetic ﬂux escapes above the target [167]. With these ﬂux from thetarget to the substrate, ions and/or electrons are conveyed toward the target, by which the ﬁlm
growth process will be controlled. Howson et al. has performed the probe measurement with thisunbalanced type magnetron plasma, and obtained the ion energy distribution near the target
as a function of gas pressure and discharge power [50]. By insulating the substrate, 30 V ofsubstrate self-bias and 100 mW/cm2 of heat supply to the substrate were observed.
In 1995, the group of Musil and Kadlec have developed the experimental system with twoelectromagnets in between / outside of two cylindrical yokes. The system can generate bothbalanced and unbalanced magnetron conﬁguration by changing the current ratio of these electromagnets [96]. With this system, they observed that the dependency of discharge ignitionpressure on the current ratio, and found its minimum at Iouter /Iinner ∼ 2. By calculating themagnetic ﬁelds, they showed that the area of the target face covered by the magnetic ﬂux became
maximum at this current ratio [58].

1.1.2Electron trajectory and ion production

As mentioned above, the ion production proﬁle is inhomogeneous in the magnetron plasma. It
is important to predict this proﬁle for the simulation of sputter deposition process, since it gives
the way to calculate the erosion formation from the magnetic conﬁguration.
Sheridan et al. [132] have developed the method to do this using Monte Carlo simulation of
the electron motion. They traced the electrons produced at the target by secondary electron
emission process. They solved the equation of motion of the electron in the applied electric andmagnetic ﬁelds. During the motion, electron collides with gas atoms, and ionizes them withsome probability. Those ions are then dragged onto the target, bombard it and produce anothersecondary electrons. Iterating this processes, they could obtain the arrival distribution of ions(=depth proﬁle of the erosion) for given magnetic conﬁguration, which agreed fairly well withexperiments.
Ido et al. at Saitama University have extended this MC simulation method, and appliedthem to the erosion track formation on the target of magnetic materials [52–54]. Shidoji et al.
also have used this method for the rectangular target, and obtained the good agreement withthe experimentally obtained erosion depth proﬁles.
Though the method of Sheridan has been used widely, it has to assume the electric ﬁeld
proﬁle before the simulation, hence is not an a priori calculation. But there have been several
proposals recently, to predict the sheath potential of plasma without magnetic ﬁelds [10] as
well as for Ar magnetron plasma [134]. Both of these are the combinations of ﬂuid dynamical
approach and of kinematical one for the electron motions. If they are merged with the Sheridan’s
method, it will be possible to predict the erosion track from the magnet conﬁgurations, and to
optimize it during the design. The probe measurement of the plasma parameters and the spatial

CHAPTER 1. INTRODUCTION

5

distribution of them has also been performed by Sheridan’s group themselves [133], which will
be compared with theories and activate the development of the theoretical considerations in the
future.

1.1.3

Sheath formation in RF plasma

In DC discharges, the potential drop at the cathode sheath is almost equivalent to the discharge
voltage, and the energy of ions bombarding the target is almost identical to sheath potential
drop times the ion charge1 . On the other hand, the ion energy distribution in radio frequency(RF) discharge is more complicated and has been extensively studied.
The study reported by Tsui [150] in 1968 was one of the pioneering work in this ﬁeld. He
solved equations of motion for ions and electrons under the electric ﬁeld with both RF and DC(self-bias) components. The resultant ion energy distribution function (IEDF) has two peaks at
top and bottom energies of the distribution.
This shape of distribution has conﬁrmed by the experiment by Coburn and Kay [16], who
measured the IEDF at the DC biased substrate exposed to the RF plasma. Recent experiment [77] and calculation [76] also supports the IEDF of Tsui’s.

1.1.4

Gas rarefaction eﬀect

For the particle transport, gas density is a critical parameter since it determines the collision
frequency of sputtered particles. It has been proposed that the gas density is reduced near
the discharge target, which is called “gas rarefaction eﬀect”. In sputter deposition, sputteringtarget can be regarded as the source of various kind of high energy particles, such as reﬂected
and neutralized ions, charge exchanged atoms, and sputtered target atoms. They heat the gas
near the target region which reduces the density of them (assuming the pressure is constant).
This eﬀect has been proposed by Rossnagel [115], based on the concept of “sputtering wind”
by Hoﬀman [45]. Rossnagel has measured this local reduction directly by inserting the tube
near the cathode, and proposed the model that the gas density is proportional to the one-second
power of the discharge current.
He also has studied the magnetron sputtering discharges by various techniques, and discussed
the ramiﬁcation of this eﬀect. It includes the saturation of the optical emission intensity at highpower region [122], and the energy distribution of ions incident on the substrate [116].

1.2

Processes at the target

Sputtering process on solid state surfaces has initially been studied mainly motivated by thesurface damages introduced by radiations. And as the sputter deposition technique came to be
used widely, energy and angular distributions of sputtered particles, which are necessary for thesimulation, have also been attracting interests.

1

if we can neglect the charge exchange process [21].

CHAPTER 1. INTRODUCTION

1.2.1

6

Sputtering phenomena

In the ﬁeld of theoretical study of sputtering, Thompson’s work for energy distribution [146]and Sigmund’s for the sputtering yield [135] have been well known and cited frequently. Since
the former is especially important for the MC simulation, it will be mentioned in detail in 2.2.2.
After their work, Monte Carlo codes (MALOWE, TRIM, ACAT, etc.) have been developed,
and been applied to obtain various properties of sputtering: angular distribution [61,103], energy
distribution [15,176], dependence of sputtering yield on primary ion energy [128], and the depth
proﬁle of ejected particles [107].
Experimental studies have also been carried out intensively (most of which has been compiled
by Mahan [72]). Especially, the angular proﬁle of sputtered particles measured by Tsuge forpolycrystalline metal [149] has actually been used by several groups in the particle transport
simulation. In view of experimental technique, the work by Whitaker et al. is also interesting, inwhich they applied the Laser Induced Fluorescence (LIF) to measure the angular distributionsof sputtered particles accurately.
Though it’s rather diﬃcult to use in the study of sputter deposition process, it should be
noted that there is a comprehensive data set of sputtering yields for elemental targets at normal
incidence of primary ions [177]. It is the compilation of experimental data, part of which arecompensated by the simulation using ACAT . Mahan has used this kind of data, and discussed
the periodicity of sputtering yield with the atomic weight, as well as the eﬀect of nonlinear
cascades [73].
Recently, several calculations have been presented for the selective sputtering phenomenaon the alloy / compound targets, including the eﬀects of atomic diﬀusion [7, 55]. They will be
important as the application of those materials becomes major (see also 1.5).

1.2.2

Target roughening during sputtering

The contents of this subsection is rather oﬀ-topic from the MC simulation, the main subject of
this study, but it has been attracting interests relating to the formation of self aﬃne fractals [6].
Since surface roughening by the ion beam sputtering aﬀects the resolution of the depth-proﬁlemeasurement in surface analysis, it is also interested by researchers of surface physics and/or
chemistry [85]. Surface roughness may also aﬀect the ejection of sputtered atoms, so we will
review some reports on this topic here.
Sigmund has also made signiﬁcant contributions in this ﬁeld. According to his theory on
the sputtering yield mentioned above [135], the primary high-energy particle causes collisioncascade at some depth under the surface, and the energy of it is spread in ellipsoidal form. If
the energy delivered to the surface by this collision cascade is larger than the binding energy ofatoms there, the sputtering ejection occurs. The probability of this ejection becomes larger at
positions where the curvature is larger. Hence the random roughness grows up, and more roughsurface forms [136].
Using his theory, Bradley and Harper have successfully explained the wavy structure observedon the amorphous materials sputtered by ions with oblique incidence angle in 1988 [13]. Theirtheory has developed to explain experimental results:• At suitable incidence angle, ﬂux and target temperature, wavy ripple structure appears

CHAPTER 1. INTRODUCTION

7

on the surface.
• The 2D wave vector of the ripple structure is parallel to the surface component of the incidence vector of ions, when the incidence angle is small (i.e. the nearly normal incidence).
• If the incidence angle becomes larger than some threshold, the wave vector of the structure
rotates in 90 degrees, and becomes vertical to the surface component of the ion’s incidence
vector.
The smaller the wave vector is, the faster the roughness growth is, because of the larger curvature.
But the smoothing eﬀect by the surface diﬀusion is also faster for smaller wave vector structure.
As a result of the competition of these two eﬀects, some characteristic wave vector is determined
and the structure with it appears on the surface. Applying this model and comparing the
growth rate of structures with the wave vector parallel / vertical to the surface component of
the incidence vector, it is shown that their relationship inverses at certain incidence angle which
agrees with the experimental results above2 . Mayer et al. have carried out the experiment whereions of H, He, Xe, etc. were irradiated onto the SiO2 surface [80], and showed that the surface
roughening and smoothing were well described by the Bradley’s model.
Barabási and coworkers have developed the model which includes both the Sigmund’s energy
transfer model and the curvature dependence of the surface diﬀusion. With the model, theysummarized the ion-induced surface structure into a phase diagram with two parameters: thepenetration depth a of the primary ion and the incidence angle θ of the ions [17, 74].
Makeev and Barabási have applied this model to the sputtering of self aﬃne surfaces, andcalculated the dependence of the sputtering yield (for secondary ions) on the fractal parametersof the surface. They showed that the yield Y depended on the ion penetration depth a, surface
saturation width w and the correlation length ξ:• Y ∝ w2 when a >> w.
• Y ∝ w−2 when a << w.
• Y is larger for larger w when ξ >> w.
• Y is smaller for larger w when ξ << w.

1.3

Particle transport process

As sputter deposition has become popular, problems of it have also been recognized, e.g. theinhomogeneity of the deposition rate on substrate positions, and the generation of convex curvature near mask edges. In addition, it has also been understood that the sputtered particles
have much larger energies compared to the particles in vacuum evaporation, and aﬀects theproperties of deposited ﬁlms in various way (see 1.4). Therefore, it has been interested how
sputtered particles are transported to the substrate, and what energy, ﬂux and incidence angle
2

There has been another report using the simulation with discrete model, protesting that the ripple structureappears only in the initial stage of the surface growth, and changes into the self-aﬃne fractals [18]

CHAPTER 1. INTRODUCTION

8

they reach there with. The Monte Carlo simulation, the main topic of this thesis, is one of thesestudies.
In this section, we introduce the studies on this particle transport process, except for thoseusing MC method. MC studies will be reviewed in separated section (1.6), and other miscellaneous subjects are covered here: thermalization of particles in the ambient gas, analytical
study of the transport process, Direct Simulation of Monte Carlo (DSMC) method for particle
transport and experimental methods developed to control them during sputter deposition.

1.3.1

Thermalization of particles

Since the eﬀect of high energy particles came to be known, the behaviour of them have attracted interests. High energy particles generated on the target are decelerated by collisions
with sputtering gases. When the energy of sputtered particle becomes comparable to that ofgases, the particle is called “thermalized”. This process —thermalization— gave one of the mostsigniﬁcant motivations to the particle transport study.
In 1981, Sculler et al. have calculated how the energy of sputtered particles, ejected with
Thompson’s EDF, were quenched as they went away from the target for the cases of Nb and Cu
sputtering by Ar [83]. After that, they have applied this result to evaluate the degradation of thesteepness of Cu/Nb multilayer interfaces induced by the milling of those high energy particles
delivered to ﬁlm surface [129].
Inspired by these theoretical considerations, experiments have been performed to measure
directly the velocity distribution of particles in the plasma. Ball and coworkers have performed
this for Cu/Ar sputtering plasma by observing the Doppler’s shift of the atomic emission linefrom Cu in 1986 [5], and Park et al. have done for Fe using laser induced ﬂuorescence in 1991 [104].
The mean energy and the proportion of the thermalized particles were given as functions of the
product of the distance and the gas pressure. These studies gave an strong inﬂuence on thedevelopment of the MC simulations.
For the origin of high energy particles and the evaluation of their energies, comprehensive
review has been submitted by Winters, Eckstein, and others [168]. When impinging ions areneutralized and recoil on the target surface, they have very high energy so that they cause majoreﬀects on the ﬁlm. Window stated that heavy gas (Xe, Kr) should be used to reduce theseeﬀects, and proposed the static system without the ﬂow of expensive gases in the production
system [166].

1.3.2

Analytical calculation

To evaluate the thickness proﬁle in sputter deposition, analytical calculations were developed
at ﬁrst ( [161] pp.115). In these, deposition rate was obtained by integrating the ﬂux from eachpart of the source of sputtered particles. The intensity distribution of the source was taken to
be proportional to the depth proﬁle of the erosion track. Then the ﬂux was given by the productof the source intensity and the solid angle of the source from the substrate, the latter of which
was calculated from the distance and the angle between them. Introduction of the computer
has made this calculation applicable to any form of the chamber and the erosion track [142].

CHAPTER 1. INTRODUCTION

9The scattering of sputtered particles by gas has come to be considered since late 1970’s. For
example, Wu et al. have elucidated the deceleration of high energy particles from their mean
free path (MFP), and evaluated the damage brought to the ﬁlm surface during the deposition of
Nb3 Ge superconducting ﬁlms [171]. Westwood has noticed that the experimental thickness proﬁle and/or the form of the curvature at the mask edge show the discrepancy from the calculation,and has proposed the model where sputtered particles move diﬀusively after they travel certain
distance from the target [163]. Keller et al. have indicated that the transfer to the diﬀusive stage
was dominated by the product P d of gas pressure P and the target–substrate distance d, andlead the expressions which determined the high-energy ﬂux and the diﬀusive ﬂux [63]. These
have been conﬁrmed by Drüsedau et al. in 1998 [29]. They used several sputtering chambers
and measured pressure dependence of deposition rate for Al/Ar and Si/Ar, which were welldescribed by Keller’s expression. They also measured deposition rates for combinations of 4gases and 7 targets, and stated that the prefactor of P d, which determined the ﬂux decrease,
was aﬀected by the atomic mass ratio between the gas and the target elements.
The MC simulation followed these studies by Wu, Westwood and Keller, and becomes dominant in this ﬁeld. However, analytical methods have still been used by several groups. For
example, diﬀusive process was treated by solving the diﬀusion equation with cylinder boundary [109], and the composition of the high-Tc superconducting ﬁlm was evaluated with diﬀerent
geometry of target and substrate, considering the angular distribution of the sputtered particles [67].

1.3.3

DSMC method

Direct Simulation of Monte Carlo (DSMC) has been developed by Bird and others to treat therareﬁed gas ﬂow [8]. In this method, the considering volume is split into cells as small as a halfof the MFP, where some thousands of representative particles are manipulated assuming theBoltzmann’s distribution. And then, temperature, density and particle ﬂux are obtained in eachcell. This method has been applied to the transport process of sputtered particles since early
1990’s, and several reports were submitted.
In 1993, Urbassek et al. have applied DSMC to Cu/Ar sputtering, and obtained proﬁles oftemperature, density, pressure and velocity for both Ar and Cu [156]. But their results showed
the higher gas density near the target, which did not agree with the experimentally observed“gas rarefaction” by Rossnagel [115, 116, 122], so the validity of their calculation was rather
questionable.
Kersch et al. have reported the simulation study on Ti/Ar in 1994, using DSMC andTEPA [64]. The latter method was the derivative of the MC simulation, assuming certain
temperature / density proﬁles of Ar. DSMC and TEPA did not agree well, and it was also
diﬃcult to determine which was valid, since experimental result was not shown. In their TEPA
calculation, the VHS (Variable Hard Sphere) and M1 interatomic potentials were used, which
was rather scarce in MC simulations.
In 1999, Dosa and coworkers have applied DSMC to the simulation of Ar gas ﬂow in the large
size sputtering chamber [28]. They only treated the ﬂow of supplied Ar, and not the sputteredparticles. But the Ar gas density proﬁle is an important parameter for MC simulation, so it isgoing to be taken into account in the future.

CHAPTER 1. INTRODUCTION

1.3.4

10

Directional sputtering and Ionization sputtering

Since 1990’s, as characteristic dimensions of semiconductor devices shrink, metalization of trench
structure of high aspect-ratio has been required with void-free dense ﬁlms. It is diﬃcult for normal sputtering technique since the incidence angle distribution of arrival particles is so scatteredthat overhangs are generated at the entrance of trenches and blocks the further deposition at
the bottom of them. To solve this, so-called “directional sputtering” and “ionization sputtering”
techniques have been developed mainly by Rossnagel and coworkers at IBM [114, 117].
Collimated sputtering is one of the directional sputtering, proposed by Rossnagel in 1991 [121].
In this method, the array of tubes (=collimator) is located between the target and substrate. In
his case, each tube was 5 mm in diameter and 20 mm in length. Though it was demonstrated
that trenches with aspect ratio as large as 1:3 were ﬁlled well, it has not widely been usedbecause the droplets grow on the collimator in relatively short deposition time [117].
Long throw sputtering is an another kind of directional sputtering, where the substrate is
located far away from the target, to converge the incidence angle of particles to surface normal
and improve the trench ﬁlling characteristics. It is known that the reduction of gas pressure
decrease the scattering of particles, and more sputtered particles directly reach the substrate [19].
Moreover, in long throw sputtering, the target–substrate (T–S) distance is increased and theangle subtended by the ring-like evaporation source (erosion track on the target) at the substrate
is reduced. By generating the discharge at as low as 10−1 ∼ 10−3 Pa and increase the MFP,
it is insisted that the sputtered particles arrive at the substrate preserving its initial direction
at the target3 . Motegi et al. at ULVAC Corp. have deposited Ti, TiN and Al alloy under thepressure of 0.04 Pa and T–S distance of 31 cm, and obtained the bottom coverage4 of 50 % for
the hole with 1:2 ratio [89]. Mayo et al. have studied this technique both experimentally andtheoretically [81]. Their conditions were T–S distance of 30 cm, gas pressure as low as 0.1 Pa,
and target / substrate diameters of 30 cm and 20 cm, respectively. They compared the side /bottom coverage of experiment and calculation, and indicated that long throw sputtering had
the advantage on side wall coverage.
Ionized sputtering [47] has proposed also by Rossnagel and coworkers in 1993. In this method,
RF coil is inserted between the target and substrate to enhance the ionization of sputteredparticles there. Then the direction of these ions are controlled by applying the bias to thesubstrate. In their ﬁrst paper, they located the thickness monitor behind the grid electrode,and evaluated the ionization percentage by comparing the deposition rate with or without the
repelling bias. When magnetron power was 1 ∼ 3 kW and the coil power was 100 ∼ 1000 W,
50 ∼ 90 % of ionization eﬃciency was observed under 4.7 Pa [119]. In this condition, thesubstrate was inserted and the bias voltage (20 ∼ 40 V) was applied to it. For contact holes withaspect ratio of 5 ∼ 6 on the substrate, bottom coverage of 50 % was obtained [120]. They also
simulated the deposition process in this method, and reproduced well the deposition proﬁle inside
the trench, considering the resputtering eﬀect, sticking probabilities, and angular distributionsof incident neutrals / ions / ambient gas ions [40]. There have been some attempts to extend
the ionization sputtering: Magnets were arranged to form a multipole magnetic conﬁguration
3

But the MFP is still some centimeters in length, and much shorter than the T–S distance in these pressures.
”bottom coverage” means the ratio of the deposition rate at the bottom of the trench or of the hole, comparedto the substrate surface outside of them. “side(wall) coverage” is also deﬁned as well.
4CHAPTER 1. INTRODUCTION

11around the RF coil which improved the ionization eﬃciency [159], or some hundreds gauss ofmagnetic ﬁeld was applied to penetrate the RF coil to form helicon wave5 in the plasma [160].
The latter is called “helicon plasma sputtering” by several groups, and has been applied to the
production of optical multilayers of Si 2 /TiO2 , which show the well deﬁned interfaces [40].

1.4Processes on ﬁlm surfaces

In sputter deposition process, high energy particles arrive at the ﬁlm surface and cause various
phenomena there, as already mentioned above several times. These eﬀects were noticed by many
researchers in early 1980’s, when the sputter deposition had become popular by the developmentof magnetron sputtering apparatus. Since then, many experimental studies have been reported
about the ﬁlm properties deposited by sputtering. These knowledge was summarized and sharedat the end of 1980’s.
In 1990’s, as the performance of the computer improved and the simulation methods (suchas molecular dynamics) were developed, simulation of the growth process of sputtered ﬁlms has
been studied by several groups. These studies were enhanced by the requirement from VLSI
industries, in which trenches and holes with high aspect ratio must be ﬁlled and metallized by
the sputter deposition.
In this study, our goal is not to simulate the detailed properties of deposited ﬁlms, but toobtain frequency, energy and angular distributions of arrival particles. However, we believe that
it is meaningful to review here the growth process of sputtered ﬁlms and eﬀects of high energy
particles on it.

1.4.1

Experimental studies

As described above, it was known in the beginning of 1980’s that high energy particles broughtsigniﬁcant eﬀects to ﬁlm properties. Moreover, it was also noticed that applying the bias voltage
to the substrate changed the ﬁlm characters: incorporation of discharge gases (e.g. Ar), ﬁlm
stress, density, electrical conductivity, and so on.
Ziemann and Kay proposed a measure “normalized energy” to explain these phenomena [179].
It is the energy delivered onto the ﬁlm surface per depositing atom. To study the eﬀect of thisparameter, they prepared the Pd ﬁlm with dual ion beam sputtering system. It was shown that
normalized energy well described the variation of ﬁlm properties: the gas (Kr) incorporation,
lattice distortion, grain size, orientation and electrical conductivity. Just after their study,
Hoﬀman have reported the dependence of ﬁlm properties on gas pressure, discharge current andsubstrate position. It was concluded that their behaviors were described qualitatively by theconcept of normalized energy.
In 1985, Kay and coworkers have also used the dual ion beam apparatus for Ag, and showedthat the change of ﬁlm orientation was dominated by the competition between strain energy
and the surface energy: the former was originated from the biaxial stress generated in the ﬁlm
volume, and the latter was dependent on the lattice face appeared on the ﬁlm surface [51]. They

5

About helicon wave in the plasma, see [108] for example.

CHAPTER 1. INTRODUCTION

12have also conﬁrmed this phenomenon for sputtered Ni ﬁlms deposited by conventional sputter
deposition using Ar discharge plasma [62].
For the structure of sputter deposited ﬁlms and their dependency on the process parameters,
Thornton’s zone model [148] has generally been accepted since 1980’s. It takes the gas pressure
and the substrate temperature6 as parameters, and categorize the region into 4 zones (Zones
1, T, 2 and 3), each of which have characteristic microstructures. Thornton has discussed therelation between this model and the eﬀects of high energy particles, and explained why thesezone structures were grown in respective conditions [147].
In 1989, reviews of these experimental results have been compiled independently by Mattox [79] and by Rossnagel [118]. Both have summarized the phenomena by normalized energy,
and become the standard way to understand the growth of sputtered ﬁlms. They mainly treatedthe sputter deposition using plasma discharges, where the high energy particles were originatedby:• ions accelerated by the plasma sheath near the substrate, which are aﬀected by the substrate bias (mostly charged positive)• atoms / molecules sputtered from the target (electrically neutral)• primary ions bombarding the target which are neutralized by the charge exchange process
there, and recoil (neutral)• negative ions generated near the target, and inversely accelerated by the cathode’s potential drop (negative charge)
These high energy particles aﬀect the deposited ﬁlms in various way including:• desorption of adsorbates on ﬁlm surface• secondary electron emission, which may aﬀect the chemical process at the surface• resputtering of surface atoms and redeposition of them• enhancement of surface diﬀusivity of adsorbed atoms• increase in nucleation density by the introduction of surface defects• implantation of high energy particles into ﬁlms
• recoil implantation of surface atoms into ﬁlms
• lattice defects caused by the collision cascade generated underneath the surface• increase surface temperature

6

more accurately, the ratio of substrate temperature and the melting temperature of the ﬁlm.

CHAPTER 1. INTRODUCTION

13

These primary events then caused changes in grain size, defect density and internal stress, whichwere summarized as a function of normalized energy. We cannot go into the detailed discussion
here, but can roughly summarize that “moderate normalized energy brings good eﬀects to theﬁlms, but too high energy aﬀects rather badly”.
It should be worth to mention here, that the stress in sputtered ﬁlms has thoroughly beenreviewed by Windischmann in 1991 [164] and 1992 [165], including the mechanisms of the stressgeneration for both tensile and compressive cases, and the compilation of vast of experimental
results. Hoﬀman has also discussed the stress transition between tensile and compressive [46],relating to elements of gases and targets, deposition rate, substrate angle, etc.

1.4.2

Simulation studies

As the general understanding for the eﬀect of high energy particles has been common, simulation
studies of the ﬁlm growth process with high energy particles appeared in late 1980’s.
Müller have applied the molecular dynamics (MD) method to the simulation of sputteredﬁlms in 1987 [95]. He arranged atoms randomly on the substrate, and supplied additional
atoms from the above of it. Lennard-Jones interatomic potential was used for the pair ofatoms composing the ﬁlm, and Moliére potential was used between the bombarding high energy
particle and the ﬁlm atom. Atoms just above the substrate were ﬁxed, which meant the substratetemperature of zero. The result of the simulation showed tensile stresses in any condition. The
strength of the stress increased as energies of depositing atoms and/or the high-energy particles
increased, took a maximum at a certain energy, and decreased after that. When the energy was
nearly zero, ﬁlm structure became columnar because of the shadowing eﬀect. Near the peak of
the stress, these columns coalesced and caused the tensile stress by the mixing of atoms.
In 1993, Fang et al. have reported the MD simulation also [34]. It was based on the Müller’s,
but more realistic modeling was developed and applied. Their model could account for the gas
incorporation into the ﬁlm, so that the generation of the compressive stress could be reproduced [33]. In the second report, they used the MC simulation to obtain energy and angular
distributions of incoming particles.
Apart from these MD simulations, a group of Smy, Tait and Brett has developed the simulation code called SIMBAD . It simulated more macroscopic ﬁlm growth. Initial model was rather
simple, where disks were incident at certain angles, rolled on the surface, found energetically
preferable place and located there. Despite of its simplicity, the ﬁlm structure near the substrate
step and the density proﬁles were well reproduced [23, 138, 143, 144]. This SIMBAD was combined with the MC particle transport program, SIMSPUD, and was modiﬁed to use the energy
and angular distributions calculated with it [22, 24]. Independently, their group has developed
the grain growth simulation program called GROFILMS applicable to sputtered polycrystalline
ﬁlms [35]. This method handled the surface / interface diﬀusion of atoms on / between thegrain(s), as well as the ﬂow of atoms by bulk diﬀusion between grains, through the motion of
the linear nodes representing each crystalline. It has been applied to the growth of sputteredCu on the trench structure of Ta/SiO2 , and successfully reproduce the cross sectional structureof experiments [36, 37].
Through these studies, the importance of the energy / angular distributions of incident
particles have been recognized to simulate the sputtered ﬁlm growth. Therefore, some attempts

CHAPTER 1. INTRODUCTION

14

have been submitted to combine the particle transport simulation with that of the ﬁlm growth
process. They include the reports by Smy and Fang mentioned above, a study by Hoshi andcoworkers [49] which applied the Henderson’s rigid sphere model [44], and that by Ohta etal. [102] which used the Leon’s model [68] where the surface motion was accounted as theviscous ﬂow to minimize the surface energy.

1.5

Stoichiometric deviation in the sputtering of compounds

Sputtering deposition has also been widely used for alloys and compounds including magnetic /dielectric / superconducting materials [161]. It is critical to control the ﬁlm composition for theproperties of these materials, but it frequently deviates from that of the target. In this subsection, the studies on this compositional deviation are summarized. Reactive sputtering method
has also been used generally to prepare compound ﬁlms, where some of elements consisting theﬁlm are supplied as a reactive gas. But it is out of the scope of this study7 .
Since the discovery of High-Tc superconducting materials in the midst of 1980’s, many groups
have reported the preparation of the ﬁlm of it by sputtering. The compositional deviation
problem has also been recognized since the early stages of them, and the “oﬀ-axis sputtering”
method has been proposed [32, 145], where the substrate did not face to the target but was
inclined or was even located vertically to it. It was found that this method actually decreases
the extent of the deviation. The reason of the deviation were thought to be electrons and/or
negative ions generated near the target, which were accelerated by the cathode potential, hit
the substrate and caused selective resputtering or desorption of particular element [174]. The
dependence of the deviation on the pressure (P ) [97] and target–substrate distance (d) [127]showed that the deviation decreases as the P and/or d increases. It was explained that thehigh-energy particles were blocked by the scattering with ambient gas atoms [38].
Apart from these discussions that ascribed the deviation to substrate reactions, there havealso been reports that consider the diﬀerence in scattering process by gases due to the massdiﬀerence between elements [172, 173]. It has been considered in the study by Rossnagel etal. [123], where they study the pressure dependence of the ﬁlm composition of metallic alloys
e.g. AlCu and WTi. They deposited ﬁlms between 5–30 mTorr, and found the minimum of
the concentration of lighter elements at the intermediate pressure. They explained this by thetransition of the transport process from ballistic to diﬀusive. In the both low and high pressure
extremes, the diﬀerence in the element does not appear, while in the midst of them, lighter
atoms are backscattered by ambient gases, and deplete more in the ﬁlm. We also have observed
the depletion of lighter elements with pressure increase at LaB6 deposition [59], and discussed
the eﬀect of the scattering of lighter boron [65].
Moreover, selective sputtering at the target [178] has also been discussed. However, it cannot make the deviation by itself if the enough time of pre-sputtering is performed before thedeposition, and if the diﬀusion of elements in the target is not so signiﬁcant [169].

7

For more description of this method, see section IV in [105], for example.

CHAPTER 1. INTRODUCTION

1.6

15

Monte Carlo simulation of particle transport processes

As described in preceding sections, particle transport process during the sputter deposition issigniﬁcantly important for deposited ﬁlms. It dominates the thickness proﬁle of the ﬁlm, andalso aﬀects ﬁlm properties through the energy / angular distributions of incident particles.
However, it is not so easy to model the transport process and to reproduce the experiment.
It comes from two reasons. Firstly, each sputtered particle collides with gases, be scattered bythem, and changes its energy and direction until it reaches the chamber wall. Therefore, simple
geometrical calculation is not valid as in the case of vacuum evaporation. Secondly, sputteredparticles ejected from the target have energies of some electron volts, which are much higher
than the thermal energy of ambient gases [146]. Hence the transport of particles has the strong
directivity initially, and cannot be treated as the simple homogeneous diﬀusion.
To simulate this problem, Monte Carlo (MC) method has been proposed in the midst of
1980’s. It assumes probability functions for various processes that particles experience during
the transport, which construct the Markov chain. It traces 105 ∼ 106 of trial particles, andstatistically obtains the distribution of arrival positions on the chamber wall and the energy /angular distributions of particles at each positions.
In this section, we review historically the development of this method. The model generally
accepted at the time of writing will be discussed in sections from 2.2 to 2.4 of the next chapter.
Applications and the latest developments of the MC method are also introduced.

1.6.1

Early developments of MC simulation

The reports by Somekh in 1984 [139] and by Motohiro et al. [91–93] in 1983–1984 may be theearliest ones which apply the MC method to the particle transport. In this subsection, we
describe their studies in some detail.
Somekh has been motivated by the damage of ﬁlms brought by high energy particles originated from the target, and simulated the deceleration process of these particles in the gas
atmosphere. His assumptions were:• Energy of sputtered particles was ruled by the Thompson’s formula [146] and that of(neutralized and) recoiled ions obeys the Kaminsky’s formula [60]. Angular distribution of
the sputtered particles was assumed to be over-cosine distribution measured by Patterson
et al. [106].
• In the collision process with gases, diﬀerential scattering cross section proposed by Robinson [113] was used to determine the scattering angle from the collision parameter. Total
cross section σ was obtained by linearly approximating the Robinson’s data. The mean freepath kTg /P σ was deduced from it, and respective free paths were yielded using Poisson’s
distribution.
• Gas atom was assumed to be stationary on collision.
• Only elastic collision was treated. The direction and the energy of the particle were
calculated from the scattering angle.

CHAPTER 1. INTRODUCTION

16

In his model, cross section of the collision was dependent on the particle energy. As the energy
increases, the cross section became smaller. As a result, the deceleration distance became longer
than the case where the cross section of the thermal energy region was assumed. It was also
pointed out that high energy particles aﬀecting ﬁlm properties were mostly recoiled ions when
the P d (gas pressure P times the target–substrate distance d) was as large as 10 ∼ 500 Pa mm.
The report by Motohiro et al. have treated sputtered particles only. In their initial studies [91–93], the Thompson’s formula was used for energy distribution, and cosine rule was used
for angular distribution. Ambient gases were assumed to be immobile as Somekh’s model, andrigid sphere model was applied for the collision between the sputtered particle and the gas atom.
Each radius was determined by interpolating the values obtained from the viscosity of their gas.
Using this model, they deduced the probability of the sputtered particle from the target toarrive the substrate facing to it (νA ) as well as to be back and redeposit onto the target (νR )as functions of P d. In the case of Ag and P d is more than 0.1 Torr · cm (1.34 × 102 Pa mm),it was indicated that νA became almost zero, νR saturated at 0.5, and the rest of particles
were thermalized in the atmosphere. The deceleration process of the particles and the angular
distribution of them at the substrate were also discussed. Their assumptions in these reports
were rather simple compared to Somekh’s.
Motohiro’s paper in 1986 [90] was greatly improved, in which he compared the simulation
with experimental results of thickness proﬁles of Ti sputtered by Ar. The disk form substrateholder was dangled above the target, and pressure dependence of thickness proﬁles of the holder,
front and back sides from the target, were measured. In MC model, both the rigid sphere modeland the potential scattering model were used. In the latter, Born-Mayer interatomic potential
was adopted, and Abrahamson’s potential parameters [1] were used. This potential decreases
exponentially as the distance increases, hence the “cross section reduction” eﬀect proposed bySomekh was incorporated implicitly. This potential was chosen since it was more similar to
the result of Gaussian 80 compared to the Robinson’s pair potential used by Somekh. The
Thompson’s distribution was used for the energy of sputtered particles. Three types of angular
distributions; cos θ, cos3 θ and sin 2θ were tested as the angular distribution. Gases were also
assumed to be stationary. In high pressures, the sputtered particle collides with gases manytimes and their energy may become lower than that of gases, where the assumption is not valid.
In this case, he assumed that the motion of the sputtered particle became the random walk.
The mean free path of this random walk was determined from the thermal radius of the gases,
considering the factor of (1 + m/M )−0.5 which came from the gas motions. Mainly because
the computer performance, this random walk was limited as many as 50000 collisions, and theparticle which could not reach the chamber wall until then were thought to be lost. As a resultof the simulation, they concluded:• The result using sin 2θ angular distribution showed the best agreement with thickness
proﬁles of the experiment.
• Signiﬁcant diﬀerence did not observed between the result with rigid sphere model and thatwith potential scattering model.
• The discrepancy between the simulation and the experiment became considerable at pressures more than 5 × 10−2 Torr.

CHAPTER 1. INTRODUCTION

1.6.2

17

Eﬀects of the thermal motion of gases

As appeared in the report by Motohiro in 1986, it had been recognized that the decelerated
particles with many collisions, i.e., thermalized particles, was dominant in high pressures, where
the stationary gas assumption was not valid. For example, the experimentally observed energy
distribution of sputtered particles by Ball et al.8 did not agree well with the theoretical one by
Westwood [163], which used the resemble method with Somekh’s simulation.
Motivated by this, Turner and coworkers proposed the method to incorporate the thermal
gas motion into the MC simulation [153]. Taking the thermal velocity distribution of gases into
account, mean free path of the sputtered particles becomes shorter (e.g., even the particle withzero speed may be collided). Note that this eﬀect overlaps with the reduction of the collisioncross section at high energies which extend the MFP. The MFP of the particle in the gas with
Maxwellian has been given by Jeans [56] as a function of the particle speed, which was used intheir study. On collision, the gas speed was determined randomly assuming the gaussian with
FWHM of 2 kTg /Mg for each component of x, y and z. Other processes were modeled almostin the same way with Somekh’s and/or Motohiro’s; Thompson’s formula and cosine distribution
were used for the ejection of particles from the target. Lennard-Jones potential was chosen
for the scattering potential. Cu-Ar, Cu-Xe and Cu-Kr cases were studied. They calculated
the distribution of velocity components both parallel and vertical to the target surface, anddiscussed their dependence on the gas pressure and the distance from the target. The mean
speed observed from the side of the chamber were also calculated, which agreed fairly well withthe experimental results by Ball et al. [5].
With this paper by Turner, the framework of the MC simulation was completed as used inthese days. As the computer performance improves, the MC simulation of sputtered particle
transport, using their models, has come to applied widely since then.

1.6.3

Applications of MC simulation

Because of the detailed description of the simulation method by Motohiro, and the establishment of the reasonable treatment of low energy region by Turner, many applications have beenreported since the beginning of 1990’s.
Several MC studies were already performed before the Turner’s study, including the discussion of the edge eﬀect of the mask in ion beam sputtering system [41], the evaluation of thedependence of particle energies on hydrogen partial pressure in the reactive sputtering of a-Si:H
using Ar-H2 gas mixture [157], and the study of the energy and angular distributions of arrival
particles at the substrate for Si/Ar sputtering [43]. Elsing has been developed the MC simulation independent of Somekh, Motohiro and Turner in 1991 [31]. His model was almost identicalwith theirs, but simpler assumptions were used for the ejection from the target. It should be
noted that he took into account the inclination of the target normal caused by the formation of
the erosion track, when ejection angle of the sputtered particles were deduced. Some results ofcalculation were shown for Ti/Ar system.
Turner et al. have also been presented some reports themselves. In 1990, they discussed
the dependence of the Cu concentration in YBCO superconducting ﬁlms on gas pressure and8mentioned in 1.3.1

CHAPTER 1. INTRODUCTION

18

T–S distance [155]. In 1992, they chose several elements which have a variety of binding energy(BE) and atomic mass, and simulated pressure dependence of their thickness proﬁle, arrival
angular distribution, mean energy and energy distribution for Ar sputtering [154]. The elements
include: C, Nb, W (high BE), Si, Pd, Au (medium BE), Na, Cd, Pb (small BE). Furthermore,in 1995, they applied the MC simulation to evaluate the gas rarefaction eﬀect [115] proposed by
Rossnagel et al.By splitting the chamber space into meshes and calculating the energy transferfrom the sputtered particle to ambient gases, they obtained and displayed the gas temperature
proﬁle and its dependence on target elements, discharge current and voltage [152].
Nanbu et al. have also used MC simulation and discussed the eﬀect of target power and gas
pressures [99], as well as the pinhole ﬁlling in case of the collimated sputtering [42].
Tuda et al. have applied MC method to interpret the experimental result of the trench ﬁllingwith (Ba, Sr)TiO3 in 1994 [151]. Though the simulation was performed only for Ti/Ar and thesimple hard sphere model was used without the change in its radius by energy, the experimental
results were reproduced fairly well where the side wall coverage in the shallow trench increased
with increase in gas pressure (0.1 → 10 mTorr).
In 1995, Eisenmenger-Sittner et al. have compared the angular distributions of the sputteredparticles at chamber wall obtained from MC simulation and experiment, for the axisymmetricmagnetron sputtering system which has the target inside [30]. They assumed that the gas
atoms were stationary. They reported that they had to modify the angular dependence ofparticle ejection (the n of cosn θ distribution) as the gas pressure changed to ﬁt experimental
data.
Though some of the contents were already mentioned in 1.4.2, we would like to show here
several reports again which combine the ﬁlm growth model on the substrate with the particle
transport process. They have appeared since the midst of 1990’s.
The group of Smy, Dew and others have developed the simulation code, called SIMSPUD/ SIMBAD . The former is the MC simulation treating the particle transport, and the latteris the one to simulate the ﬁlm growth. SIMSPUD has been proposed in 1994, based on the
Somekh’s model. They used this SIMSPUD / SIMBAD, and indicated that good agreements
were obtained with experiments for trench ﬁlling [24] and for collimated sputtering [22]. Lin
and Cale also used their method, and calculated the thickness proﬁle on the substrate which
located behind the hexagonal collimator [70].
Ohta et al. have proposed the simulation method which uses the viscous ﬂow model ofLeon [68] for the atomic diﬀusion process on the substrate, and calculated the thickness proﬁleinside the contact hole embedded in the substrate [102]. Agreement with experiments was goodunder 1 Pa, but became poorer at higher pressures than that.

1.6.4

Extensions of MC simulation

Proposals to extend the MC method have also been reported since the study by Turner. In thissubsection, we would like to review them in the historical order.
In 1991, Myers, Ruzic and coworkers have used their fractal TRIM simulation code [124,125]
to calculate the energy and angular distributions of sputtered particles, instead of using theThompson’s formula and cosine distribution [98]. Van der Waals potential with r −6 dependence
was used for the scattering potential, and the motion of Ar gases were neglected. They applied

CHAPTER 1. INTRODUCTION

19

the simulation to the Si/Ar sputtering with the disk target and the rectangular chamber. As
the gas pressure increased, thermalized particles increased and the high energy ﬂux of Si atoms
decreased at the substrate. They did not discuss the diﬀerence between the simulation withfractal TRIM and that with Thompson’s formula.
Avaritsiotis et al. have applied the MC simulation to the reactive sputtering in 1992 [3].
They incorporated the absorption process of the reactive gas, and deduced deposition rates andﬁlm compositions at various positions on the substrate.
The report by Yamamura et al. in 1995 have also used the simulation inside of the target, toobtain the ejection parameters of sputtered particles [175]. For Cu/Ar sputtering, they applied
their ACAT code [176] and calculated angular / energy proﬁles for both sputtered particles
and the neutralized / recoiled high energy Ar atoms. It was indicated that the result of ACAT
code diﬀered signiﬁcantly from those of Thompson’s formula and cosine rule. For the scattering
potential, they proposed the Thomas-Fermi-Morse potential which was obtained by adding the
exponential term to the Thomas-Fermi potential to ﬁt the Morse potential.
In the same year, Stache has proposed a method to treat the thermalized transport with
diﬀusion equation, which were solved by ﬁnite element method [140] to reduce the computational
complexity. In his method, both target and substrate were assumed to be the inﬁnite plane, andthe thermalized position of particles were recorded by the relative coordinate from the ejection
position. They were superposed with the weight of the depth of the erosion track at the target,
and the proﬁle of thermalized position of particles were obtained. It was treated as the sourceof particle in the diﬀusion equation. The calculation results were compared with experiments
for Si/Ar sputtering, and good agreements were obtained for pressures of 1 Pa and 10 Pa.
Bogaerts have studied the proﬁle of thermalized position of sputtered Cu atoms for Cu/Ar
sputtering in 1995 [12]. He intensively discussed the eﬀect of interatomic potential, and compared
the result using Molier potential with those using power law potential with orders of 1 ∼ 4. At
100 Pa, the peak of the thermalized proﬁle came to about 0.05 cm from the target with Molier
potential, while it went much further in the case of the power law potential with order 4.
In 1996, Malaurie et al. have compared the thickness proﬁle of simulation and experimentfor Cu, Ti, Al, etc. in the axisymmetric sputtering chamber [75]. In their simulation, they
introduced the ellipsoid-like ejection angle distribution:g(θ) = b2α2 cos θ,
1 + (α2 − 1) cos2 θ(1.1)

where b was the major axis of the ellipsoid, and α was the ratio of b with the minor axis a. By
comparing with the experiment, α = 1 (i.e. cosine rule) was optimal for Cu, and the undercosine
values of 1.75 and 2.5 were good for Ti and Al, respectively. These optimal α’s did not dependon the gas pressure nor the T–S distance. In the case of Ti-Al alloy, α = 2.5 reproduced wellthe experiment for both elements.
In the report by Serikov et al. in 1996 [130], they have proposed a method in which theﬂow ﬁeld of the ambient gas is deduced by DSMC9 at ﬁrst, and the MC simulation based on
the Turner’s method was performed. They have developed this model further, and proposed
the particle-in-cell / Monte Carlo (PIC/MC) method in 1997 [131]. In this hybrid method,
9

mentioned in 1.3.3

CHAPTER 1. INTRODUCTION

20

the chamber was split into cells, and the energy brought into them by the sputtered particles
were deduced using MC simulation. Then the energy ﬂow between cells was solved by the ﬁnite
diﬀerent method. This procedure was iterated until the thermal proﬁle converged.
In 1997, Smy et al. have discussed the treatment of sputtered particles redeposited onto theerosion track [137]. In high pressure environment, sputtered particles experience many collision
nearest the target, so the redeposition occurs frequently. They explained that the net of theejected particles, including these redeposited ones, should be proportional to the depth proﬁleof the erosion track, or the match between the simulation and the experiment became worse inhigh pressures. They also showed the result treating the particle transport only by FEM where
the ejection ﬂux of particles were proportional to the erosion depth, but the result did not agreewith the experiment well.
Bogaerts et al. have proposed the method to combine the MC with the FEM treatment of
the diﬀusion process in 1998 [11]. They applied it to Li/Ar. Using laser absorption experiment,
spatial proﬁle of Li atoms were measured. And the simulation was performed by changing the
sticking coeﬃcients on the chamber wall to ﬁt the result with the experiment. However, theagreement between them was rather poor.

1.7

The goal of this study

In previous studies, gas atoms / molecules were treated as immobile, since the sputtered particles
have energies much higher than that of gases initially. Nevertheless, these sputtered particles
are decelerated as low as thermal energy of gases (i.e. “thermalized”) with less than 10 timesof collisions. Hence as the pressure increases and the mean free path becomes shorter, most ofparticles are thermalized before they reach the chamber wall. We believe this is the reason why
the accuracy of the simulation becomes worse in high pressures.
Turner’s model is properly incorporate the eﬀect of the gas motion for mean free path, but
still assumes that the velocity distribution is independent on the speed of the sputter particle. It
is not valid as you can see by considering the coordinate on the sputtered particle in the velocity
space. For example, when the particle speed is high enough, most of the gas collides from theforward of the particle. Since the scattering direction is signiﬁcant in the MC simulation, it isdesirable to treat this eﬀect accurately.
The further pressure increase brings more scattering events, hence much more calculation
time is needed. After the thermalization, however, the motion of the particle is well describedby the random walk, so the ensemble of them can be treated by the diﬀusion equation with thesource term constructed by thermalization points of particles. By combining the MC simulationand the diﬀusive treatment, the calculation time should greatly be reduced.
Summarizing, we would like to develop the models that can be applied to the MC simulation:• We deduce the distribution of colliding gas velocity and its dependence on the particle
speed by considering the kinetic theory of gases.
• We propose the method to account thermalized particles with the diﬀusion equation, andtry to reduce the computational complexity. It will compared with the MC method tocheck the validity, and combined to it after that.

CHAPTER 1. INTRODUCTION

21

The developed model will be applied to the high pressure sputter deposition process. We
will check how well the simulation reproduce the experimental results, and discuss the origin ofsome phenomena observed experimentally in high pressure environment.

Chapter 2TheoryIn this chapter, we discuss physical backgrounds of the particle transport simulation. We skip
the detailed description about the implementation of the simulation program here, which will
be given in appendix A. In the following, we assume that the sputter deposition system isaxisymmetric for simplicity, while our discussion can be extended to 3-dimensional system easily.
In the Monte Carlo simulation of the sputtering process, we assume that the “life” of sputtered particles can be treated as the Markov process1 . The life includes ejection from thetarget, the collision and scattering by the ambient gas atoms/molecules, and deposition on thechamber wall. To execute this simulation, we need probability functions for these steps. We
prepare them (e.g. f (x) for variable x) before the simulation, and use them to determine x withthe pseudo-random number w generated by the computer. Normally, w distributes uniformly
between (0 ≤ w < 1).
Using these functions, we track down 104 ∼ 106 particles, and stochastically obtain the particle ﬂux onto the chamber wall (including substrates), as well as energy and angular distributionsof those particles.
In the ﬁrst section of this chapter, we brieﬂy explain how to determine the parameters (i.e.,x) in each Markov process of the Monte Carlo (MC) simulation. Three sections follow, where thephysical models of MC methods are discussed in the sequence: ejection, collision and scattering,and deposition. These have been proposed by previous works, reviewed in sections 1.6.1 and1.6.2. The last two sections are the discussion of methods we developed in this study. In
section 2.5, we extend the collision model to take the thermal gas motion into account. And in
the ﬁnal section, we discuss the treatment of diﬀusive transport of thermalized particles using
the diﬀusion equation and boundary element method.

2.1

Concept of the Monte Carlo method

In this section, we discuss how to determine the parameter x at each Markov process, fromits given probability function f (x). Two methods are generally used. The ﬁrst one is called“acceptance-rejection” technique. This is generally used when the analytical integration of f (x)1

In this stochastic process, the probability fi , for which the system in status q0 at time t0 becomes qi at t + δt,
does not depend on the history of the system before t0 but only on the status q0 itself.

22

CHAPTER 2. THEORY

23

is impossible. Before the simulation, you have to determine the value fsup , which is larger thanthe maximum of f (x) between xmin ≤ x ≤ xmax . On simulation trials, two independent pseudorandom numbers w1 and w2 are generated between the range [0, 1), and are tested whether theysatisfy the following conditions at the same time:
xr = xmin w1 + xmax (1 − w1 )w2 fsup ≤ f (xr )

(2.1)
(2.2)

If both of them are true, xr is used as the result of this trial. If any of them is false, anothertwo random numbers are generated and tested.
The second method uses the relation: xr
f (ξ)dξ.
(2.3)
w =  xxmin
max
xmin f (ξ)dξ
If f (x) is analytically integrable, and if the primitive function of f (x) can also be solved analytically, you can directly calculate xr from eq. 2.3. It is preferable generally from the viewpointof the calculational complexity.
You can still use the second method even if the analytical root of the primitive function off (x) cannot be obtained. You can solve the eq. 2.3 numerically using Newton method, bisection
method, etc. It’s diﬃcult to determine which method to use in this case. It depends on theprobability of the rejection in “sampling and rejection” method and on the expected number of
iteration required to solve the equation 2.3 with satisfactory precision.

2.2

Particle ejection from the target

Each MC trial starts from the ejection of the sputtered particle. In this section, methods todetermine its position, energy and angle are explained in corresponding subsections.

2.2.1

Ejection position

We consider the axisymmetric system in this study, and the target has disk form. Ejection pointof the sputtered particle from the target is determined so that the probability of the ejection is
proportional to the removed volume of the target there. The removed volume is calculated fromthe depth proﬁle of the target erosion track actually measured after sputtering experiments. If
we measure the depth dT (r) of the target at radius r from the center of the target, probability
density function of the ejection point P (r) is:
P (r) ∝ r dT (r)(2.4)

At higher pressures, re-deposition of the sputtered particles onto the erosion track becomes
more signiﬁcant. In the simulation, these particles should be re-ejected immediately from there,or the depth proﬁle of the simulation will be diﬀer from that of the experiment [137].

CHAPTER 2. THEORY

2.2.2

24

Energy distribution

For the energy distribution function (EDF) of sputtered particles, we use the Thompson’s formula, which has been widely accepted and used in various sputter-related simulation studies.
He actually measured the EDF for several poly-crystalline metal targets, and found that thehigh-energy tail of those EDF’s decreased inversely proportional to the square of the particle
energy. He then proposed this EDF formula to explain these trends. In the following of thissubsection, we try to follow his paper [146] and deduce this formula.
Firstly, we consider the collision cascade occurring in the bulk of the sputter target, andderive the ﬂux of particles through an arbitrary plane. The ﬂux Φ (E  , r  ) with direction r  andenergy E  is:cos θ  dΩ 
dE ,
(2.5)
Φ (E  , r  )dE  dΩ = v  ρ(E  , r  )4πwhere v  is the speed of the particle with energy E  , and θ  is an angle between the normal of
the plane and r  .
There should be particles which act as the source of such a collision cascade. We now
introduce the energy density q(E2 ) of such particles with energy E2 . We also introduce thecollision frequency ν(E2 , E  ) of these, where their energy decreases from E2 to less than E  .
Using these, the spatial density of particles which are decelerated below E  in unit time is given
by:

∞

E

q(E2 )ν(E2 , E  )dE2

(2.6)

The reduction rate of the kinetic energy of the particle at E  in unit time can be expressed
as dE  /dt. This can be rewritten as v  (dE  /dx) using the particle speed v  = dx/dt, hence we
can deduce the relation

 −1

 dE
.
(2.7)
dt = dE v
dx
Substituting 2.6 and 2.7 into 2.5 leads: ∞q(E2 )ν(E2 , E  )cos θ  dΩ dE 
  

dE
(2.8)
Φ E , r dE dΩ =2
dE  /dx
4π
EAt this stage, Thompson introduced several assumptions:• ν(E2 , E  ) = ηE2 /E  (η is a constant with the order of unity)• dE  /dx ∼ E  /D (D is the nearest neighbor distance of the bulk atoms)• q(E2 ) is independent of the location (uniform in space).
• Inter atomic potential between the primary ion and the bulk atom is:
2ER(Z1 Z2 )5/6 (a0 /r)2 ,ewhere ER is the Rydberg energy and a0 is the Bohr radius. Using this,
0(E2 > ΛE1 )q(E2 ) =π 2 a2 nEa Λ1/2 Φ1(E
2 < ΛE1 )
8E 1/2 E 3/2
V (r) =

1

2

(2.9)

(2.10)

CHAPTER 2. THEORY
E1
Φ1
n
a
Λ
Ea

25

energy of the primary ion
ﬂux of ions crossing unit area normal to their pathdensity of atoms
screening radius of interatomic potential (= a0 /(Z1 Z2 )1/6 )energy transfer coeﬃcient between the ion and the atom(= 4M1 M2 /(M1 + M2 )2 )the value of E1 that gives a distance of closest approach of
a in a head-on collision (= 2ER (Z1 Z2 )7/6 (M1 + M2 )/(eM2 ))Table 2.1: notations in equation 2.10

(please refer to table 2.1 for notations).
With these assumptions, equation 2.8 becomes:
2 a2 nE Φ ηDπa1Φ (E  , r  )dE  dΩ =1−
16E  2

EΛE1

cos θ  dΩ
4π

(2.11)

Now we progress into the case where the bulk is cut by the surface, and the potential
diﬀerence (=binding energy) Eb exists between the particles inside/outside of the bulk. Whenthe inside particle with kinetic energy E  , speed v  and angle θ  emits outside of the surface
and has E, v and θ there, these quantities should fulﬁll the following relations because of thepreservation of energy, tangential momentum, etc.:v  sin θ  = v sin θ
1
1
2
M2 v  cos2 θ  =M2 v 2 cos2 θ + Eb
22
E  = E + Eb .

(2.12)(2.13)(2.14)The relation between the density function of the ﬂux inside/outside of the bulk is:Φ (E  , θ  )dE  sin θ  dθ  = Φ(E, θ) sin θdθ,(2.15)

that is,dE  dθ  sin θ .
dE dθ sin θ
Therefore, we can rewrite the equation 2.7 using eqs. 2.12–2.14 and obtain: ∞q(E2 )ν(E2 , E  )
cos θ
dE2 .
Φ(E, θ) =4π(1 + Eb /E) E+Eb
dE  /dx
Φ(E, θ) = Φ (E  , θ  )(2.16)(2.17)

Now we can use the same assumption and procedure, which are used to obtain eq. 2.11 fromeq. 2.7. The result is:
1 − (Eb + E)/ΛE1
dEdΩ,(2.18)
Φ(E, θ)dEdΩ = P cos θ
E 2 (1 + Eb /E)3

CHAPTER 2. THEORY

26

where P = πa2 nΛEa ηDΦ1 /16.
As shown in eq. 2.18, energy distribution is independent of the ejection angle θ. So on applying this to the MC simulation, we determine the energy independently of the angle. Neglecting
the constant term, EDF can be written as:
1 − (Eb + E)/ΛE1(2.19)
f (E) ∝ E
(E + Eb )3
The region of the energy E is 0 < E < ΛE1 − Eb . Since ejection event is basically only at once
for each particle, it does not eﬀect signiﬁcantly on the calculational complexity. So we use theacceptance-rejection procedure to determine E from this function in this study.

2.2.3

Angular distribution

In case of collision cascade as described in the preceding subsection, the particle ﬂux ejecting
through the unit area of the surface obeys the cosine rule:Φ(θ  ) ∝ cos θ  .

(2.20)

Note that it is lead by the assumption that the distribution of particle speed in the collisioncascade is independent of the direction. θ is the angle between the direction of the ﬂux and thesurface normal, and is in the range 0 ≤ θ < π. We assume that the cosine law is valid throughout
this study. Though it should be noted that a few reports have stated the distribution does notagree with the cosine law (as introduced in 1.2.1).
In the simulation of the particle transport process, including ours, you have to determine
the 3D vector of the particle direction on the ejection event. Hence, it is more convenient toknow the distribution of cos θ than that of θ itself. So we try to transform the ﬂux distribution
Φ(θ) into f (x) with x ≡ cos θ. Then we can obtain the relation:dx−f (x) √1 − x20 ≤ θ ≤ π → 1 ≥ x ≥ 0.
√
If the primitive function F (x) of f (x)/ 1 − x2 can be obtained analytically, x is yieldedfrom the pseudo-random number w by solving:
Φ(θ)dθ

ω==F (x) − F (0).
F (1) − F (0)(2.21)

If the inverse function F −1 of F is also given, x is obtained by:
x = F −1 ((F1 − F0 )w + F0 )(2.22)√In the case of cosine distribution, where f (x) is simply x, F (x) becomes − 1 − x2 . Hence,using eq. 2.22, cos θ is obtained from the random number w by:√(2.23)cos θ = 1 − w

CHAPTER 2. THEORY27

At the same time, sin θ is also given with the same w:√
sin θ = w

(2.24)

On ejection, you also need the polar angle φ in addition to the θ. In this study, we simply
assume that it distributes uniformly between 0 ∼ 2π.

2.3

Collision and scattering with gases

In this section, we consider the collision of the sputtered particle with ambient gases. The
treatment of the thermal motions of gases is described in detail in section 2.5, while some of therules in this section are still applicable there.

2.3.1

Mean free path

In this subsection, we treat the “free path” of the sputtered particles. It is the distance inwhich the particle goes without the collision event. In the current discussion, we assume thatthe sputtered particle is traveling through the gas environment, in which each gas is stationary.
If the collision occurs when the distance between the center of the sputtered particle and thatof the gas becomes less than bmax , gases can be treated as targets with the cross section ofσ = πbmax 2 . In this case, when the sputtered particle goes the distance dx, the collision occurs
if the gas is in the volume σ dx. Using the gas density n, the probability of this collision event
becomes nσdx. As a consequence, the distribution p(x) of the free path x of the sputteredparticle obeys the diﬀerential equation:dp= −pnσdx
dx(2.25)

If we deﬁne the quantity λ which has the unit of length as:
λ≡

1
,
nσ

(2.26)

we can write p(x) in the form of Poisson distribution
p(x) =

x
1
exp −
.
λ
λ

(2.27)

Note that this function is normalized. Also note that the λ is the mean value of x since: ∞x p(x)dx = λ.
(2.28)
0This λ is called “mean free path (MFP)”.

CHAPTER 2. THEORY

2.3.2

28

Potential scattering problem

We use the potential scattering for the scattering process between the sputtered particle andthe gas. First of all, we introduce the mass, position and the velocity of these two particles: mp ,mg , r p , r g , v p , v g . Subscript p and g denote the particle and the gas, respectively. Boldfaces
denote the vector quantity. We also introduce the interatomic potential U (r). It is assumed
that U (r) is spherically symmetric (i.e., dependent only on the distance r), and that U (r) → 0
when r → ∞.
We deﬁne the center-of-mass coordinate r G and the relative coordinate r R as:
mp r p + mg r g
mp + mg= rp − rgrG =(2.29)

rR(2.30)

Velocities v G and v R are also deﬁned as r˙G and r˙R , respectively. Then the kinetic equations ofr G and r R are:(mp + mg )r¨G ≡ (mp + mg )v˙G = 0(2.31)

μr¨R ≡ μv˙R = gradU (rR )(2.32)

where we introduce the eﬀective mass μ = mp mg /(mp + mg ). We can see from eq. 2.31 that vG
does not change throughout the collision process.
Since U (r) is central potential, the r R is always included in a single plain which also includes
the origin r = 0. So we assign r R in this plane with polar coordinate (ρ, χ) as in ﬁgure 2.1. The
equations of motion (2.32) then become:
d
U (ρ)
dρ
μ (2ρ̇χ̇ + ρχ̈) = 0.
μ (ρ̈ − ρχ̇) =(2.33)(2.34)

The former gives the conservation of energy:1μ ρ̇2 + (ρχ̇)2 + U (ρ) = const.
2(2.35)and the latter gives the conservation of angular momentum:ρ2 χ̇ = const.

(2.36)b in ﬁgure 2.1 is called “scattering parameter”, which is deﬁned by the component of r R
vertical to v R at inﬁnity. In the MC simulation, the collision event occurs if this is less thanbmax , which was appeared in 2.3.1. If we deﬁne v0 as the absolute value of v R at inﬁnity, theRHS constant terms of eqs. 2.35 and 2.36 are determined:

1
1  2
μ ρ̇ + (ρχ̇)2 + U (ρ) =
μv0 2
2
2
ρ2 χ̇ = v0 b

(2.37)(2.38)

CHAPTER 2. THEORY29(ρ 0 ,χ 0 )
θ

ρ

v0
χ

Scattering Centerb

Figure 2.1: Potential scattering in the center-of-mass system.
We can combine eqs. 2.37 and 2.38:
χ̇
b/ρ2
dχ
= = ±.
dρρ̇
1 − (b/ρ)2 − 2U (ρ)/μv0 2(2.39)

Now we introduce the point (ρ0 , χ0 ) as shown in ﬁg. 2.1 where ρ̇ = 0. Then the integration
between ρ∞ → ρ0 corresponds to 0 → χ0 , so we obtain:
 ρ0
b/ρ2

dρ.
(2.40)
χ0 = −
1 − (b/ρ)2 − 2U (ρ)/μv0 2∞
ρ0 can be given by the relation:1−b 2 2U (ρ0 )
−
=0
ρ0
μv0 2

(2.41)

which is obtained by substituting eq. 2.38 and ρ̇ = 0 into eq. 2.37. It means that the integrationof eq. 2.40 has the singularity at ρ0 . However, as described in A.1.1, this singularity can be
avoided by the numerical integration using the double exponential type transformation.
Finally, the scattering angle θ is given by:θ = π − 2χ0(2.42)as shown in ﬁg. 2.1.
Remember that the ﬁgure 2.1 represents the plane including r R and the origin (hence thev R also), and note that the angle of this plane around the direction of initial vR should be
randomly distributed between 0 ∼ 2π. Also note that the absolute value of v R does not change
since we are considering the elastic collision only. Therefore, if we take the z axis parallel to thedirection of vR before collision, v R is:v R = T (0, 0, v0 ),where T means the transposition of the vector. Then, velocity v R  after the collision becomes:v R  = T (v0 sin θ cos φ, v0 sin θ sin φ, v0 cos θ)CHAPTER 2. THEORY

30

where θ is given from eq. 2.42 and φ = 2πw with the pseudo random number w between [0, 1).
Now we replace r → v in eqs. 2.29 and 2.30 (which also gives the valid expressions) andsolve them. Then we obtain:mgvR .
(2.43)
vp = vG +mp + mg
Equation 2.43 is also valid after the collision, namely v p  and v R  . Remember again that r G
does not change by the collision. If the gas is stationary before the collision (i.e., v g = 0), vG
is:
mp
vp.
(2.44)vG =mp + mg
In summary, the velocity v p  of the sputtered particle after the collision is given by:
vp =

2.3.3

mp v p + mg v R 
mp + mg(2.45)

Scattering potential

In the calculation of eq. 2.40 in this study, we use the Born-Mayer type interatomic potential 2 :U (r) = A exp(−br).

(2.46)

It has a simple form, and the parameters of many elements have been given by Abrahamson [1].
They are also tabulated in appendix B of this thesis.
Abrahamson has only given the potential parameters between the atoms of same element,
but he also showed that the potential between the diﬀerent elements could be approximated.
The potential Ua−b (r) between the atoms of element a and b is:
ba + bb
1/2
1/2
r ,(2.47)Uab (r) = (Ua−a Ub−b ) = (Aa Ab ) exp −2
where Aa and ba are the parameters between a–a, and Ab and bb are those for b–b.
For example, we show the potential between Cu and Ar in ﬁgure 2.2 (a). Using this potential,
the scattering angle θ can be obtained from eq. 2.42, which is shown in ﬁgure 2.2 (b). θ is afunction of the initial energy E0 = (1/2)μv0 2 and the collision parameter b. You can see fromthe ﬁgure that the smaller E0 gives larger θ at constant b, and requires larger b for constant θ.

2.3.4

Collision parameterIn the simulation, you need to choose the maximum collision parameter bmax to calculate the
MFP as described in 2.3.1. The collision parameter b for each event in the simulation can be
determined such that it distributes uniformly in the area of the target σ = πbmax 2 , i.e.:f (b)db ∝ bdb.

2

It means that we assume all the sputtered particles, as well as gases, are atoms. See section 3.1.2(2.48)

CHAPTER 2. THEORY

31

Cu−Ar

2

10

U(r) (eV)Scattering Angle (rad)

3

0

10

−2

10

1

2

3

2

r (angstrom)

E0=0.1eV
1eV

1

10eV
0
0

4

Cu−>Ar(CoM)

2

4

6

Collision Parameter (angstrom)

Figure 2.2: (a) Born Mayer potential between Cu and Ar. (b) Dependence of the scattering
angle on collision parameter in center of mass system, in case of Cu → Ar collision.
Therefore, using the pseudo-random number w between [0, 1), b is determined by:√b = bmax ω.

(2.49)For the potential-scattering case, there is no standard way to determine bmax . It depends on
the threshold you will give for θ (with some E0 ). Larger bmax can treat the subtle change in θ,
while it causes the shorter MFP, hence the larger computation time.
Nevertheless, the reproducibility of the MC calculation is not aﬀected if bmax is large enough.
We show it in the following.
We consider two cases, one is bmax = b1 and the other is bmax = b2 (b1 < b2 ). We assume
that b1 is already large enough, and that collision does not actually occur between b1 < b < b2
in the latter case.
When bmax = b1 , mean free path λ1 is given by eq. 2.26:λ1 =1nπb1 2

And the probability p1 (x)dx where the collision occurs between x and x + dx is:


1
x
exp −
p1 (x) =λ1
λ1(2.50)(2.51)On the other hand, when bmax = b2 , MFP λ2 is:
λ2 =

1
.
nπb2 2(2.52)

CHAPTER 2. THEORY

32

Now we consider the probability p2 (x)dx, where the collision of b < b1 occurs between x andx + dx. It is given by the sum of the probabilities that “collision of b1 < b < b2 occurs n times”
AND “collision of b < b1 occurs between x and x + dx”, from n = 1 to inﬁnity. Hence:







 x
x−ξ
b1 2 b1 2
1
x b1 2
1
ξ
exp −
1− 2
p2 (x) =exp −+dξ 2 exp −λ2
λ2 b2 2
λ2
λ2
λ2
b2
b2 2
0






2 2
 x
 x
1
b1
ξ1
ξ2 − ξ1x − ξ2b1 2+
dξ1
dξ2 3 exp −
exp −
exp −
1− 2
λ2
λ2
λ2
λ2
b2
b2 2
0
ξ1
+ ....
(2.53)
Note that the terms of exp’s in the integral can simply be rewritten as exp(−x/λ2 ), so using therelation
b1 2
λ2
2 = λ
b2
1
 x

 x

and

0

dξ1

ξ1

dξ2 · · ·

 x
ξn−1

dξn =

xn
,
n!

p2 (x) becomes:
p2 (x) =
=







1
x
1
1
λ2
λ2 λ2
λ2 2 λ2 x2+ . . . ..
exp −+x+ 2 1−
1−
λ2
λ2
λ1 λ2
λ1 λ1
λ1
λ1 2!
λ2




 

x
1
1
x
1
λ1
1
exp −
−
=
exp −
exp
x
,
(2.54)
λ2
λ2
λ2 λ1
λ2
λ1
λ1

and is proved to be identical with p1 (x).

2.4

Trapping on the chamber boundary

In this study, we simply assume that the particle reached at the chamber boundary deposits therewith the probability of unity. The estimation whether the particle path crosses the boundary israther complicated in the axisymmetric coordinate, hence we summarize the procedure here.
We use the index p for the position of the particle, and b for the boundary. The componentsof the 3D Cartesian coordinate is expressed by x, y and z, and its projection to the axisymmetriccoordinate by r and z. With these, the following relation comprises:r 2 = x2 + y 2(2.55)

After a collision event, the direction and the free path of the particle is determined. The
start and the end of this path are deﬁned as (xp1 , yp1 , zp1 ) and (xp2 , yp2 , zp2 ). In the following, we
consider the method to determine whether this path crosses with any of the (linear) boundary
element from (rb1 , zb1 ) to (rb2 , zb2 ).

CHAPTER 2. THEORY

33

At ﬁrst, we consider necessary conditions of the trap event. The path of the particle is:
x = (1 − t)xp1 + txp2
y = (1 − t)yp1 + typ2z = (1 − t)zp1 + tzp2
where (0 ≤ t ≤ 1). By projecting these to the axisymmetric coordinate, we obtain:22+ t2 rp2+ 2t(1 − t)(rp1 · rp2 )r 2 = (1 − t)2 rp1(2.56)z = (1 − t)zp1 + zp2 ,(2.57)

where (rp1 · rp2 ) is deﬁned by:(rp1 · rp2 ) = xp1 xp2 + yp1 yp2 .

(2.58)

Hence the path of the particle becomes hyperbola in (r, z) plane.
Substituting eq. 2.57 into eq. 2.56 and eliminating t, we obtain:22 − (r · r )
 2

rp1
p1
p2
2
2
rp1 + rp2 − 2(rp1 · rp2 )t− 2r =2 − 2(r · r )rp1 + rp2p1
p2+

2 r 2 − (r · r )2
rp1
p1
p2
p2
2 + r 2 − 2(r · r )
rp1
p1
p2
p2

Hence, if0≤.

2 − (r · r )
rp1
p1
p2
2 + r 2 − 2(r · r )
rp1
p1
p2
p2

(2.59)≤1(2.60)

is true, the minimum rpmin of r component through the path is given by:
2
=
rpmin

2 r 2 − (r · r )2
rp1
p1
p2
p2
2 + r 2 − 2(r · r )
rp1
p1
p2
p2

.

(2.61)

Otherwise, rpmin is the smaller one of rp1 or rp2 . As well, rpmax , zpmin and zpmax are:rpmax = max(rp1 , rp2 )zpmin = min(zp1 , zp2 )zpmax = max(zp1 , zp2 )respectively. We also prepare the deﬁnitions of rbmin , etc. for the boundary node (rb1 , zb1 )(rb2 , zb2 ).
Using these parameters, the necessary condition that the the particle path crosses the boundary is the AND’s of:
rpmax ≥ rbmin
rpmin ≤ rbmax
zpmax ≥ zbmin
zpmin ≤ zbmax(2.62)

CHAPTER 2. THEORY

34

Next, we consider suﬃcient conditions when eqs. 2.62 are true. For this, all of:22+ t2 rp2+ 2t(1 − t)(rp1 · rp2 )r 2 = (1 − t)2 rp1(2.63)z = (1 − t)zp1 + tzp2

(2.64)r = (1 − s)rb1 + srb2

(2.65)z = (1 − s)zb1 + szb2

(2.66)

must be satisﬁed. Moreover, (0 ≤ t ≤ 1) and (0 ≤ s ≤ 1) are required. When such t’s are given
from several nodes, the node which gives the smallest t is the one where the particle is trapped
actually. We will discuss how to obtain such t’s for several cases where the ends of the node
have diﬀerent conditions.
CASE1: rb1 = rb2
From equation 2.63, we obtain:
 2

2
2
2
2
(rp1
+ rp2
− 2(rp1 · rp2 ))t2 − 2 rp1− (rp1 · rp2 ) + (rp1
− rb1
) = 0.

(2.67)

If the discriminant D of this quadratic satisﬁes the condition D ≥ 0, then check:• Solution t satisﬁes 0 ≤ t ≤ 1.
• Substituting the t into eq. 2.64, and the resultant z satisﬁes zbmin ≤ z ≤ zbmax .
If both of these conditions are satisﬁed, the particle is trapped there. When both solutions of
the quadratic satisﬁes these, the smaller one is adopted.
CASE2: zb1 = zb2
If zp1 = zp2 is also satisﬁed, the path of the particle becomes the line parallel to the r axis,
and the trapping point r is:• r = rbmin if rpmin ≤ rbmax .
• r = rbmax if rpmax ≥ rbmin .
If zp1 = zp2 , t becomes:t=zb1 − zp1
zp2 − zp1

(2.68)

from eq. 2.64. This t automatically satisﬁes 0 ≤ t ≤ 1 from eq. 2.62. So if r, obtained by the tand eq. 2.63, satisﬁes rpmin ≤ r ≤ rpmax , the particle is trapped.

CHAPTER 2. THEORY35

CASE3: the other
Eliminating r, z and s from eq. 2.63 gives:
rb2 − rb1 2
2
2(zp2 − zp1 )2 − rp1− rp2+ 2(rp1 · rp2 ) t2
0 =zb2 − zb1
rb2 − rb1 2
zb2 − zb1
2+
(zp2 − zp1 ) zp1 − zb1 +rb1 + rp1
− (rp1 · rp2 ) 2t
zb2 − zb1
rb2 − rb1
 
2
zb2 − zb1
rb2 − rb1 2
2
rb1 − rp1
zp1 − zb1 +.
(2.69)+
zb2 − zb1
rb2 − rb1
If the discriminant D of this quadratic satisﬁes D ≥ 0, then solve the equation and determinet. And then check:• Solution t satisﬁes 0 ≤ t ≤ 1.
• s deduced from eqs. 2.64 and 2.66 satisﬁes 0 ≤ s ≤ 1.
If both of these are satisﬁed, the particle is trapped there. If both of the solution of the quadraticsatisﬁes these, the smaller one is adopted.

2.5

MC calculation with the thermal motion of gases

In this section, we consider the particle p traveling in the gas atmosphere with the speed vp ,and discuss its interaction (i.e., collision and scattering) with environmental gas particles. We
propose the method to simulate this process, which can be applied to the Monte Carlo simulation.

2.5.1

Collision frequency and mean free path

We consider the particle p traveling in the gas atmosphere, and show how the gases collide to
the particle. The discussion in this section is partly based on the chapter X of the reference [56]by Jeans.
We assume that the distribution of the gas speed is the Maxwellian:
4 3
f (vg ) = √ α 2 vg 2 exp (−αvg 2 )π(2.70)

where vg is the speed of the gas and α = mg /2kTg . k is the Boltzmann’s constant, mg is massof the gas particle, and Tg is gas temperature. We also assume that the distribution is only
dependent on the speed, and not on the direction.
Now we consider the contribution of gases from [vg ∼ vg + dvg ], [θ ∼ θ + dθ], [φ ∼ φ + dφ] tothe collision. θ and φ are the azimuthal and the polar angles of the velocity of g corresponding
to the direction of p, respectively. Then the contribution dν from this part of the velocity-space
to the total collision frequency ν is given by:dν = πbpg 2 vR ng f (vg )dvg

sin θdθdφ
4π(2.71)

CHAPTER 2. THEORY

36

where bpg is a maximum collision parameter, ng is gas density. vR is the relative speed between
the particle p and colliding gas g, which is deﬁned by (see ﬁgure 2.3):
2
vR= vp2 + vg 2 − 2vp vg cos θ.

VR(2.72)Vg

Vp

θ

Figure 2.3: vp , vg , vR and θ.
ν is obtained by integrating 2.71 within whole velocity space.
 ∞  π  2π
dν
ν=
vg =0

θ=0

(2.73)

φ=0

Assuming that the density and the temperature of the gas are uniform, integration with φ simply
becomes 2π:
 ∞
 π
1
ν = 2π 2 bpg 2 ng α3/2
dvg
sin θdθvR vg 2 exp (−αvg 2 ).
(2.74)
00

We now transform the integration variable from (vg , θ) to (vg , vR ). Range of the vR thenbecomes from |vp − vg | to vp + vg , while the Jacobian is: ∂vg ∂vg ∂θ
∂(vg , θ)vR ∂vg ∂vR =  ∂θ
.
(2.75)=
=∂θ∂(vg , vR )  ∂vg ∂vR  ∂vRvp vg sin θ
Using these, equation 2.74 can be rewritten as: ∞ vp +vgvg
ν=Advg exp (−αvg 2 )vR 2 dvRvp0|vp −vg |(2.76)

1

where A = 2π 2 bpg 2 α3/2 .
Integration with vR is given as follows:
 2 3
 vp +vg
2(vg > vp )2
3 (vp + 3vp vg )vR dvR =2
3
2|vp −vg |3 (vg + 3vg vp ) (vg < vp ).

(2.77)

CHAPTER 2. THEORY

37

By substituting eq. 2.77 into eq. 2.76, we obtain:
∞
1 vp 4

2
2 2
3
2
2
(vp vg + 3vg ) exp (−αvg )dvg +(vg + 3vp vg ) exp (−αvg )dvg ,
ν=A
vp 0
vp

(2.78)

where A = 23 A.
The ﬁrst term of equation 2.78 is:
 ∞
 ∞
2
2
vg exp (−αvg )dvg + 3
vg 3 exp (−αvg 2 )dvg
vp
vp

vp

∞
vg 2 exp (−αvg 2 )=vg exp (−αvg )dvg − 3
α
2α
vp
vp


∞


∞
2
2
2
1 exp (−αvg )3vg exp (−αvg )
−
=
vp2 +
−3
α
2
α
2α
vp
vp


2
exp (−αvp )3
2vp 2 +
.
=
α
2α


3
vp2 +

 ∞



2

(2.79)

Likewise, the second term is: vp
1 vp 42
vg exp (−αvg )dvg + 3vp
vg 2 exp (−αvg 2 )dvg
vp 0
0
v


  vp1
3
vg 3 exp (−αvg 2 ) p
=
+
+ 3vp
vg 2 exp (−αvg 2 )dvg
−
vp
2α
2αv
p
0
0

v
 

 vp
2
2
3
vp exp (−αvp )1
vg exp (−αvg 2 ) p
2
+
= −
+ 3vp
+
exp (−αvg )dvg
−
2α
2αvp
2α
2α 0
0
√


 αvp
3
1
vp 2 exp (−αvp 2 ) 3 exp (−αvp 2 )
−
+
+ 3vp
exp (−y 2 )dy, (2.80)
= −2
α
4α2
2αvp
2α3/2 0√where y is deﬁned by y = α vg . Combining 2.79 and 2.80 gives:
3 exp (−αvp 2 )
+
4α2



3
+ 3vp
2αvp



1
2α3/2

Now we deﬁne the function Ψ:
Ψ(x) = x exp (−x2 ) + (2x2 + 1) √αvpexp (−y 2 )dy.

(2.81)0 x
exp (−y 2 )dy√0πErf (x).
= x exp (−x2 ) + (2x2 + 1)
2
Erf(x) is the Gaussian Error function which is deﬁned by: x2
exp (−y 2 )dy.
Erf (x) = √π 0(2.82)(2.83)

CHAPTER 2. THEORY

38

The value of this function is in between zero to unity.
By using 2.82, equation 2.81 becomes:√
3
Ψ(vp α).
5/24vp α(2.84)

Finally, the collision frequency ν and the mean free path λ are obtained as follows.
√3Ψ(vp α)
5/2
4vp α

√
2 πbpg 2 ng kTg
Ψ vp
mg vp

ν = A
=

λ≡

vp
ν

mg
2kTg



=

mg vp 21√
2
2 πbpg ng kTg Ψ vp mg
2kTg=

mg vp 21√
2
2 πbpg Pg Ψ vp mg
2kTg(2.85)(2.86)

Pg = ng kTg is the gas pressure.

Mean Free Path (m)

10-2

10-3
Gas: Ar
Pg = 1Pa
Tg = 400 K
bpg = 0.6 nm
10-4 110102103Particle Speed (m/s)

104

Figure 2.4: Dependence of the particle mean free path on its speed.
Figure 2.4 shows the dependence of the particle MFP on its speed, calculated by eq. 2.86.
It is the case for an Ar gas of 1 Pa and 400 K. The latter condition leads the mean speed of Ar
atoms about 460 m/s. The maximum collision parameter is set to 0.6 nm. It is rather larger

CHAPTER 2. THEORY

39

than the one used in the hard-sphere model, since we use the potential scattering so that eventhe large collision parameter may cause the deﬂection of the particle which cannot be neglected
at lower energies (see ﬁgure 2.2). The value of 0.6 nm is chosen in this study so that the cosineof the scattering angle θs in equation 2.42 is less than 10−5 when the relative kinetic energy
2 /2 is 10−3 eV.
μvR
The MFP dependence on the particle speed vp can be separated into two regions by the
relation between vp and the thermal speed of gases. When vp is much faster, gas particles canbe regarded as stationary targets. In this case, the MFP becomes the constant value which has
been used in the MC simulation ignoring gas motion. On the other hand, when vp is muchslower, the particle can be taken as immobile. In this region, the collision frequency approaches
to constant, hence the MFP becomes proportional to vp . From ﬁgure 2.4, it is clearly understoodthat ignoring the gas motion leads to the overestimation of MFP in lower energy region.

2.5.2

Determination of the velocity of colliding gas

In this subsection, we describe how to determine the colliding gas speed vg and the relative speedvR from the pseudo-random number w, based on the discussion in the previous subsection. Oncevg and vR are given, collision angle θ (see ﬁg. 2.3) can easily be deduced.
We assume that the particle p has the speed vp just after some collision (or the ejection fromthe target). Firstly, we calculate the mean free path λ with eq. 2.86, and then the free path ofp using the method described in 2.3.1.
Now we deﬁne the function ν  of the variable vg by (see equation 2.76): vg

ν (vg ) = A
0

vg
dvg
exp
vp

2
−αvg

 vp +vg|vp −vg |

2
vR
dvR .

(2.87)

With this, speed of the colliding gas vg can be determined by the condition:w = ν  (vg )/ν

(2.88)

where w is a pseudo random number in (0, 1].
Once vp and vg are given, another w is generated and used to determine vR , from the following
relation: vR
2|v −v | v dv
w =  vpp+vgg
2|vp −vg | v dv=

3 − |v − v |3
vRpg
(vp + vg )3 − |vp − vg |3(2.89)In the simulation program, distribution function of vg , which depends on vp and random
number w, is approximated by two dimensional spline function so that the calculational complexity is reduced. The detail of the program used in this study will be described in appendix A.
The examples of the vR distribution are shown in ﬁgure 2.5. In these, environmental gas is
assumed to be an argon with the temperature of 400 K. Each subﬁgure represents the cases of
vp = 300 m/s, 1000 m/s and 3000 m/s, respectively. Z-axis of these ﬁgures are taken to be theCHAPTER 2. THEORY

40

Vp = 300m/s

Vp = 1000m/s

VRZ(m/s)

Vp = 3000m/s

VRZ(m/s)

VRZ(m/s)

2000

2000

0

0

2000
0

-2000

-2000

-2000

2000
-2000

0

VRX(m/s)

2000

0
-2000

2000
-2000

VRY(m/s)

0

VRX(m/s)

2000

0
-2000

2000
-2000

VRY(m/s)

0

VRX(m/s)

2000

0
-2000

VRY(m/s)

Figure 2.5: Dependence of vR distribution on vp .
direction of the particle p. vR is randomly determined 5000 times using the eqs. 2.87, 2.88 and
2.89, and plotted.
Centers of distribution seem to be at about (0, 0, vp ), but they are not actually. As can beseen in eq. 2.71, the distribution of vR is given by the multiplication of the Maxwellian and therelative speed vR . Hence, if vg is equal, the probability of the collision becomes higher wherevR is large. In other words, the particle p tends to be collided more by the gases whose velocity
is opposed to it. Also note that the colliding gas velocity cannot be determined simply by theMaxwellian, which was the method used in the report of Turner et al. in [153].

2.5.3

Particle velocity after the scattering

Since our treatment includes the motion of gas particles, the determination of the particle
velocity after the scattering is slightly complicated. Since we prepare this subsection to discussit.
First of all, the following two deﬁnitions are introduced.
1. The velocity vector of the sputtered particle p before the collision has the form:
vp = vp (sin θp cos φp , sin θp sin φp , cos θp )(2.90)

in the laboratory system. vp is an absolute value of the speed. θp and φp are the azimuthal/polar angle of vector vp in the polar coordinate.
2. Three dimensional rotational transforms RY (θ) and RZ (φ) are deﬁned as follows:⎛⎞cos θ 0 sin θ
0 1
0 ⎠RY (θ) = ⎝− sin θ 0 cos θ⎞
cos φ − sin φ 0
cos φ 0 ⎠ .
RZ (φ) = ⎝ sin φ00 1(2.91)⎛(2.92)

CHAPTER 2. THEORY

41

RY (θ) denotes the clockwise rotation around (0, 1, 0) direction by θ, and RX (φ) does around(0, 0, 1) by φ.
Using these, following relation can be deduced:
RY (−θp )RZ (−φp )T vp = (0, 0, vp )(2.93)

where T means the transposition of the vector.
Now we consider the following coordinates:
I The laboratory system.
II The coordinate whose z axis is parallel to the velocity of the particle p.
II’ The coordinate that includes the gas velocity vector within its xz plane, in addition to thecondition of II.
III The coordinate whose z axis is parallel to the relative velocity vector between p and g, andthe velocity vectors of both p and g are included in the xz plane.
In the following, we will obtain conversion rules between these coordinates. We use the format

v I to denote the Cartesian representation of the vector in the coordinate I, and likewise for II,
II’ and III. We also deﬁne the operator TI→II , etc., which transform these vector representations:v (II) = TI→II v (I)(2.94)

Remember that the velocity of the particle p in coordinate I is given by 2.90. As apparently
from equation 2.93, one of the operations which makes this vector parallel to the z axis is givenby:(2.95)
TI→II = RY (−θp )RZ (−φp ).
Hence the reverse operation is:
TII→I = RZ (φp )RY (θp ).

(2.96)In the coordinate II, the colliding gas speed vg and the relative speed vR can be determinedusing the method described in 2.5.2. Now we introduce the velocity vectors vp and vg forparticles p and g, respectively, and vR for the relative velocity between them. If we express
these vectors in II, we obtain:
vp (II) = vp (0, 0, 1)(2.97)

vg(2.98)(II)= vg (sin θg cos φg , sin θg sin φg , cos θg )vR (II) = vR (sin θR cos φR , sin θR sin φR , cos θR )(2.99)

θg and θR are given so as to fulﬁll the following three conditions simultaneously:
vp 2 + vg 2 − vR 22vp vg2
vp + vR 2 − vg 2=
2vp vR= −vR sin θRcos θg =(2.100)cos θR(2.101)

vg sin θg(2.102)

CHAPTER 2. THEORY

42

φg and φR are the identical value, which is obtained by the random number uniformly distributing within 0 ∼ 2π.
Using θR and φR , we can deﬁne the TII→II , TII →III and their inverses.
TII→II= RZ (−φR )(2.103)

TII →II = RZ (φR )(2.104)

TII →III = RY (−θR )(2.105)

TIII→II= RY (θR )(2.106)Note that the components of vp , vg and vR in coordinate II’ are:
v p (II ) = vp (0, 0, 1)
vg(II )

vR(II )and that of vR in III is:(2.107)= vg (sin θg , 0, cos θg )(2.108)= vR (sin θR , 0, cos θR )(2.109)v R (III) = vR (0, 0, 1).

(2.110)In this coordinate III, we can obtain the relative velocity vector vR  after the potential
scattering (We use ´ for the values after the scattering) by using the method described in 2.3.2.
Once vR  is yielded, we simply trace back the conversions.
Summarizing, the velocity vector vp  of the particle p after the collision is, in the laboratory
system, determined by solving:
v p

(I)

(II)
v p
(II )v pv R(II )= TII→I vp 
=

(II)


(II )
TII →II vp= vG(II )

mg(II )+v R
mp + mg= TIII→II v R(III)

(2.111)(2.112)(2.113)(2.114)

where vG (II ) is given using θg appeared in equation 2.100:

mg vg sin θg
mp vp + mg vg cos θg(II )=, 0,vG
mp + mg
mp + mg

2.6(2.115)

Particle transport after the thermalization

Until the previous section, we have discussed the Monte Carlo simulation of the particle transport. As pressure increases, however, computational complexity of the MC calculation becomesvery large. For example, under the condition of 10 Pa and 300 K, if you set the maximum collision parameter to be 7 Å, mean free path λ (ignoring the gas motion) becomes 2.69 × 10−4 m,
which is much smaller than the size of the ordinary vacuum chamber.

CHAPTER 2. THEORY

43

After sputtered particles lose their initial energy and are thermalized 3 , their motion can be
described by the 3-dimensional random walk. In this stage, the number of collisions before thedispersion of its position spreads x is approximately (x/λ)2 . Hence, the calculation amount of
the MC simulation increases proportional to the square of the pressure.
The random walk type transport can be described with the diﬀusion equation, if the mean freepath is much smaller than the chamber size. On applying it to the sputter deposition problem,
we have to consider only the steady state (time-independent) condition, where the diﬀusion
coeﬃcient is not necessary as discussed later. In this section, we discuss the description of thediﬀusion process of the sputtered particles using Poisson equation, based on the idea reported
by Petrov et al. [109].

2.6.1

Diﬀusion equation

The diﬀusion equation of particles in time-dependent form is:∂η(x, t)= DΔη(x, t) + ρ(x, t),∂t(2.116)

where η is spatial density of particles, D is the diﬀusion coeﬃcient, and ρ is a distributionfunction of particle sources. In the current problem, ρ can be regarded as the thermalization
points of the sputtered particles where they start the diﬀusive transport. Since ∂n/∂t = 0 in
the steady state, equation 2.116 becomes the Poisson’s equation.
If we assume that the sticking coeﬃcient of the particle at the chamber boundary is unity,
the ﬂux from outside of the boundary is always zero regardless of the value of particle density
η near the wall. Hence, we can use the boundary condition for ordinary chamber boundary:η = 0.

(2.117)The value to be calculated is the particle ﬂux F onto the boundary:F = Dgrad η · n,

(2.118)where n is the normal vector of the boundary surface. Since we have the partial diﬀerential equation of second order with Dirichlet boundary condition, F at each boundary can be determinedby solving the equation.

2.6.2

Relationship with MC simulation

In this subsection, we consider the case where we trace NM C particles with MC simulation. By
this MC simulation, we see the life of particles; the ejection from the target, the collision andscattering with the gas atoms, etc. The particles may be reached at the chamber boundary
before they thermalized in the chamber. For these particles, we obtain the arrival number ofthem on each chamber boundary segments, as well as their angular/energy distributions. With
the MC calculation, we also obtain the distribution of the thermalized positions of the particles,
as described above. We used it as the source term of the Poisson’s equation.
3

see 1.3.1 for the deﬁnition of “thermalization”

CHAPTER 2. THEORY

44

Now we consider the time TM C corresponding to the trial of NM C particles in actual deposition experiment. If we express the thermalized position of the particle i by xTi h , the source
term ρ in the Poisson’s equation can be written as:
ρ(x) =

1TM C

δ3 (xTi h − x)(2.119)δ3 (xTi h − x) = 0(2.120)

Thermalized

And the diﬀusion equation in steady state becomes:
DΔη +



1
TM C

Thermalized

We must treat two kinds of particles in a same time scale: the ﬁrst includes the particles
which reached the chamber wall during the Monte Carlo step (i.e., before the thermalization),and the second includes the particles arrived at the wall via diﬀusive treatment. Now let us
consider the node Sj of the chamber boundary. The number of particles of the ﬁrst kind, whichreach at this boundary during the time TM C , can be obtained directly from the MC calculation.
We write it as NjEn . On the other hand, the number of diﬀusive particles reached at Sj during
TM C is:(2.121)
NjT h = TM C Sj F = TM C Sj Dgrad η · n.
If we deﬁne θ as:θ ≡ TM C Dη,eq. 2.120 becomesΔθ +

δ3 (xTi h − x) = 0,(2.122)(2.123)

Thermalized

and eq. 2.121 becomes:

NjT h = Sj gradθ · n,

(2.124)

which gives the number of particles NjT h deposited at the node j after the thermalization.
In above, we described that the boundary condition of the chamber wall is Dirichlet type;
η = 0. However, it causes a problem if we apply this condition to the erosion track of thetarget. As discussed in 2.2.1, the particle deposited on the erosion track must immediately beejected again to reproduce the depth proﬁle of the erosion track, which would be obtained by
the experiment.
In this study, we use the Neumann condition (∇η = 0) at the boundary nodes representing
the erosion track. Ideally, we should record the re-deposition proﬁle of the diﬀusive particles,
feedback it to the MC calculation, and iterate this procedure until the erosion proﬁle converges.
This strategy, however, will make the calculation procedure too complicated, so we decided touse the former simpliﬁed model.
Summarizing, the problem in this section is to solve the equation 2.123, under the boundary
conditions ∇η = 0 (Neumann) for the erosion track, and η = 0 (Dirichlet) for others. The
solution is the ﬂux of the particles (in eq. 2.124) at the boundaries where the Dirichlet condition
was applied.

CHAPTER 2. THEORY

2.6.3

45

Boundary integral equation

In the following part of this section, we discuss the boundary element method (BEM) to solve theproblem submitted at the end of the previous subsection. The BEM method has been proposedby Brebbia and Wrobel [14]. We would like to note that the doctoral thesis of Wrobel [170] hasbeen particularly consulted on applying this method to the axisymmetric coordinate.
We show the equation to be solved again:∇2 θ(x) +h(x)= 0 : inΩk
θ = 0 : atΓ

(2.125)∇θ · n = 0 : atΓ
where Ω denotes the considering volume (3D), Γ and Γ denote its boundaries (2D). To convert
this problem into boundary integral equation, we introduce the function θ ∗ called “fundamental
solution”, which satisﬁes the following relation in Ω:k∇2 θ ∗ (x, y) + δ3 (x − y) = 0(2.126)

In eq. 2.126, x and y are points in Ω. δ3 is the Dirac’s 3D delta function. Operator ∇2 may
eﬀect either x or y. The solution of eq. 2.126 is:θ ∗ (x, y) =

1
,
4πrk

where r = |x − y|.
Next, we integrate the equation 2.125 using this fundamental solution in Ω, e.g.: h(x )2θ ∗ (x, x )dΩ(x ) = 0.
∇ θ(x ) +k
Ω

(2.127)(2.128)dΩ(x ) means that the integration is performed with x . Since the term in the parenthesis isalways zero in Ω, this equation is always satisﬁed.
If we execute the partial integration twice and use the relation in eq. 2.126, we obtain:



1
θ(x)
 ∗


 ∗


+ θ(x )q (x, x )dΓ(x ) − q(x )θ (x, x )dΓ(x ) = −h(x )θ ∗ (x, x )dΩ(x ), (2.129)
k
k Ω
Γ
Γ
∂
∂ ∗
∂
θ and q ∗ = ∂n
θ . ∂ndenotes the derivative by x along the boundary normal
where q = ∂nwhich heads for the outside of the chamber. In other words, it is:∂= n · ∇ .
(2.130)∂nwhere n is a surface normal heading for the outside of the boundary. Equation 2.129 meansthat once θ and q are determined on the boundary Γ, θ at any x in Ω can also be determined.
Now we apply the equation 2.129 for xs on the boundary. In this case, θ ∗ and q are singular
at xs = x , so we must introduce the coeﬃcient c(xs ):θ(xs )
1
 ∗


 ∗


+ θ(x )q (xs , x )dΓ(x ) − q(x )θ (xs , x )dΓ(x ) = −h(x )θ ∗ (xs , x )dΩ(x )c(xs )kkΓΓΩ
(2.131)CHAPTER 2. THEORY

46

If the boundary at xs is smooth, c(xs ) is 1/2. If not, c(xs ) is given by 1/2π times the solid angleof Ω from xs .
Equation 2.131 gives the restriction to θ and q on the boundary in functional sense. In the
problem of the second partial derivative equation, usually the value of either θ or q is given on
the whole boundary, and the counterpart will be derived using equation 2.131. When both θand q are determined on the whole boundary in such way, the value at any point inside theregion can be derived using the equation 2.129.

2.6.4

Axisymmetric caseIn the problem we are considering in this study, the system is axisymmetric. So we rewrite theequation 2.131 here using the polar coordinate. The radial, polar, and azimuthal coordinates ofpoint x are denoted as R(x), φ(x) and Z(x), respectively.
θ(xs )c(xs )
k



 2π

+

θ(x)


0

Γ̄

−

 2π
q(x)


=

q ∗ (xs , x )dφ(x )R(x )dΓ̄(x )θ ∗ (xs , x )dφ(x )R(x )dΓ̄(x )0Γ̄h(xs )kΩ̄ 2π

θ ∗ (xs , x )dφ(x )R(x )dΩ̄(x ).

(2.132)0

in which only θ ∗ and q ∗ are dependent on φ.
When we express the points x and x in polar coordinate, the distance r between these twopoints is:
(2.133)r = R(x)2 + R(x )2 − 2R(x)R(x ) cos{(φ(x) − φ(x )} + Z(x) − Z(x ) ,Hence, from eq. 2.127, the integration is:
θ̄

∗

 2π
=

θ ∗ (xs , x )dφ(x )

0

=

4K(m)
√
,
a+b

where K(m) is the complete elliptic function of 1st kind, and m is:2ba+b
a = R2 (xs ) + R2 (x ) + [Z(xs ) − Z(x)]2m =b = 2R(xs )R(x )Note that m has the value between −1 ≤ m ≤ 1.

(2.134)

CHAPTER 2. THEORY

47

Similarly, we also obtain: 2π∗¯q ∗ (xs , x )dφ(x )
q =
0

 2

1
R (xs ) − R2 (x ) + [Z(xs ) − Z(x )]2
4√
E(m) − K(m) nR (x )
=
)
2R(x
a
−
b
a+b

Z(xs ) − Z(x )
E(m)nZ (x ) ,(2.135)+a−bwhere E(m) is the complete elliptic function of 2nd kind, nR and nZ are R and Z componentsof the normal unit vector at the boundary x , respectively.
Substituting eqs. 2.134 and 2.135 into eq. 2.132, we ﬁnally obtain the expression:θ(x )q¯∗ (xs , x )R(x )dΓ̄(x )c(x )θ(xs ) +Γ̄h(x ) ∗ ∗θ̄ (xs , x )R(x )dΩ̄(x ). (2.136)
q(x )θ̄ (xs , x )R(x )dΓ̄(x ) +
=
k
Γ̄
Ω̄

2.6.5

Boundary element method

Boundary element method (BEM) is a technique to solve the boundary integral equation inwhich you discretize the boundary into nodes, approximate θ̄ and q̄ on it with simple functions
including a few parameters, build and solve simultaneous linear equations of these parameters,
and determine the approximate solutions for θ̄ and q̄. With BEM, we ﬁrst split the boundary Γ
into nodes. In this study, we assume that each node has line form for simplicity.
For the “simple functions” on the nodes, there are some candidates depending on the numberof parameters, e.g.:• constant element (with one parameter)• linear element (with two parameters)• parabolic element (with three parameters)
Note that the “element” denotes the type of approximate function deﬁned on each boundary
node, and has nothing to do with the form of the node. Also note that you can select theconstant element for θ̄ while the linear one for q̄, and the like. In the following of this study, we
use the constant element for both θ̄ and q̄.
We split the boundary Γ̄ into n linear nodes, and denote them as Γ̄i (i = 1 . . . n). The both(i)(i)
ends of the Γ̄i are denoted by x1 and x2 . 1 and 2 are chosen such that the region Ω̄ comes(i)(i)
left side of the vector x2 − x1 . As the representative point of Γ̄i , we choose:
(i)

x

(i)

(i)

x + x2= 1
2On Γ̄i , θ̄ and q̄ are approximated by the constant value θi and qi , respectively.

(2.137)

CHAPTER 2. THEORY48

The integration on the boundary (eq. 2.136) can be split into these nodes:(i)θj q¯∗ (x(i) , x )R(x )dΓ̄(x )c(x )θi +
j

=

Γ̄j


j

qj θ¯∗ (x(i) , x )R(x )dΓ̄(x ) +Γ̄j

h(x ) ¯∗ (i) θ (x , x )R(x )dΩ̄.
kΩ̄

Since θj and qj can be moved outside the integration, equation 2.138 becomes:


∗
∗
qij
θj =
θij
q j + Bici θi +
j(2.138)(2.139)jwhereci = c(x(i) )∗q¯∗ (x(i) , x )R(x )dΓ̄(x )
qij =
Γ̄j

∗
θ¯∗ (x(i) , x )R(x )dΓ̄(x )
θij =

Bi =

(2.140)(2.141)(2.142)Γ̄jh(x ) ¯∗ (i) θ (x , x )R(x )dΩ̄.
kΩ̄

(2.143)In the current problem, either θi or qi is known on every node. Hence eq. 2.139 can be regarded
as the simultaneous linear equations.
Now we deﬁne:∗
ci + qij(i = j)Hij =∗(i = j)
qij
Gij

∗
= θij
.

(2.144)Using these, equation 2.139 can be expressed in the matrix form:
HΘ = GQ + B.

(2.145)Since N1 of Θ and N2 of Q should be known (N1 + N2 = N ), this simultaneous equations canbe solved.
In the diﬀusion problem of the sputtering deposition, source term is given by the sum ofdelta functions, as shown in 2.119. It can be regarded as the sum of ring sources with constantintensity of 1/2πR(x) in the axisymmetric coordinate. Therefore, h(x )/k becomes:h(x )=kThermalized1δ2 (xp − x ).
R(xp )(2.146)

Hence the integrations in Ω̄ in equation 2.143 are simply replaced by the function value at xp ,which gives:(2.147)θ̄(x(i) , xp ).
Bi =ThermalizedCHAPTER 2. THEORY

2.6.6

49

Comparison with MC method

To check the validity of this diﬀusive approach, comparison with the MC method was performed.
Results are shown ﬁgure 2.6. As described above, axisymmetric chamber was assumed in both
calculations, and its boundary was split into nodes as shown in the bottom-left subﬁgure inﬁg. 2.6. The form of this chamber boundary is basically the same with the one actually usedin the experiments (see section 3.1.2), where the nodes from 1 to 9 represent the target surface,and nodes from 31 to 45 represent the substrate holder.
A source of particles was set in this chamber, and both simulations were performed. Several
source positions were tested, and three of them are shown in ﬁgure 2.6. By taking the center of
the target surface as an origin, point A is at (R, Z) = (0, 4), B is at (4, 2.5), and C is at (2, 1) (incm unit). Please note again that we are considering the axisymmetric case. Therefore, that thesource is actually a ring form, and its strength is proportional to R−1 if we see it in 3D system.
MC calculation was performed for 105 of Cu atoms which started at the source position withzero kinetic energy. The argon gas environment was assumed with the temperature of 400 K
and the pressures of 2, 5 and 10 Pa. Born-Mayer potential was used for the scattering potential,
and the collision/scattering events was treated as described in section 2.5. Since the gas motion
is also included in this scattering model, the Cu atoms would start the motion by the collision
impact by an ambient argon gas atom.
Results of the calculation are shown in subﬁgures of A, B and C corresponding to the locationof the particle source. The deposition ﬂux given by BEM calculation is represented in red-broken
line, while that of MC is in marks. The agreement of both is fairly good, but is relatively poor at
nodes where the distance from the particle source is near. It may reﬂect the “constant element”
model of the BEM we used, in which the ﬂux is assumed to be constant on each node, andrepresentative point of each node (xi in eq. 2.137) is located at the center of them. At nodes
near the source, the relative diﬀerence of the distance between the source and the boundary
becomes large, hence the ﬂux should be over- or under-estimated. If we would split the nodeﬁner, or would use “linear element model” or “parabolic element model” for the ﬂux function,
this problem might be overcome.
The edge eﬀect of the BEM [14], which is originated by the singularity of the diﬀusion
equation there, also makes worse the agreement of the BEM result at corner nodes (c.f. nodes 9,10, 37 and 38 in this case). But in most cases where the MC simulation is applied, the interest is
mainly in a uniformity of the proﬁle in the midst of the substrate holder (i.e., away from edges
or corners), hence we probably don’t have to get into this problem so seriously.
It is also observed that the higher the gas pressure the better the agreement becomes. It isbecause of the shorter mean free path in high pressures, which is inﬁnitesimal in BEM case. The
error caused by this eﬀect will be reduced by setting the lower thermalization threshold energy,
where more particles reach the boundary in the MC step. It will be discussed later in the next
section.
Both calculations were executed by a PC with Intel Pentium III (750 MHz). The computation
time is summarized in table 2.2. The reduction of the computation time is obvious. In thesimulation of actual sputter deposition process, you have to sum up the right hand side of theequation 2.146 for every thermalized particle. Therefore, when we obtain the elapsed time of
BEM, θ̄(x(i) , xp ) is re-calculated 105 times (same with the particle number in MC) and summed

CHAPTER 2. THEORY

50

1500

Flux (arb.)

A
1000

2Pa
5Pa
10Pa
BEM

500

0

0

10

20

B

60

30

40
20
0
0

40

Node Number

10

20

30

40

Node Number
300

45

37

A

B
C
1

18

Flux (arb.)

C
200

100

10

0
0
10cm

10

20

30

40

Node Number

Figure 2.6: Comparison between the BEM simulation with MC. Bottom-left shows the crosssection of axisymmetric chamber, where three points of A, B and C denotes the particle source.

CHAPTER 2. THEORY

51

SourceBEM

MC (2 Pa)MC (5 Pa)MC (10 Pa)

A
B
C

18.3
18.6
16.2

93
142
89

401
689
386

1375
2519
1306

Table 2.2: Times consumed to calculate the result of ﬁgure 2.6. The values are in seconds.

up at each nodes.

2.6.7

Thermalization threshold

The BEM procedure described above is used combining with the MC method in the actual
simulation program. At ﬁrst, the particle ejected from the target is treated by MC, where itis gradually decelerated by the collision with gases. When its energy becomes lower than the“thermalization threshold” speciﬁed before the simulation, the MC step for that particle is overand the particle position is recorded as xTi h , which is regarded as one of the particle sources
and used in the right hand side of eq. 2.119. After the MC calculation for all the particles isover, BEM calculation is performed.
The problem arisen here is how to determine this thermalization threshold. When you selecta large value, MC step will be quickly over and the computational time can be reduced. But the
motion of the sputtered particle may not go into the random diﬀusion at such high energy, so
the accuracy of the calculation result may be degraded. Hence we tried several calculations bychanging this threshold, and evaluated the calculation time and the proﬁle of the particle ﬂux
on the substrate holder. For the MC step, procedures described through 2.2 –2.5 are applied(They will be summarized in section 3.1.2 in the next chapter). Copper and argon are selected
as the target material and the sputter gas, respectively. As a sputter chamber, the left bottom
of ﬁg. 2.6 is used.
Figure 2.7 (a) shows the dependence of the computation time to calculate the thickness
proﬁles for 105 of sputtered atoms on the thermalization threshold. Horizontal axis is thethreshold energy in the unit of kT , where T is the gas temperature set to 400 K here. Casesof gas pressures of 2 Pa and 10 Pa are presented. The computer used for this calculation is
IBM PC clone with Intel Pentium III (733 MHz) CPU, the same one used in the previous
section 2.6.6. The computation time was about 5–20 min, and increased as the the thresholdenergy decreased. On comparing the diﬀerent pressures, the 10 Pa case showed steeper increaseat lower threshold energy region. At high threshold energy as high as 1 kT , most of the particles
are thermalized shortly after the sputtering ejection. It only needs a few times of collision, so
the computation time does not diﬀer so much. On the other hand, as the threshold becomes
smaller, more collisions are required for the particle to have —accidentally— smaller energy
than the threshold.
In lower pressures, where the mean path of the particle becomes longer, signiﬁcant part ofthe particles may reach the chamber wall before they thermalized. Therefore, less calculation

CHAPTER 2. THEORY[x103]
1.5

Comp. time (sec)

(a)

52

2Pa
10Pa

1.0

0.5

0

10-1

1

Thremal Threshold (in kT)(b)Particle Flux (a.u.)

80

2 Pa
60
40
10 kT
1 kT
0.1 kT
0.05 kT
0 kT

20
0
0

1

2

3

4

R (cm)
Figure 2.7: Eﬀect of the thermalization threshold on the particle transport simulation combining
MC and BEM. Upper ﬁgure is the calculation time dependence on the threshold, while the lower
ﬁgure gives the resultant ﬂux dependence at substrate.

CHAPTER 2. THEORY

53

time is needed. Moreover, at higher pressures, sputtered particles are more likely to stay near
the erosion track because the diﬀusion coeﬃcient becomes small. So the particles tend to go
back and stick to the erosion track, where they will be re-ejected. It may also result in thecomputation time consumption. When threshold energy was set to 0 K (in other words, BEM
treatment was used and only MC calculation was performed), computation time was 674 ± 3 (s)for 2 Pa and 15355 ± 65 (s) for 10 Pa.
In ﬁgure 2.7 (b), ﬂux of sputtered particles at the boundary nodes on the substrate holder
face is displayed. The substrate holder face corresponds to the node 37 to 45 in ﬁgure 2.6 (a)(the mesh ﬁneness has been modiﬁed, though). The horizontal axis is the distance from thesymmetrical axis. When the threshold energy is greater than 1 kT , the ﬂux becomes smaller.
On the other hand, for the smaller energies, the ﬂux values stay within the calculation error,
and are also identical to the result of zero threshold energy case (i.e., no BEM case). The ﬂux
at the node near the corner edge (R = 3.75 cm) shows a slightly higher ﬂux than its neighbors,which may be caused by the edge eﬀect originated in the BEM calculation as mentioned earlier.
The smaller ﬂux at higher threshold energy (kT ∼ 10kT ) is because of the following reason: thesputtered copper atoms are considered to be thermalized near the target (and treated as theparticle source of random diﬀusion), before they completely lose the initial momentum toward
the substrate. In other words, the centroid of the particle source proﬁle becomes (mistakenly)
far away from the substrates in higher threshold energy cases.
Summarizing, it is preferable to set the thermal threshold to lower value for the accurate
simulation, but more computation time is necessary in this case. In this study, we determined
the threshold at 0.1kT as the compromise of these two requirements.

Chapter 3

Results and Discussions
In this chapter, we perform the simulation based on the theory described in the previous chapter, and show its results. By comparing them with experiments, the validity of models andassumptions is checked out. At the same time, dependence of ﬁlm properties on deposition
conditions (especially gas pressure) are studied with both experiment and simulation, and theorigin of these behaviors is discussed from the viewpoint of the particle transport.
In 3.1, the experimental apparatus used in this study is brieﬂy described. At the same time,
methods and assumptions used in the simulation are summarized. For the latter, the discussion
in the previous chapter is referred from each item.
In 3.2, we simulate the behaviour of high energy copper atoms injected into argon gas. Timeevolutions of copper atoms in both velocity and real spaces are visualized. Especially, the eﬀectof the thermal motion of colliding gas atoms is studied. By introducing the gas motion using themethod in section 2.5, it is shown that the velocity distribution of sputtered atoms approaches
to the Maxwellian with the same temperature of the ambient gas after some period of time. On
the contrary, if the gas motion is neglected, the speed of the sputtered atom is monotonicallydecreased.
In 3.3, sputter deposition process of the copper ﬁlm is simulated. Deposition experiment is
actually performed, and the reproducibility of the thickness proﬁle of the copper ﬁlm by thesimulation is studied. The dependence of the proﬁle on the distance between the chamber walland the substrate is investigated. The pressure dependence of the proﬁle is also studied by meansof the thickness ratio between two substrate holder faces facing toward / away from the target(BFR), which has been proposed by Motohiro [90]. He has found that the relative thickness on
the back face shows a maximum in the middle pressure region by experiment. This behaviour
is reproduced by our simulation model, and the reason of the decrease (and the approach toconstant) in high pressure region is discussed.
In 3.4, we investigate the compositional deviation of LaB6 sputter deposited ﬁlms from thetarget, and its dependence on gas pressure. Its behaviour is resemble to the previous studies
mentioned in section 1.5; as the gas pressure increases from the very low pressure region, theﬁlm content of the lighter element (boron in this case) decreases at ﬁrst, passes a minimum,
and recovers to the stoichiometry at high pressures. It is reproduced both by the experimentand the simulation. Using the simulation, it is shown that this behaviour is originated from the54

CHAPTER 3. RESULTS AND DISCUSSIONS

55

mass diﬀerence between La and B, which aﬀects the scattering process of these atoms with theargon gas atom.
Finally in 3.5, we introduce a novel method to calculate the spatial density of sputteredatoms using the MC simulation. As will be shown in section 3.2, by introducing the thermal
motion of colliding gas atoms, the sputtered atom does not stop and its motion can be traced
until they reach the chamber wall. With this method, the residual time of each sputtered atom
is calculated. Combining this residual time with the deposition rate (i.e. the particle ﬂux on
the wall) obtained by the experiment, atomic density can be deduced. We show the calculation
results for the Cu sputtering plasma, and discussed its correspondence with the optical emission
from copper neutral atoms in the discharge plasma.

3.1

Experimental Apparatus and Simulation policies

3.1.1

Sputtering chamber

Sputter chamber used in this study is shown in ﬁgure 3.1. This is based on a commercial sputter
deposition system (Nichiden Varian: SPF-210B), with some extra features being added. Inner
diameter of the chamber is 20 cm, and the height is 16.5 cm. Target is a disk form whose
diameter is 8 cm. Substrate holder is located above the target, and the distance between themis variable. Because of the shutter motion, the symmetrical axis of the target and the substrateholder displaces from that of the chamber by 3 cm.
A SmCo magnet and a steel yoke are located under the target to form a slightly-unbalanced
magnetron ﬁeld. The ﬁeld strength just above the target is about 0.2 T. With this magnetic
ﬁeld, erosion track was formed on the target as a ring with 2 cm in radius. In this study using Cu
and LaB6 as target materials, the track had about 2 cm width in both cases. The depth proﬁleof them are shown in ﬁgure 3.2. Note that they were formed by the deposition experiments invarious discharge conditions, hence they should be regarded as some kind of average.
Argon was used as a discharge gas, and its pressure during the sputtering was changed from
1 to 20 Pa. The gas pressure was monitored by a Schultz gauge. Flow rate of the argon gas
was kept at 10 sccm with a mass ﬂow controller (Nippon-Tylan SG-3Al-000), and the pressure
was set by throttling the conductance of two evacuation valves connected in parallel to the large
ﬂow-rate type turbo molecular pump (Osaka Vacuum TG200). Both DC and RF discharges
can be generated by switching the DC power source (Advanced Energy MDX-1.5k) and the RF
power source (Tokyo High-Power model A300).
On comparing experimental results with simulations, we should note that the measure ofdischarge gauge is not the pressure itself but the particle density [101]. Between the gauge valueof the pressure PG and the gas density nG , the following relation can be applied:
nG =

PG
,
kTG

(3.1)

where TG is the speciﬁc temperature given by the maker of the Schultz gauge, which is 296 K
in the model we used.
Though the gas density is more essential than the pressure for the particle transport, it israther unfamiliar in experimental ﬁelds. In the following discussion, therefore, we use the eq. 3.1

CHAPTER 3. RESULTS AND DISCUSSIONS

56

x

20 cm
Figure 3.1: Sputter chamber for ﬁlm deposition experiments.

CHAPTER 3. RESULTS AND DISCUSSIONS

57

0

Depth (mm)

0.1
0.2

Cu

0.3
0
0.1
0.2
0.3
0

LaB6
1

2

3

4

Radial Position (cm)Figure 3.2: Radial depth proﬁle of erosion tracks formed on Cu and LaB6 targets.

CHAPTER 3. RESULTS AND DISCUSSIONS

58

inversely, and display the value in pressure unit. So if the pressure value appears in the text inthis chapter, please regard it as the density implicitly.

3.1.2

Summary of simulation policies

In this section, models, methods and assumptions of the simulation used in this study aresummarized. Most of them have already been described in the previous chapter. Citation tothem are added on each item if possible.
1. The sputter chamber (i.e., boundary) is assumed to be axisymmetric to simplify the simulation. Note that in the experimental apparatus, symmetric axis of the target–substrate
system is 3 cm apart from that of chamber.
2. All the sputtered particles are atoms. It is partly because that the collision and scattering
between molecules and clusters are diﬃcult currently. But for the cases covered in thisstudy (Cu and LaB6 ), we believe this assumption is adequate.
3. The ejection position of sputtered atoms on the target is determined so that the depth
proﬁle of the erosion track is reproduced (section 2.2.1). These depth proﬁles used in thisstudy are shown in ﬁgure 3.2, which have been actually measured on each target. If the
sputtered atom comes back to the erosion track and sticks there during the Monte Carlo
step, it will be re-ejected from there immediately.
4. On ejection from the target, the angle and energy distributions of sputtered atoms are
ruled by the cosine distribution and the Thompson’s formula, respectively (sections 2.2.2
and 2.2.3).
5. The pressure and the temperature of the ambient gas are uniform in the chamber.
6. The mean free path (MFP) of the sputtered atom is determined using the method describedin section 2.5.1, and the actual free path on each event is deduced by the method in
section 2.3.1 with this MFP value.
7. The speed and the direction of the colliding gas atom to the sputtered atom are determinedusing the method in section 2.5.2.
8. The scattering between the sputtered atom and the gas atom are calculated based on thepotential scattering framework. The Born-Mayer potential with the parameters given by
Abrahamson is used as the scattering potential (section 2.3.3).
9. The speed and the direction of the sputtered atom after the scattering is determined as
described in section 2.5.3.
10. When the sputtered atom arrives on the chamber boundary, it sticks there with the probability of unity (section 2.4). However, if the boundary belongs to the erosion track region,
the atom will be re-ejected as mentioned above.

CHAPTER 3. RESULTS AND DISCUSSIONS

59

for LaB6

for Cu10 cm
Figure 3.3: Axisymmetrical chamber boundary used with MC simulation. Location (and shape)of the substrate holder diﬀers for LaB6 (3.4) and for Cu (3.3 and 3.5) depositions. Radius of
the chamber wall is changed to study the eﬀect of it in 3.3.

CHAPTER 3. RESULTS AND DISCUSSIONS

60

11. When the kinetic energy of the sputtered atom becomes less than the speciﬁed threshold,it is regarded as “thermalized” there, and treated as the particle source with the strengthof unity. These sources are summed up, and the diﬀusion equation corresponding to thissource term is solved after the MC calculation is over. The particle ﬂux deduced by thisprocedure is added to the ﬂux of the atoms arrived to the wall before the thermalization(section 2.6.2).
12. The thermalization threshold energy is assumed to be 0.1 kT , where k is the Boltzmann’sconstant and T is the gas temperature (section 2.6.7).
13. The boundary condition of the diﬀusion equation on normal boundary is the Dirichletcondition of zero density, and the one on the erosion track is the Neumann condition withzero ﬂux (section 2.6.2).
14. The diﬀusion equation is solved by the Boundary Element method (BEM) as shown insection 2.6.5.

3.2

Deceleration and diﬀusion of high energy atoms in gases

In this section, we simulate the behaviour of high energy atoms injected into gas environment;
deceleration by the collision with gas atoms and diﬀusion. Especially, we discuss the diﬀerence
between two cases: whether the thermal motion of gas atoms is included or not. The situation
with copper atoms and argon gas was studied.
On simulating the time evolution of the velocity distribution, 5000 copper atoms were injected
into the gas at 0 s, to the positive direction of the z axis with 5 eV of kinetic energy. Then,the snapshots were taken after some periods of time. On calculating the elapsed time, we
assumed that the copper atom changed the speed just instantly at collision points, and ﬂew at
the constant speed between two collisions (in other words, integration of the time in potential
scattering framework was not made). The gas density is set to 10 Pa in the unit of vacuum
gauge indication (see eq. 3.1).
Figure 3.4 shows the result. Snapshots taken at 2, 10 and 50 μs are plotted in the velocity
space. Figures (a), (c) and (e) are the results considering the gas motion with temperature
Tg = 400 K, while (b), (d) and (f) are the ones calculated with Tg = 0, i.e., neglecting thegas motion. The range of each axis (Vx , Vy and Vz ) in ﬁgure 3.4 is taken from −3000 to 3000
m/s. The speed of the copper atom at kinetic energy of 5 eV is about 3900 m/s. And if
copperatoms are ruled under the Maxwellian with the temperature 400 K, their mean speed is
8kT /πm = 365 m/s.
Cu atoms are put at Vz = 3900 m/s at 0 s, and their directions are quickly scattered. If you
estimate the eﬀective collision radius between Cu and Ar atoms to be 3 Å (see ﬁgure 2.2), themean free path of copper in 10 Pa (gauge value) of Ar gas is deduced to be about 1.4 mm from
eqs. 2.26 and 3.1. Therefore, in 2 μs, copper atom of 3900 m/s speed will collide 5 times withargon gas atoms. Actually, since the speed of the copper atom decreases by each collision, thismay be reduced to 3 times or so.
After that, the velocity distribution converges toward the origin of the velocity space as thetime goes by. The cases with Tg = 400 K and with 0 K look similar at 10 μs, while these are

CHAPTER 3. RESULTS AND DISCUSSIONS

61

3000

V Z (m/s)

(m/s)

3000

V

Z

0

-3000
3000

3000

0

-3000
3000

-3000

V Z (m/s)(m/s)

V X (m/s)

3000

V

Z

0

-3000
3000

3000

0

0

-3000
3000

-3000

V

X

3000

0

0
-3000

(m/s)

V

(m/s)

Y

0
-3000

(m/s)

(c) T=400K 10 μs

-3000

V

X

(m/s)

(d) T=0K 10 μs

3000

(m/s)

3000

(m/s)

-3000

(b) T=0K 2 μs

3000

Y

0
-3000

V Y (m/s)

V X (m/s)(a) T=400K 2 μs

V

3000

0

0
-3000

V Y (m/s)

0

Z

0

V

V

Z

0

-3000
3000

3000

0

V Y (m/s)

-3000
3000

-3000

V X (m/s)(e) T=400K 50 μs

3000

0

0
-3000

V Y (m/s)

0
-3000

-3000

V X (m/s)(f) T=0K 50 μs

Figure 3.4: Time evolution of the velocity space distribution of 5000 Cu atoms injected with
5 eV into Ar at 10 Pa. Left subﬁgures ((a), (c) and (e)) are the cases where the gas temperature
is 400 K, while the right ones ((b), (d) and (f)) are the cases with the gas temperature is 0 K(e.g., gas motion is neglected).

CHAPTER 3. RESULTS AND DISCUSSIONS

62

clearly diﬀerent at 50 μs. In the former case (Tg = 400 K), the velocity distribution keeps somesize as large as the thermal speed (385 m/s), while in the latter case (Tg = 0 K), the distribution
continuously converges into a very small volume. With ﬁnite gas temperature, the velocity ofcopper should approach to the equilibrium with the gas. On the other hand, if the gases areimmobile, copper atoms always decelerated on each collision event and their speed approaches
to zero.
Figure 3.5 shows the elapsed time dependences of the speed distribution of these copper
atoms (“speed” means the absolute value of the velocity). Figure 3.5 (a) is the case where theargon gas temperature is 400 K, corresponding to the ﬁgure 3.4 (a), (c) and (e). On the contrary,ﬁgure 3.5 (b) is the zero gas temperature case, corresponding to the ﬁgure 3.4 (b), (d) and (f).
Times when the snapshot is taken are in the range of 2–50 μs, which are also the same withﬁgure 3.4. On calculating the proﬁles shown in ﬁgure 3.5, however, 105 of particles are traced.
The function drawn in a solid line in ﬁg. 3.5 (a) is the Maxwellian appeared in eq. 2.70:
4 3f (vg ) = √ α 2 vg 2 exp (−αvg 2 ),
π
m
α=
2kT

(2.70 )with m to be the atomic mass of copper (63.54 × 1.661 × 1027 kg) and with T to be the same
temperature with gas (400 K). You can see that the speed distribution of copper atoms is
reproduced very well by the simulation considering the thermal motion of gases.
The velocity distribution f0 (vp ) in Tg = 0 K environment seems to obey the scaling law:
f0 (γvp )= constγ(3.2)

with the scaling parameter γ(t), which is found to have the time dependenceγ ∝ t1.2∼1.3 .

(3.3)

What determines this time exponent is left unclear at this time of writing. If argon atom isstationary, the rate of energy decrease is deduced from the scattering angle. As the copper
speed (and kinetic energy) decreases, tail of the Born-Mayer potential becomes important, andthe same scattering angle is obtained with larger collision parameter. This reduces the mean
free path of the copper atom, and enlarge the time exponent. On the contrary, the time between
the collision becomes larger as the speed decreases. This diminish the exponent. These two
eﬀects may compete each other and determine the time exponent.
The propagation process of copper atoms in real space has also been investigated, and are
displayed in ﬁgure 3.6. The conditions are the same with ﬁg. 3.4, 5000 of 5 eV copper atoms are
injected at 0 s, along the positive direction of the z axis from the origin. Argon gas temperature
is 400 K and the pressure is 10 Pa. x, y and z axes in ﬁg. 3.6 are the coordinates in real space.
Their range is from −5 to 5 cm, which is comparable with the size of the conventional sputter
chamber. The allocation of subﬁgures is also identical with ﬁg. 3.4; left sequence (a), (c) and(e) are the cases where the argon gas temperature is 400 K, while the right (b), (d) and (f) are

CHAPTER 3. RESULTS AND DISCUSSIONS(a)dN/dv [s/m]10-2 CuAr(400K)

10-3

63

2 μs
5 μs
10 μs
20 μs
50 μs
Maxwellian

10-4
10-5
0

1000

2000

Speed [m/s](b)
dN/dv [s/m]10-2

CuAr(0K)

2 μs
5 μs
10 μs
20 μs
50 μs

10-3
10-4
10-5
0

1000

2000

Speed [m/s]
Figure 3.5: Time evolution of the energy distribution of Cu atoms injected at 0 s with 5eV into
Ar gas of 10 Pa.

CHAPTER 3. RESULTS AND DISCUSSIONS

64

5

Z (cm)

Z (cm)

5

0

-5
5

5

0

-5
5

-5

(a) Tg =400 K 10 μs

Z (cm)

Z (cm)X (cm)

5

0

-5
5

5

0

0

-5
5

-5

(c) Tg =400 K 50 μs

0-5

Y (cm)X (cm)

500-5

Y (cm)-5

X (cm)(d) Tg =0 K 50 μs

5

Z (cm)

5

Z (cm)-5(b) Tg =0 K 10 μs

5

0

-5
5

5

0

Y (cm)

0-5

Y (cm)X (cm)

500-5

Y (cm)

0

0

-5
5

-5

-5

X (cm)(e) Tg =400 K 200 μs

5

0

0

Y (cm)

0-5-5

X (cm)(f) Tg =0 K 200 μs

Figure 3.6: Time evolution of the real space distribution of 5000 Cu atoms injected into Ar at
10 Pa.

CHAPTER 3. RESULTS AND DISCUSSIONS

65

the one with 0 K. The times to take the snapshots in this ﬁgure are 10, 50 and 200 μs. Notethat they are diﬀerent from those in ﬁg. 3.4 to make clear the diﬀerence between two cases.
You can see in Tg = 0 K case that the position proﬁle at 200 μs (f) is not so diﬀerent fromthe one at 50 μs (d). On the other hand, in Tg = 400 K case, diﬀusion of copper particles
proceeds until 200 μs. As shown in ﬁg. 3.5, the diﬀerence between these two temperatures
becomes signiﬁcant since 50 μs or so. The velocity distribution approaches to the Maxwellian
in 400 K case, while it converges to vp = 0 by neglecting the gas motion. In the latter case,therefore, copper atoms are frozen at some point, and do not move anymore from there.
From the standpoint of the particle transport simulation, this “freeze” of copper atoms may
not be the major problem if only the atomic ﬂux on chamber wall is necessary, since the diﬀusion
stage can be treated with the framework described in 2.6 using the Poisson’s equation. In thisframework, copper atoms, which have high energy initially, are decelerated during the ballistictransport process, and turn into the random diﬀusion. If the positional distribution, whichstarts from this frozen proﬁle and evolutes under the diﬀusion equation in some period, would
be similar to the one obtained in the ﬁnite temperature case, the ﬂux proﬁle in both cases shouldlook alike each other.
It should be noted, however, that the gas motion must be incorporated when you have toconsider the motion of sputtered atoms including the time parameter (as will be discussed in
section 3.5). If you neglect it, sputtered atoms will stop inside of the chamber, hence it iscompletely impossible to deduce the time corresponding to the transport and/or diﬀusion.

3.3

Thickness proﬁle of sputtered Cu ﬁlms

In this section, we simulate the actual sputter deposition process, and compare the thickness
proﬁle with experiments. As in the case for general sputtering apparatus, the sputter target andthe substrate holder is placed facing to each other in this study. It is shown that the thickness
proﬁle is aﬀected by the distance between this target–holder system and the chamber wall. The
reason of this eﬀect is discussed from the standpoint of the particle transport.
During the sputter deposition at pressures not so low, ﬁlms also deposit on chamber walls
which cannot be seen directly from the target, by the diversion of sputtered atoms deﬂected bycollisions with gas atoms. It has been reported by Motohiro that the thickness ratio on the backand front faces of the substrate holder (BFR) shows a interesting pressure dependence [90]. He
observed that this BFR took a maximum at midst pressure and decreased in higher pressures.
But in his study, the reason was not fully understood because of the limitation of simulation
model and the computer performance. In this section, we apply our novel methods to thisproblem, and discuss the reason of it.

3.3.1

Experiment and simulation

The deposition experiment was performed with the apparatus described in section 3.1.1. Target
was a copper metal disk with 4N purity. Argon was used as the sputter gas, and the discharge
was generated with the DC power source. Corning 7059 glass substrates were used, whichwere set on both faces (front/back) of the substrate holder shown in ﬁg. 3.1. Film thickness
was measured by the multiple interferometer. Step was prepared by the cover glass mask with

CHAPTER 3. RESULTS AND DISCUSSIONS

66

0.12 mm thickness. To evaluate the eﬀect of the neighbor wall, thickness proﬁle was taken along
the line passing both the symmetric axis of the substrate holder and that of the chamber. The
former was taken to be the origin, and the positive direction of the measurement axis (x axis
in ﬁgure 3.1) was chosen heading toward the nearer wall. Since these centers are eccentric, thethickness proﬁle along the x axis is also asymmetric by x = 0.
The simulation was carried out in a similar manner with the previous section 3.2, where wehave seen the behavior of copper atoms injected into the free space ﬁlled with argon gas. In addition to it, we used the methods described in section 3.1.2 to determine the position/energy/angleof the sputtered copper atoms from the target. The BEM treatment of the diﬀusive transport
was also incorporated. We also considered the chamber boundary as shown in ﬁgure 3.3. To
investigate the eﬀect of the distance to the wall, three cases with chamber radii of 7, 10 and
13 cm were tried. These correspond to the nearest distance to the wall, actual radius of thechamber, and the furthest distance to the wall in the experimental apparatus, respectively.
Figure 3.7 is the snapshots of the position of copper atoms taken during this simulation. In
this ﬁgure, argon gas pressure is assumed to be 10 Pa 1 , and the gas temperature is 400 K. The
outer radius of the boundary is 10 cm (not drawn). 5000 copper atoms are plotted in the ﬁgure;red points denote the atoms still ﬂying in the chamber, while the blue ones denote the atoms
already deposited on the wall. As you can see in this ﬁgure, most of the atoms are still ﬂying
even 200 μs has elapsed. At this time, all copper atoms are thermalized and its velocity proﬁleobeys the Maxwellian. It means that the thermal motion of gases is important to simulate theparticle transport in high pressure sputtering.

3.3.2

Thickness proﬁleFigure 3.8 shows the experimental and numerical results of the thickness proﬁle of Cu ﬁlms, forgas pressures of 2 Pa and 10 Pa. Horizontal axis is the position along the x axis, as described
previously. Nearer wall of the chamber is located at the positive side x = 7 cm, while the further
wall is at x = −13 cm. The proﬁle is normalized by the largest thickness, which was not always
at x = 0. The results on both front and back faces are plotted simultaneously, the much smaller
data are those at back face.
The marks in the ﬁgure are the experimental results, while the lines are simulation results.
Since our program code can only account for the axisymmetric cases, we performed the simulation for two cases, with chamber radius of 7 cm and 13 cm, which were the nearest and thefurthest distances between the target center and the chamber wall, respectively. The result of
13 cm is plotted in a solid line in the region of x < 0, while that of 7 cm is plotted in a broken
line in x > 0.
It is observed on the front face that the ﬁlm thickness takes the maximum at above the target
center, and becomes thinner as going away from there. By comparing the both halves, nearer
and further from the chamber wall, the decrease of the thickness is steeper on the nearer side(x > 0). It can be explained by considering the diﬀusion process. If the sticking coeﬃcient of
the particle is unity on the chamber wall, it can be regarded as the boundary with zero particle
density, as discussed in 2.6.1. Hence, nearer the chamber wall, smaller the density of particles. It
1

Again, this is the vacuum gauge value and indicates the gas density (see eq. 3.1).

CHAPTER 3. RESULTS AND DISCUSSIONS

67(a) 5 μs(b) 10 μs(c) 20 μs(d) 50μs(e) 100 μs(f) 200μs

Figure 3.7: Time evolution of the real space distribution of 5000 Cu atoms: Actual deposition
condition is assumed where the argon gas pressure is 10 Pa and primary ion energy is 300 eV.

CHAPTER 3. RESULTS AND DISCUSSIONS

68

1.0

Normalized Flux (a.u.)

0.5

Pressure: 2Pa
Calc (r=13cm)
Calc (r= 7cm)
Experiment

0
-4

-2

0

2

4

1.0

0.5

Pressure: 10Pa
Calc (r=13cm)
Calc (r= 7cm)
Experiment

0
-4

-2

0

2

4

Substrate Position (cm)Figure 3.8: Thickness proﬁle of sputter-deposited Cu ﬁlms.

CHAPTER 3. RESULTS AND DISCUSSIONS

69

is reproduced by the simulation also. You can see that the dotted lines show the rapid decreaseas they come far from the center on both 2 Pa and 10 Pa cases.
While the agreement between the experiment and the simulation was fairly wall on the nearerside (x > 0), it was rather poor on the further side (x < 0). The experimental result shows
more moderate decrease there. It may be originated from the density gradient reﬂecting thesurrounding wall conﬁguration or the ﬂow ﬁeld of the argon gas, but we have no clear explanation
for this yet anyway. If full-3D simulation would become possible in the future, it should be worth
to evaluate this gradient with it.

3.3.3

Diversion of sputtered atoms onto the back face

As mentioned in the ﬁrst part of this section, sputtered atoms are deﬂected by the collision
with gases, hence they may arrive at the chamber wall which cannot be seen directly from thetarget. Now let us consider the thickness proﬁle on the back face of the substrate holder, as
schematically shown in ﬁgure 3.9. On the contrary to the front face, the proﬁle on the back facebecomes thinner as entering inside from the edge of the substrate holder.
In the low pressure limit, where no collision occurs, atoms cannot arrive at the back face.
So you may think that the higher the pressure, the thicker the backface ﬁlm is. But it is nottrue. On extremely high pressure, the relative thickness on the back face becomes thinner.
Actually, the result of nearer half (x > 0) at 10 Pa lacks the experimental data on the back face,because they were too thin to be measured with the interferometer. In the following, we deﬁnethe BFR (back-front ratio) as the thickness ratio of back/front faces at same x. We obtain the
BFR on further wall side (x < 0) experimentally, and compare its pressure dependence with thesimulation.
The pressure dependence of the BFR is shown in ﬁgure 3.9. In this ﬁgure, the experimental
and simulation data are plotted with ﬁlled and open marks, respectively. It was taken on thefurther half (x < 0), 30 mm and on 35 mm from the center of the substrate holder (holder radius
is 40 mm). The position x = −35 mm is nearer at the holder edge, hence the ratio is larger
than at x = −30 mm. Experimental scatter was rather large because of the thin ﬁlm thickness
at back face, but the agreement of experiment and simulation on both location is satisfactory.
And it also should be noted that the tendency of BFR in this ﬁgure agrees well with the results
reported by Motohiro [90].
This BFR behavior can be explained considering the transition of the particle transport,
from the ballistic one to the diﬀusive one. As mentioned above, at the low pressure limit,
sputtered atoms cannot deposit on the back face. As the pressure increases, the particles cometo be able to arrive there by the scattering, and the arrival rate increases gradually. At higher
pressures, the sputtered atoms tend to be thermalized before they deposit on the wall. They
are thrown away some distance from the target through ballistic transport, and are thermalized.
From there, they start diﬀusion, which is actually the random walk in the MFP scale. In thisregion, the thermalization proﬁle approaches nearer to the target as the gas pressure increases.
It means that the diﬀusion source of the particles come closer to the substrate holder 2 . It resultsin the decrease of the atomic ﬂux at the back face, since the distance from the source becomes2

see also section 2.6.

CHAPTER 3. RESULTS AND DISCUSSIONS

70

40mm
R/F Thickness Ratio

0.10

-35mm

0.05

-30mm
0

1

5

10

Pressure (Pa)
Figure 3.9: Ratio of deposition rates at front/back face of the substrate holder. Filled and open
marks denote the experimental and simulation results,respectively.

CHAPTER 3. RESULTS AND DISCUSSIONS

71

relatively larger. At the high pressure limit, where all atoms are thermalized just after theejection from the target, erosion track itself can be regarded as the particle source. Therefore,
the BFR becomes constant in this region.
In ﬁgure 3.10, we show the gas pressure dependence of the thermalization proﬁle obtainedby the simulation for 106 copper atoms. The density of the thermalization atom is displayedby the color map in linear scale of arbitrary unit. Thermalization threshold was set to 0.1kTg ,as described in section 3.1.2. Sputtered atoms may start the random diﬀusion at more higher
energy, but ﬁg. 3.10 should still give the information we need. As can be seen from the ﬁgure, thethermalization comes to signiﬁcant at around 1.5–2 Pa, where the BFR takes the maximum in thesimulation. Above this pressure, it is clearly shown that the thermalization proﬁle approaches
toward the target. Thus the BFR decreases in high pressures, and becomes constant.

3.3.4

Summary

In this section, we have simulated the sputter deposition process of copper ﬁlms. The thickness
proﬁle on the substrate holder was changed by the distance from the chamber wall. When thewall was near, the proﬁles was thinner. It is also conﬁrmed by the simulation, and can be
explained by considering the chamber wall as the zero density boundary.
We then deﬁned the BFR ratio (thickness ratio on back and front faces of the substrateholder), and investigated its pressure dependence. The simulation reproduces its behavior;
BFR increases with increasing the gas pressure in low pressure region, and then decreases andapproaches to constant. It also coincides with the experimental result reported by Motohiro.
The pressure dependence of the BFR can be explained as follows. In lower pressures, where
the ballistic transport is dominant, the increase of the gas pressure causes more collisions andraise the BFR. On the other hand, at high pressures where the diﬀusive transport occurs, thermalization proﬁle concentrate into the neighbors of the erosion track as the pressure increases.
Thus the BFR decreases and approaches to the constant value, where the erosion track itself
can be regarded as the particle source of the diﬀusion.

3.4

Composition of LaB6 sputtered ﬁlms

In this section, we study the compositional deviation of sputtered LaB6 ﬁlms and its gas pressure
dependence. As mentioned in section 1.5, sputter deposition technique has been applied widely
not only to elemental materials but also to compounds and alloys. In the sputter deposition ofthese multi component materials, it has been frequently observed that the composition of the
ﬁlm deviates from that of the target, which critically aﬀects the physical properties of the ﬁlm.
In these investigations it has also been pointed out that the extent of the deviation depends onthe gas pressure.
The reason for the compositional deviation may be classiﬁed into three categories corresponding to the stages of the sputter deposition; (1) processes on the target, (2) particle transport
process through the ambient gas, and (3) processes on the substrate. As for the ﬁrst stage(i.e. when the target material is sputtered), there may be the diﬀerence in angular- and energydistributions and/or sputtering yield between composing elements. However, these should notdepend strongly on the gas pressure. It has also been pointed out that the sputtering yield

CHAPTER 3. RESULTS AND DISCUSSIONS

Figure 3.10: Distribution of thermalized positions of sputtered particles.

72

CHAPTER 3. RESULTS AND DISCUSSIONS

73

diﬀerence should be compensated by the resultant deviation of surface composition [169]. As
for the second, the diﬀerence in scattering behavior of sputtered particles has been pointedout [123, 173], which is also the main subject in this section. As for the third, the diﬀerence
in sticking coeﬃcient or re-sputtering rate may cause the compositional shift. For example,
the eﬀects of negatively ionized oxygen molecules have been intensively studied for oxide ﬁlm
deposition [32,127]. But in case of non-reactive sputtering, where the generation of negative ion
hardly occurs, this eﬀect can be neglected.
In the present study, we concentrate on the second stage. In a standard sputter-deposition
system, inert gas of 0.1–10 Pa is introduced to generate the plasma discharge. In this case, sputtered atoms experience several times of collision by ambient gas atoms. As shown in section 3.2,
they gradually lose their energy and change directions with these collisions, and their motion
become more diﬀusive rather than ballistic. Flux of sputtered atoms arrived on substrate should
be strongly aﬀected by the number of collision events, which is dependent on the gas pressure.
In this section, we have sputtered LaB6 and studied the compositional dependence of deposited ﬁlms on the gas pressure by both experiments and MC simulation. LaB6 is the material
which has preferable features for hot cathode; e.g., high melting point, low vapor pressure, andlow work function. Preparation of LaB6 ﬁlms by sputtering technique has already been reported
by several authors [94, 126], in which the compositional shift problem is also pointed out.
To study the compositional deviation problem, LaB6 happens to be suitable because of thefollowing reasons. Diﬀerence in the mass between lanthanum and boron is extremely large, sothat the diﬀerence in the scattering behavior with argon gas atom is emphasized. In addition,
La–B system is in solid phase independent of the composition as high as 2000 K [57]. Therefore
the desorption of deposited atoms can be neglected. And ﬁnally, it has already been conﬁrmedthat only single atoms (no clusters) were observed with mass spectroscopy measurement in
sputtering of LaB6 [59]. Hence we can apply a simple simulation model.

3.4.1

Experiments

Lanthanum boride (LaBx ) ﬁlms were deposited using the sputtering apparatus described in
section 3.1.1. Target was a sintered LaB6 disk of 99.9 % purity, and the target-substrate distance
was set to 7cm. For LaB6 deposition, RF power source was used. RF power applied to thetarget was 100 W. Argon was used as a discharge gas, and its pressure during the sputtering
was changed from 2 to 15 Pa. For conditioning the target surface, pre-sputtering was carriedout for an hour before opening the shutter and starting deposition.
Films of 200–300 nm thickness were deposited on the 22 × 15 mm2 of (boron free) Soda-Limeglass substrate. The composition of the ﬁlms were determined by ICP-Mass (Inductively coupledplasma with mass spectrometer) method where the ﬁlms were dissolved by nitric acid, and theconcentration of La and B in the solutions were measured. By dividing them with the ﬁlm area
and the deposition time, atomic ﬂuxes of two elements were deduced.
Figure 3.11 is the result of the ICP-Mass measurement. A process pressure dependence of
the atomic ﬂuxes is shown in (a), and the B/La ratio of the ﬂux is shown in (b). The atomic
ﬂux dependence on the gas pressure can be separated into two regimes at around 5 Pa. In low
pressure region, atomic ﬂux of both elements gradually decreased with increasing the pressure.
Decrease of boron was slightly faster than La, hence the composition of the ﬁlms gradually

CHAPTER 3. RESULTS AND DISCUSSIONS

74

B
1016

La
1015

(b)

7

B/La at.ratio

Atomic fluxl (cm-2min-1)(a) 1017

6
5
4
3

5

10

Pressure (Pa)
Figure 3.11: Result of ICP-Mass measurements of LaB6 sputtered ﬁlms: (a) atomic ﬂux (b)B/La ratio.

CHAPTER 3. RESULTS AND DISCUSSIONS

75

shifted to La-rich. In higher pressures, on the other hand, the ﬂux sharply dropped, but theratio of the two elements recovered to the target composition. Deposition rate of the ﬁlms
was also measured by the multiple interferometry. It was decreased constantly with pressure
increase; 7.3 nm/min at 0.5 Pa and 0.5 nm/min at 15 Pa. These values agreed with the one
calculated from the ﬂux in ﬁg. 3.11 (a) and atomic density of LaB6 .
At the process pressure less than 3 Pa, deposited ﬁlms had purple-gray color, as reported
previously [94] for stoichiometric LaB6 ﬁlm. In these pressures, target also exhibited fresh purple
surface as before the sputtering. On the contrary, at higher pressures, the ﬁlms looked like inmetallic gray, and target were covered with griseous-colored overlayer except for the erosion
tracks, on which the target surface still looked like that at lower pressures.

3.4.2

Simulation

The Monte Carlo simulation has been carried out as described in section 3.1, except for thefollowing point. In this simulation, we supposed that all sputtered particles had the identicalenergy when they ejected from the target. It is because, as far as we know, there is no establishedtheory for the energy distribution of the sputtered particles from compounds. We also assumed
that the ejection energy was identical for B and La. In the sputtering theory, energy transfer isgoverned by the value Γsg deﬁned by:Γsg =4Ms Mg(Ms + Mg )2(3.4)

where Ms and Mg are mass of colliding two atoms. Here the Γsg for La with Ar (0.69) and thatfor B with Ar (0.67) are nearly equal. That’s why we assumed the same ejection energy for both
B and La. The simulation was carried out with ejection energies of 5 and 10 eV. As describedabove, all of sputtered particles were assumed to be atoms.
The chamber boundary shown in ﬁgure 3.3 was used. The distance between the target andthe substrate was set to 7 cm in accordance with the experiment. Reﬂecting the form of substrateholder used in the experiment, the thickness of the holder and the rod hanging the holder werethinner in this simulation. One million (106 ) atoms were traced for both lanthanum and boron,
and the numbers of them arrived at the substrate position were counted. Actually the 10 trials
with 105 atoms were performed, and the calculation error was evaluated. The resultant B/La
ratio was multiplied by 6 and plotted to compare with the experimental results.
Figure 3.12 shows the composition of sputtered LaB6 ﬁlms obtained by the simulation.
For comparison, experimental result is also shown in (a) again, while (b) is the simulation
result. B/La values are plotted for diﬀerent ejection energies (5 and 10 eV) and on diﬀerent
positions away from the center of the substrate holder (nodes of 0–5 mm, 10–15 mm, and
20–25 mm). Reﬂecting the large distance between the target and the substrate, signiﬁcantdependence on the substrate position is not observed. Simulation does reproduce the dependenceof ﬁlm composition on gas pressure; B/La is stoichiometric at low pressure, decreases withincrease in gas pressure, passes a minimum, and recovers to 6 again in high pressure region. The
minimum of the B/La is about 4 both in experiment and in simulation. However, the pressuresat which the B/La becomes minimum are diﬀerent. In experiment, it appears at around 10 Pa,

CHAPTER 3. RESULTS AND DISCUSSIONS767(a) Experiment6

B/La at.ratio

5
4
3
7

(b) Simulation
65
Pos 5eV 10eV
0- 5
10-15
20-25

4
3
0.5

1

5

10

Pressure (Pa)
Figure 3.12: Composition of sputtered LaB6 ﬁlms: (a) experimental result, (b) simulation result
for ejection energies of 5 and 10 eV’s. For (b), composition at diﬀerent radii are plotted (seetext).

CHAPTER 3. RESULTS AND DISCUSSIONS

77

while it does at 2–3 Pa in the simulation, and slightly shifts to the higher pressure with highejection energy (10 eV).

3.4.3

Mass diﬀerence eﬀect on energetic particle transport

In this subsection, we discuss the origin of the “U-shaped” dependence of B/La ratio on gas
pressure, in view of the particle transport process of B and La atoms.
Neglecting the motion of argon gas atom, boron atom may be scattered backward by ambient
argon atoms, while no backscatter occurs for lanthanum. This is shown in ﬁgure 3.13, which isthe scattering angle of B and La atoms in laboratory system calculated using the formulationin section 2.3.2. This applies for atoms in the initial ballistic transport stage, since sputteredatoms from the target initially have energy as high as several electron volts, hence much fasterthan Ar atoms.
The diﬀerence can be understood as follows: When the sputtered atom collide with the gas
atom, its velocity after the collision is given by the eq. 2.45:vs =ms v s + mg v R 
ms + mg(2.45 )where ms and mg are the mass of sputtered atom and gas atom, respectively, and  denotes thevalues after the scattering. If the sputtered atom is much faster than the gas atom, the absolutevalue of relative velocity v R (= v g − v s ) is almost identical with v s . Furthermore, in the elastic
scattering, the absolute value of v R does not change, but only the direction changes. Hence, theabsolute values of v s and v R  in eq. 2.45 are equal. If the sputtered atom is heavier than gas
atom (ms > mg ), initial velocity v s aﬀects vs  more strongly than v R  , and vise versa.
In the center of mass system, the angle between v R and vR  is 2π (i.e. opposite direction)
when the collision parameter b is zero. And as b becomes larger, the angle decreases toward
zero at b = ∞ (see ﬁgure 2.1). Therefore, for the lighter atom, backscattering can occur in thelaboratory system at small b. For the heavier atom, on the other hand, it does not occur andonly forward scattering is allowed.
As carried out in section 3.2, time evolution of the velocity distribution of B and La atoms are
displayed in ﬁgure 3.14. (a), (c) and (e) are for boron, while (b), (d) and (f) are for lanthanum.
The atoms are injected into the argon environment at 0 s along the positive direction of z axis,
with the initial kinetic energy of 5 eV. The temperature of argon gas is assumed to be 400 K,
and the eﬀect of thermal motion of gas atoms is taking into account. The gas pressure is set
to 3 Pa, where the B/La ratio shows the minimum in the simulation (see 3.12). Snapshots aretaken at 10, 20 and 50 μs, which is slightly later as compared to ﬁg. 3.4, reﬂecting the lower
pressure of the gas (and larger MFP). Since their masses are quite diﬀerent, the scale of thevelocity axes also diﬀers for B (−10000 ∼ 10000 m/s) and for La (−3000 ∼ 3000 m/s).
The diﬀerence between the elements reﬂects the diﬀerence in scattering angle with argon
gases shown in ﬁg. 3.13. In boron case, their direction is quickly scattered and the velocity proﬁlelooks centrosymmetric already at 10 μs. After that, the proﬁle converges into the Maxwellian
keeping the symmetry. On the other hand, in La case, the proﬁle preserves the initial momentum
to the z direction, and seems to “drift” into the Maxwellian as a whole. The energy distributionof boron atoms converges to Maxwellian earlier than that of lanthanum, because the speed of

CHAPTER 3. RESULTS AND DISCUSSIONS

78

3

B
0.1eV

Scattering Angle (rad)21eV

1 10eV
0
3

La

2

1
10eV 1eV 0.1eV

0
0

1

2

3

4

5

Collision Parameter(Å)
Figure 3.13: Scattering angle of B and La atoms when they collide with an stational Ar atom.

CHAPTER 3. RESULTS AND DISCUSSIONS

79

4

x 10
1

V Z (m/s)

V Z (m/s)

3000

0

-1
1
4

x 10

V

Y

0
-1

(m/s)

-3000
3000

1

0
-1

0

0
-3000

V Y (m/s)

V X (m/s)

3000

0

4

x 10

-3000

V X (m/s)(a) B 10 μs(b) La 10 μs

4

x 10
1

V Z (m/s)

V Z (m/s)

3000

0

-1
1
x 10

0
-1

V Y (m/s)

-3000
3000

1

0

4

-1

0

V

V X (m/s)

3000

0

4

x 10
Y

0
-3000

(m/s)

-3000

V

(c) B 20 μs

X

(m/s)

(d) La 20 μs

4

x 10
1

V Z (m/s)

V Z (m/s)

3000

0

-1
1
4

x 10

V

Y

(m/s)

-3000
3000

1

0

0
-1

-1

3000

0

4

x 10

V X (m/s)(e) B 50 μs

0

V Y (m/s)

0
-3000

-3000

V X (m/s)(f) La 50 μs

Figure 3.14: Evolution of the velocity distribution of B and La atoms in an Ar environment with
3 Pa and 400 K. At 0 second, particles are injected along the z axis to the positive side, with akinetic energy of 5 eV. (a), (c) and (e) are for boron and (b), (d) and (f) are for lanthanum. In
these graphs, 5000 particles are displayed.

CHAPTER 3. RESULTS AND DISCUSSIONS

80

boron atoms is faster than that of lanthanum atoms when their kinetic energies are identical.
Boron atoms experience more collision than lanthanum in a same interval, and are decelerated
more.
Next we show the transport of these atoms in the real space. Figure 3.15 shows how thesehigh energy particles spread in the gas environment. Conditions are the same with ﬁgure 3.14.
B or La atoms are injected at 0 second from the origin toward the (0,0,1) direction with theenergy of 5 eV. The pressure and the temperature of Ar gas are 3 Pa and 400 K, respectively.
Two clear diﬀerences are observed between boron and lanthanum in their transport properties; one is whether the centroid of the distribution shifts from the origin or not, and theother is the rate of dispersion. The propagation of lanthanum atoms is strongly biased to theregion z > 0 while the area where they spread is rather limited. On the other hand, the boron
atoms propagate into not only z > 0 but also z < 0, and are dispersed promptly in the argon
atmosphere.
These behaviors in real space reﬂect the time evolution of the velocity distribution shownin ﬁg. 3.14. Boron atoms are scattered in its direction before they are decelerated, hence they
spread in real space homogeneously. On the other hand, lanthanum keeps the initial momentum
along the z axis, they seem to be impinged into the z > 0 initially, and diﬀuse from there. Since
we ﬁx the initial energy in both cases, the lighter boron is dispersed more quickly. It also should
be noted that the mean free path becomes larger for faster particles as shown in section 2.5.1,
which may also aﬀect the faster dispersion of boron.

3.4.4

Discussion and summary

Now we consider the reason for the compositional deviation of LaB6 sputtered ﬁlms based on
the discussion in the previous subsection.
In low pressures, where the ballistic transport is dominant, the boron content decreases
gradually with increasing the gas pressure. This is originated from the diﬀerence in the scattering
behavior with argon gas atom; boron may be backscattered, while lanthanum is not. As thepressure increases and more collisions occur in this region, arrival probability of boron onto thesubstrate decreases more rapidly, because they tend to be deﬂected more by the collision thanlanthanum.
On the other hand, at higher pressures, B/La ratio recovers to six (the stoichiometric value)as the gas pressure increases. It is explained by considering the thermalization of sputteredatoms. As described in section 3.3, sputtered atoms are more likely to be decelerated as thegas pressure increases. When their energy decreases down to the thermal energy of the ambientgas (thermalized), their transport is no more ballistic but diﬀusive. If, in this stage, the thermalization point (i.e., source of the diﬀusion) is same, the arrival probability at the substrateis also same independent of the diﬀusion coeﬃcient (see section 2.6.2). In other words, themass diﬀerence does not aﬀect the transport anymore in the diﬀusive transport. However, thethermalization proﬁle of lanthanum atoms distributes further from the target (nearer from thesubstrate) than boron, La ﬂux is still larger than B. But at high pressure limit, where the thermalization of both elements occur just above the erosion track, there is no diﬀerence between
the thermalization proﬁles of La and B, hence the ratio of atomic ﬂuxes on the target recovers
to stoichiometry.

CHAPTER 3. RESULTS AND DISCUSSIONS

81

20

Z (cm)

Z (cm)

20

0

-20
20

20

0

-20
20

-20

X (cm)

20

Z (cm)

Z (cm)-20(b) La 10 μs

20

0

-20
20

20

0

0

-20
20

20

0

0
-20

Y (cm)-20

0-20

Y (cm)X (cm)-20

X (cm)(c) B 50 μs(d) La 50 μs

20

Z (cm)

20

Z (cm)

0-20

Y (cm)X (cm)(a) B 10 μs

0

-20
20

20

0

Y (cm)

20

0

0
-20

Y (cm)

0

0

-20
20

-20

X (cm)(e) B 200 μs

20

0

0
-20

Y (cm)

0-20-20

X (cm)(f) La 200 μs

Figure 3.15: Propagation of the B and La high energy atoms in Ar gas. Conditions are the samewith ﬁgure 3.14. (a), (c) and (e) are for boron, and (b), (d) and (f) are for lanthanum. Also5000 particles are displayed.

CHAPTER 3. RESULTS AND DISCUSSIONS

82

The gas pressure where the B/La ratio shows the minimum was inconsistent between experiment and simulation. This pressure corresponds to the transition of the transport from ballistic
to diﬀusive. Thus in the simulation, the deceleration of atoms should be overestimated. As
listed in section 3.1, the gas pressure is assumed to be uniform in the simulation. But at highpressure and high discharge power, it has been proposed that the gas density becomes smaller
near the target (see section 1.1.4). Assuming this eﬀect, the gas density near the target is lower
than the gauge value in the experiment, which raises the minimum pressure. In the simulation,
the minimum pressure became higher when ejection energy was set to higher value (10 eV). It isbecause more collision is required before the thermalization. Summarizing, the ejection energyof sputtered atoms and the gas pressure may reﬂect the gas pressure at the B/La minimum. In
the simulation in this section, we simply assumed the identical (and constant) ejection energy
for both B and La. If more proper model for this ejection energy would be developed in the
future, more quantitative discussion should be possible.
In actual sputter deposition, backscattered boron atoms will deposit onto the target. Whenthey do on the erosion track, they are re-ejected again by the following (re)sputtering. But ifthey deposit on the non-erosion track part of the target, they will be accumulated there unlessthe sputtering rate is faster than the deposition rate. As predicted by ﬁgure 3.14, boron is more
likely to go back to the target surface, hence boron content may be higher there. The experimentally observed target color degradation at high pressures also suggests this backscatteringand deposition of sputtered particles, mainly boron maybe.
In summary, we have shown the strong dependence of B/La atomic ratio on process pressurein sputter deposition process of LaB6 ﬁlms from both experiments and Monte Carlo simulation
study. The ﬁlm composition is almost stoichiometry in both low and high pressure limits, butshows a minimum as low as B/La=4 in the midst pressure. The depletion in boron in lower
pressure was explained by the diﬀerence in scattering angle distribution between B and La inballistic transport stage. Since boron is more likely to be scattered backward by the collision
with argon gas atom, the arrival of it on the substrate decreases more than lanthanum as the gaspressure increases. Recovery of B/La ratio at higher pressures stems from the thermalizationand following diﬀusive transport of sputtered atoms. When the pressure is not high enough,lanthanum tends to thermalize further from the target (i.e., nearer from the substrate), so arrival
ﬂux of it is larger. On the contrary, at high pressure limit, the thermalization proﬁle of both B
and La atoms approaches to the target neighbor. Therefore, the ﬂux ratio of both elements atthe substrate is identical to the target composition.

3.5

Density of sputtered atoms in the plasma

In this section, we describe a novel method to evaluate the spatial density of sputtered particles
using Monte Carlo simulation.
While the eﬀect of energetic particles on deposited ﬁlms has intensively been studied bymany researchers, the spatial density of them has rarely been investigated. It is probably
due to the unimportance of it to the processes of conventional sputtering (in contradiction
to CVD). However, it has recently attracted attentions because of the recent development ofionized sputtering [117], as well as the fundamental interests to the magnetron plasma. Some

CHAPTER 3. RESULTS AND DISCUSSIONS

83

experiments using optical absorption [9, 25], and its variants [78] have been reported, but thetheoretical studies are rather scarce. One of such eﬀorts includes the versatile model of IPVD
process by the group at Univ. of Illinois [39]. They elucidated the particle density combining
the MC approach and diﬀusive approach. In their case, it is necessary to assume the diﬀusion
coeﬃcient of Cu in Ar gas environment.
From purely particle-oriented viewpoints, it is possible to deduce the density of sputteredparticles from the residual time of each particles and the ejection frequency of them from thetarget. For example, if 10 particles are ejected from the target per one second, and if eachparticle stays for 10 seconds in the chamber, 100 particles should be found there on the average.
Similarly, if you split the chamber volume into cells and evaluate the residual time in each ofthem, you can obtain the atomic density proﬁle. Also note that the ejection ﬂux of the sputteredatoms from the target should be proportional to the incident ﬂux on other chamber boundaries,
including substrates. Therefore, if this proportional constant is known, you can calculate thedensity from the residual time and the deposition rate.
The residual time of sputtered particles can be calculated by the Monte Carlo simulation
method. To do this, the MC method must be extended so that it can treat properly the motionof particles after they are “thermalized”, where their energy become comparable to the thermalenergy of gases. In other words, the colliding gas motion must be incorporated into the MC
framework, or the sputtered atom is monotonically decelerated and stops its motion after somecollisions (as discussed in section 3.2). Our model described in 2.5 is well suited for this purpose;the motion of sputtered atoms in this case approaches the Maxwellian distribution in this case(see ﬁgs. 3.4 and 3.5).
The introduction of the eﬀect of thermal gas atoms has partially been done by Turner et
al. in 1989 [153] (see section 1.6.2). They applied their method to the gas density calculation,and showed that the gas density increased with increase in gas pressure. However, the absolutevalues of them were not shown since the ejection frequency is not taken into account.
This method is applied to the DC magnetron plasma of Cu/Ar. The calculation result iscompared with the previous reports as well as with our OES measurement of this magnetronplasma. For the latter, the “self absorption” eﬀect of a certain line from the Cu neutral is
discussed. This eﬀect is caused by the reabsorption of the emitted photon by another copper
neutrals, and reﬂects the spatial density of them in the plasma environment.

3.5.1

Density calculation

To calculate the residual time of sputtered particles in the chamber, MC simulation describedin 3.1.2 is used, except for the BEM treatment of diﬀusion process. The motion of each atomis throughly followed until it reached the chamber wall. The gas temperature is assumed to be
400 K, and gas density is assumed to be uniform in the chamber.
Chamber boundary used in this MC simulation is shown in ﬁgure 3.16. It is basically thesame with the one used in section 3.3. The radius of the chamber is set to 10 cm. In addition,
the region between the target and the substrate holder was also split into meshes at intervals of
5 mm, in each of which the atomic density was calculated.
At ﬁrst, the residual time of sputtered particles are calculated. Since the speed of the particle
between two collision events is constant, the time can be deduced by the segment length of theCHAPTER 3. RESULTS AND DISCUSSIONS

84z

Substrate

Targetr
Figure 3.16: Chamber model used in the MC simulation for atomic density evaluation: upward
arrow at the left is the symmetric axis. Chamber boundary is split into nodes, whose ends areshown as black dots. The region between the target and the substrate holder is also split into
meshes at intervals of 5 mm for both r and z axes.

CHAPTER 3. RESULTS AND DISCUSSIONS

85

particle path, cut by the cell boundary. The time in each meshed cell was summed up for all
the particles. Now let us pick one of these cells up, and denote the residual time of it as TC .
If we can elucidate the time TM C corresponding to the ejection of whole particles in the MC
simulation, the atomic density nC at the cell is given by:
nC =TC
TM C VC(3.5)

where VC is the volume of the cell.
With the same simulation, in which TC is calculated, numbers of arrival particles are recordedon each boundary nodes. They are divided by the area of these nodes, and the areal density ofarrival particles ρMC are obtained. By comparing ρMC with the experimentally obtained atomic
ﬂux there, the TM C can be derived. If the atomic ﬂux per unit time and unit area is measured
to be FEXP , TM C is given by:(3.6)
TM C = ρM C /FEXP .
Substituting this relation into eq. 3.5, nC can be obtained.
The simulation was performed for 104 ∼ 105 particles, and was repeated several times with
diﬀerent seeds of pseudo random number sequence to evaluate the scatter of the simulation
result. It was found that typical value of the error in the residual time was less than 1%, andthe error in the deposition ﬂux at the substrate was 6% at most.
The experiment to measure the deposition rate was identical with the one described in
section 3.3. Target was 8 cmφ disk of 4N Cu. Substrate holder was positioned 5 cm above
the target. Argon was used as a discharge gas, and its ﬂow rate was kept constant at 10 sccm.
Pressure was changed between 2 ∼ 10 Pa by throttling the evacuation valve. DC power wasapplied to the Cu target cathode to generate a glow discharge. DC power was ﬁxed 100 W
throughout the experiment. Film thickness on the substrate holder was measured by multiple
interferometry, and deposition rate was obtained by dividing it by the deposition time.
Deposition rate is converted into the atomic ﬂux as follows. Atomic density of copper bulk
is given using copper bulk density ρ = 8.93 g/cm3 , copper atomic weight U = 63.546, and theatomic mass unit μ = 1.66054 × 10−27 kg:
nCu = ρ/(U μ)= 8.93 (g/cm3 )/ 63.546 × 1.66054 × 10−24 (g)= 8.46 × 1022 (atoms/cm3 ).

With this, deposition rate given in nm/s becomes:
1(nm/s) = 10−7 (cm/s) × 8.46 × 1022 (atoms/cm2 s)= 8.46 × 1015 (atoms/cm2 s)(3.7)

Figure 3.5.1 shows the pressure dependence of the deposition rate obtained by the experiment. Plotted data were measured at the center of the substrate holder. Deposition rate wasalmost constant until 10 Pa, and decreased at 20 Pa. It was as large as 0.1 nm/s, which corresponded to about 1015 atoms/cm2 s. It should be noted that the atomic ﬂux proﬁle by thissimulation has reproduced the thickness proﬁle by experiment fairly well as shown in section 3.3.

Deposition Rate (nm/s)

CHAPTER 3. RESULTS AND DISCUSSIONS

86

0.2

0.1

0
1

5

10

Pressure (Pa)
Figure 3.17: Pressure dependence of deposition rate at the center of the substrate holder (experimental result).

CHAPTER 3. RESULTS AND DISCUSSIONS

87

Figure 3.18: Density proﬁle of Cu atoms. Black lines denote the target and the substrate holder.
Scale is shown at the top in logarithm.

CHAPTER 3. RESULTS AND DISCUSSIONS

88

Using the deposition rate in ﬁg. 3.5.1 and the residual time / the arrival atom density
obtained by the simulation, spatial proﬁle of copper atomic density was obtained at four diﬀerent
gas pressures. The results are shown in ﬁgure 3.18. Black lines in the ﬁgure represent a partof chamber boundary shown in ﬁgure 3.16. The Dot-dash line at the left of each proﬁle is thesymmetrical axis. Ejection of sputtered particles occurs, as shown in ﬁg. 3.2, at around 2 cm
away from the axis on the target surface, which is placed at the bottom of the proﬁle. The color
scale of the density is located at the top, with which the logarithm of atomic density is mapped
to the color between 1016.5 ∼ 1019 m−3 .
It is observed that the atomic density increases as the pressure of the discharge gas increases.
At higher pressures, mean free path of the sputtered particle becomes shorter. Moreover, sputtered particles are quickly decelerated by collisions with gases and thermalized at high pressures,
which also contributes to reduce the mean free path (see section 2.5.1). It reduces the diﬀusion
coeﬃcient of sputtered atoms 3 , which result in the higher density. At each pressure, the atomic
density is largest at the region nearest to the erosion track, which is the source of particles.
On the other hand, the density decreases steeply as the chamber boundary approaches. We
have assumed in the simulation that the particle arrived at the boundary adsorbs there with the
sticking coeﬃcient of unity. In the diﬀusion transport model, it corresponds to the boundary
condition of zero density. The rapid reduction of the density near the wall at higher pressuresmeans that most atoms are in the diﬀusive stage, rather than in a ballistic one. If, as an extreme
case, particles do not experience any collision until they reach the wall, the atomic density proﬁlestem just from a spread of particle ﬂux as they come away from the target, and should not show
the rapid decrease near the chamber wall.
The atomic density of 1017 ∼ 1018 m−3 agrees in order with the experimental results by Donyet al. [25]. They sputtered aluminum by argon in an RF magnetron sputtering apparatus. At
argon pressure of 10.6 Pa, they have yielded 2.2 × 1017 m−3 with the RF power of 20 W, and
4.8 × 1017 m−3 with 30 W. Our result seems slightly larger than theirs, even considering thediﬀerence in discharge power. It may be due to the larger sputtering yield of Cu than Al (seereference [73] for example), or the gas rarefaction eﬀect [115,152] we neglected in the simulation,
which reduces the gas density near the target and hence reduces the atomic density of sputteredparticles.

3.5.2

Optical emission measurements

Optical emission (OE) from the plasma discharge reﬂects many important characteristics ofit, such as the density of neutral/ion species, electron density and its energy distribution [82].
Therefore, Optical Emission Spectroscopy (OES) measurement has been recognized as a suitable
method for plasma diagnostics.
There have also been many reports applying OES to the magnetron sputtering plasma. For
example, dependence of emission intensities on discharge power [82] or gas pressure [26,27] havebeen studied, and the excitation process of neutral / ion species were discussed. OES has also
been applied to the observation of magnetron plasma, where the proﬁle of the power injectionis hardly predicted because of the spatially inhomogeneous magnetic ﬁeld [84, 132]. As for a3

see chapter XIII of reference [56].

CHAPTER 3. RESULTS AND DISCUSSIONS

89

rather special case, self-sustained mode in magnetron sputtering where discharge is maintained
without argon supply has been observed by OES [110]. Since the enhancement of ionizationof depositing materials has recently attracted attention in the ﬁeld of plasma processing for
microelectronics [117], OES will be used more in the future.
In this study, we have developed an OE acquisition system based on a commercially available
reﬂex camera, and performed the OE measurement with spatial resolution. A single plano-convex
lens made of synthesized quartz with a focal length of 55 mm is used for imaging. The lens canfocus from 20 cm to inﬁnity with this system. An iris is placed before the lens, whose diameter
was ﬁxed at 30 mm in this experiment. At the ﬁlm plate of the camera, a receptacle of an opticalﬁber is mounted, and a quartz optical ﬁber is connected to it. With this optical system, we canset the focal point easily as we look the object through a view screen of the ﬁnder. The otherend of the quartz ﬁber is connected to the entrance of the “polychromator (Soma Optics Ltd.,
S-2210)” which consists of a grating and 1024 channel CCD linear array sensors, which enables
an uptake of whole spectrum simultaneously. The optical emission can be observed with a rangeof 185∼525 nm with 0.34 nm resolution.
With this measurement system, Cu sputtering plasma was observed through a quartz viewport mounted on the sputter chamber as shown schematically in ﬁgure 3.19. The data wereacquired by opening the shutter, the time of which was set to 1 ∼ 2 s in this study. The degradation of the transmittance by the view-ports could be neglected in this experimental condition.
OES spectra were measured by focusing the camera on the position 5 ∼ 10 mm above the centerof the target.
Typical OE spectra from the argon DC magnetron plasma with a copper target are shownin ﬁgure 3.20. The unit of the vertical axes of these ﬁgures is output of A/D converter of CCD
array sensors. Emission lines have been identiﬁed using the ATMLINE, atomic line databaseprovided by Kyoto University Data Processing Center [48]. Several emission lines could be
assigned to Cu neutrals, Ar neutrals and Ar ions, and a part of them is listed in Table 3.1. Thelines chosen as representatives of respective species are marked as (*), whose dependences ondischarge parameters will be discussed in the following. Though it has been reported that Cu
ion lines are observed at 213.6 nm [112] and at 404 nm [110], such peaks could not be separated
clearly from other peaks around them in our system.
The Ar pressure dependence of intensities of (*) lines in Table 3.1 are shown in Fig. 3.21.
The upper ﬁgure shows the results at the focal point of 5 mm above the center of the target,
while the lower one shows that of 10 mm. DC discharge power was ﬁxed at 100 W. The
intensities of emission lines from Ar neutral/ion (420.1, 488.0 nm) are not strongly inﬂuenced
by the Ar pressure, while those of the lines assigned to Cu neutral (324.8, 510.6 nm) increases
with increasing the Ar pressure. These behaviors of emission lines agree with the reports on themagnetron plasma [26, 27], and on the ECR plasma [112].
In general, the electron density of the glow discharge plasma (including magnetron discharges) increases with increasing gas pressure. On the other hand, the electron temperature
decreases with increasing pressure. Since it requires more than 10 eV to excite Ar atoms, theintensity of Ar emission does not increase even if electrons with medium and lower energies
are generated steadily at high pressures. On the other hand, Cu neutrals can be excited withenergies of 3 ∼ 4 eV, their emission intensities become stronger with increasing gas pressure. In
addition, the spatial density of Cu atoms in the plasma becomes larger at high pressures, which

CHAPTER 3. RESULTS AND DISCUSSIONS

Figure 3.19: OES observation system.

90

CHAPTER 3. RESULTS AND DISCUSSIONS

2000

91

1 Pa

Intensity (counts)

1000

0

2000

20 Pa

1000

0

200

300

400

500

Wavelength (nm)Figure 3.20: Optical emission spectra from DC sputter plasma with Cu target and Ar gas. Left:
Ar pressure = 1 Pa, DC Power = 100 W, DC current = 0.15 A. Right: Ar pressure = 20 Pa,
DC Power = 100 W, DC current = 0.26 A.

CHAPTER 3. RESULTS AND DISCUSSIONS

92

103

Intensity (counts.)

102

10

103

5 mm
324.8 nm (Cu I)
510.6 nm (Cu I)
420.1 nm (Ar I)
488.0 nm (Ar II)

102

10

10 mm
0.5 1
5 10
Pressure (Pa)

Figure 3.21: Pressure dependence of optical emission intensities from Cu/Ar magnetron sputter
plasma.

CHAPTER 3. RESULTS AND DISCUSSIONSType
Cu I (*)
Cu I
Cu I (*)
Cu I
Ar I (*)
Ar I
Ar I
Ar II
Ar II(*)

Wavelength (nm)
324.8
327.4
510.6
521.8
420.1
425.9
427.2
442.6
488.0

Upper Level (eV)
3.817
3.786
3.817
6.192
14.49914.73514.52519.54919.68093
Lower Level (eV)
0
0
1.389
3.817
11.54811.828
11.62316.74917.140

Table 3.1: Emission lines assigned by ATMLINE. (*) lines are displayed in the following ﬁgures.
also contributes to the increase of Cu emission intensities. The Cu density increase also brings
the consumption of the electrons of higher/medium energies, which may result in the “colder”
plasma at higher pressures [112].
In ﬁgure 3.21, it is also shown that the intensity of 510.6 nm copper line increases much more
steeply than 324.8 nm line with increasing pressure. And by comparing the results of 5 mmwith those of 10 mm, ratios of argon intensities at 5 mm and at 10 mm are almost constantat any pressure, while the copper intensities at 10 mm increase more rapidly than at 5 mm as
the pressure increases. Normally the plasma should be conﬁned nearer to the target at higher
pressures, hence it is interesting why the emission from 5 mm becomes less stronger in this case.
The reason of these phenomena will be discussed in the next subsection.
Power dependences of OE intensities for several gas pressures are shown in Fig. 3.22. The
observation position of these is 5 mm above the target. OE intensities of Ar lines are roughly
proportional to the DC power P , while those of Cu are proportional to P 2 . The latter may
reﬂect: (1) ﬂux of sputtered Cu atoms are proportional to P and (2) the excitation probability
in the plasma is proportional to electron density n, which is also proportional to P . As consistent
with Fig. 3.21, the peak intensity of 510.6 nm Cu line is much smaller than 324.8 nm Cu line at
1 Pa, but becomes comparable at 20 Pa. At 20 Pa, the 510.6 nm line intensity increases more
rapidly than the 324.8 nm line with increasing the discharge power.

3.5.3

Self absorption eﬀect of optical emission

In this subsection, we discuss the discharge parameter dependence of two emission lines fromCu neutrals, namely the lines at 324.8 nm and at 510.6 nm. The diﬀerence in pressure, power
and positional dependence of them are considered relating to the atomic density of the copper.
As already shown in table 3.1, these two lines share the initial state at 3.817 eV. More

CHAPTER 3. RESULTS AND DISCUSSIONS

94

324.8 nm (Cu I)
510.6 nm (Cu I)
420.1 nm (Ar I)
488.0 nm (Ar II)

103
102

Intensity (counts)

10

1 Pa

1
103
102
10

5 Pa

1
103
102
10
1
10

20 Pa
50

100

Power (W)
Figure 3.22: Discharge power dependence of optical emission intensities from Cu/Ar magnetron
sputter plasma.

CHAPTER 3. RESULTS AND DISCUSSIONS

95

concretely, their initial and ﬁnal states are 4 :324.8nm : 4p 2 P◦1 1/2 → 4s2 S1/2
510.6nm : 4p 2 P◦1 1/2 → 4s2 2 D2 1/2
It should be noted here that the initial state of these two lines are identical, and that theﬁnal state of the former line is the ground state so that it is absorbed again by ambient Cu
atoms at ground state. Assuming that the temperature of the electron system of copper atoms
is comparable to the electron temperature of the plasma, most of the copper atoms should beat ground state. Hence, 324.8 nm emission can be reabsorbed by another copper atoms (atground state) existing in the plasma, which is known as “self absorption eﬀect”. This meansthat the plasma with high atomic density of Cu is opaque for the light of 324.8 nm, while it isstill transparent for the 510.6 nm lines since Cu atoms at the excited level 4s2 2 D2 1/2 are muchless than those at ground state. Since the probabilities of these transition should be unchanged
unless the induced emission occurs, the relative intensity of 510.6 nm light emission increases as
the copper density increases.
According to the atomic absorption theory [162], the integration of the absorption spectrum
κ(ν) for atoms with density N is given by:

πe2
N f.
(3.8)
κ(ν)dν =mc
where f is an absolute oscillator strength, which is about 0.434 for 324.7 nm spectrum line ofCu neutral. Since the pressure is low in this case, the collisional broadening which gives theLorentz proﬁle should be weak. Hence the Doppler broadening rules the absorption spectrum,and it becomes a gaussian with the width Δν of:
Δν =

ν0
c

2 ln 2 kT
.
M

(3.9)

If we substitute this gaussian proﬁle into eq. 3.8, the inverse of the absorption coeﬃcient at
the peak of the spectrum is:
6.64 × 1014
1=[m]
(3.10)
κ(ν0 )
N [m−3 ]at 300 K. It gives the decay length of the 324.7 nm emission. At the Cu density of 1017 ∼ 1018 ,which has been obtained previously by the simulation, decay length is 6.6 ∼ 0.66 mm. It agrees
the more gentle increase of 324.7 nm emission intensity at higher pressures where the Cu density
is high, compared to the increase of unabsorbed emission at 510.6 nm.
This optical reabsorption process can also explain why the Cu intensities at 10 mm above
the target increase more rapidly than at 5 mm with increasing the Ar pressures. The emission
originated from nearer to the target must travel the thicker layer of dense Cu atmosphere sothat the escape probability is smaller. In addition, copper atoms at 10 mm above the target
may be excited by 324.8 nm photons originated from the dense plasma near the target, which
also reduce the diﬀerence of emission intensity.
4

Identiﬁcation of the initial/ﬁnal states is based on the ATMLINE database at Kyoto University Data Processing Center [48] and NBS databook compiled by Moore [86].

CHAPTER 3. RESULTS AND DISCUSSIONS

3.5.4

96

Summary

In summary, a novel method to evaluate the density proﬁle of sputtered particles has beendeveloped using Monte Carlo simulation. It has come to possible by the introduction of the gas
motion into the MC, which is one of the achievements of this study. The method was applied
to the Cu/Ar magnetron plasma, and the Cu atomic density of 1017 ∼ 1018 m−3 was obtainednear the target, with DC discharge power of 100 W and gas pressures of 2 ∼ 20 Pa. The density
increases with increasing the Ar pressure, and the spatial proﬁle of Cu density can be explained
considering the diﬀusion process of sputtered Cu atoms.
The same Cu/Ar plasma was observed by OES. It was found that the two emission lines from
Cu atom (324.7 nm and 510.6 nm) showed diﬀerent dependence on discharge parameters. The
former is found to be reabsorbed by the ambient copper atoms at ground state, and the decay
length of the 324.7 nm emission is deduced using the density value obtained by the simulation.
As a result, some millimeters to submillimeters of decay length is obtained, which can explain
the behavior of optical emission lines from Cu atom.

Chapter 4

Conclusion
In this thesis, we have studied the transport process of sputtered atoms from target to substrate,with both experiment and with computer simulation. Two new models have been proposed toextend the conventional Monte Carlo method, for more realistic simulation in high gas pressure
environment. These models were implemented into the FORTRAN program code, and thecomputer simulation was performed. By comparing the simulation results with experiments,
the validity of the model could be conﬁrmed. The reasons of the experimentally observed
phenomena have also been discussed from the viewpoint of particle transport process, and newunderstandings have been obtained.
In chapter 1, the background of this study was described. In sections 1.1–1.5, various problems relating to the sputter deposition process were reviewed. Theoretical reports and considerations were mainly treated, and relationships of them with the particle transport process was
also discussed. In the next section (1.6), historical review of the studies on Monte Carlo (MC)simulation of particle transport process was presented. Latest achievements in this ﬁeld werecovered as well as possible, and the position and the signiﬁcance of this study were shown.
In chapter 2, physical backgrounds of the simulation model were described. In sections 2.1–
2.4, conventional MC models for particle transport problem were presented, and those whichused in this study were mentioned in detail. The following two sections are the ones originallydeveloped in this study. In section 2.5, we have introduced the thermal motion of ambient gas
atoms into collision problem of the sputtered particles, and deduced the velocity (i.e., speed andangle) proﬁle of the colliding gas atoms and their dependencies on the particle speed. And insection 2.6, we have developed the method to treat the random diﬀusion transport of sputteredparticles using Poisson’s equation. The numerical method to solve this equation using theBoundary Element method (BEM) has also been described. By combining the MC and BEM
methods, it was proved that the particle transport process could be solved in applicable periodsof time with commercially available personal computers.
In chapter 3, we have applied the simulation to several problems in sputter deposition andshown the results. Firstly in section 3.1, we described experimental setup, with which thedeposition experiment was performed and was compared with the simulation. Subsequently, thesimulation methods, which had been explained in chapter 2, were summarized again.
In section 3.2, we applied the simulation to the deceleration process of high energy copper

97

CHAPTER 4. CONCLUSION

98

atoms injected into argon gas environment without boundary. Especially, eﬀect of the thermal
motion of gas atoms were investigated by comparing the cases with gas temperatures of 400 K
and 0 K. In the former, energy distribution of Cu atoms well reproduces the Maxwellian, whilein the latter, gas atoms were immobile, it approaches zero since Cu atoms lose their energy on
collision monotonically. It governs the time evolution of the spread of Cu atoms in real space;
if the gas motion is neglected, the dispersion of Cu atoms eﬀectively stops at some size. For
example, this size is the order of centimeters when Ar pressure is 10 Pa. As a result, it has been
proved to be important to consider the thermal motion of gases when you simulate the particle
transport process in high pressure environment.
In section 3.3, thickness proﬁle of the sputtered copper ﬁlm has been investigated. It was
shown that the proﬁle changed by the distance between the substrate holder and the chamberwall. When the distance is small, ﬁlm thickness decreases more steeply as apart from the holder
center and close to the wall. It could be explained by the model introduced in section 2.6, where
the chamber wall is regarded as the boundary with zero atomic density. When you prepare theﬁlm by sputter deposition in high pressure, it also deposits on the faces from which the target
cannot be seen directly. It is because of the diversion of sputtered atoms caused by the scattering
by gas atoms. It had been reported that the ratio of ﬁlm thickness on substrates facing toward/ away from the target (BFR ratio) showed the maximum at midst gas pressure [90]. It wasalso conﬁrmed with our experiments that this BFR ratio decreased and approached constant at
higher pressures where the diﬀusive transport became dominant. Using the simulation developedin this study, this tendency was reproduced well, and the reason was elucidated as follows. With
the second extension of the simulation using the Poisson’s equation, thermalization points ofthe sputtered atoms are treated as the particle sources. At pressures where the BFR ratio ismaximum, these points distribute around the substrate holder so that the ﬂux at the backface
increases relatively. On the other hand, at very high pressures, thermalization points concentrate
to the neighbors of the erosion track resulting in the decrease of the relative ﬂux at the backface.
In section 3.4, the pressure dependence of the ﬁlm composition was considered in case of
LaB6 sputter deposition. When the process pressure changed from 0.2 to 20 Pa, boron content
showed the minimum at around 7 Pa. It is because of the diﬀerence in scattering angle with
Ar gas atom. Boron, the lighter atom, is more likely to be backscattered by the argon atom,
hence the pressure increase in the ballistic transport region (i.e., low pressures) results in thedecrease of boron content. On the contrary at high pressures, where the diﬀusive transport isdominant, mass diﬀerence scarcely aﬀect the particle transport, and the composition recovers
to the stoichiometry. That is to say, the competition between the ballistic transport and thediﬀusive transport is the cause of the peak in pressure dependence.
In section 3.5, we have proposed the method to evaluate the atomic density in the sputtering
plasma using MC simulation. This has been achieved by the ﬁrst extension of MC in this study,
with which the motion of sputtered atom can be traced even after the many times of collisions
with gases. The proposed method was applied to the same situation with section 3.3 (Cu deposition), and copper atomic density was calculated to be 1017 ∼ 1018 m−3 . The sputter plasma
in this condition was measured experimentally using optical emission spectroscopy equipment.
It has been observed that the copper emission line, whose transition includes the ground stateas the ﬁnal one, becomes small when Cu atomic density is high (e.g., high gas pressure and/or
high discharge power). We have assumed that this is because this emission is re-absorbed by theCHAPTER 4. CONCLUSION

99

ground state copper atoms in the plasma, and have estimated the decay length of this emission
using the atomic density value obtained by the simulation. With it, we could successfully explain
the emission intensity reduction observed by the experiment.

4.1

Future works

By closing this thesis, we would like to point out the limitations in the simulation developed here,
and to suggest the desirable researches and developments in the future. Among the conditions
and assumptions summarized in section 3.1, we consider that four of these are not always
adequate. These are listed in the following,• In this simulation, the ejection position of the sputtered particle is determined so as toreproduce the depth proﬁle of the erosion track actually measured after the deposition
experiment (2.2.1). Therefore, it cannot be considered as the a priori model applicableto the design of sputter deposition equipments. Moreover, the experimentally measured
erosion track depth is actually the results of many times of deposition experiments with
diﬀerent discharge conditions. If, for example, the gas pressure changes, the ion production
proﬁle (hence the ion ﬂux onto the target) may also change. To overcome this problem,the simulation should incorporate such models that can calculate the ion production fromthe magnetic ﬁeld conﬁguration, as mentioned in section 1.1.2.
• The angle and energy distributions of sputtered atoms are simply assumed to be the cosine
distribution and the Thompson’s formula, respectively. This is partly because that there
are no satisfactory models for roughened targets after some time of the sputtering, nor fortargets of compounds or alloys. As referred in section 1.2, there have been some calculation
codes treating the collision cascade inside of the target with Monte Carlo method. You
should check out and validate them by experiments, and combine them with particle
transport simulations.
• The density and the temperature of the sputter gas are assumed to be uniform in thechamber. As mentioned in section 1.1.4, gas density reduction eﬀect, caused by the highenergy ﬂux near the target (neutralized and recoiled gas ions, sputtered atoms, etc.), hasbeen proposed. To treat this, gas dynamics has to be taken into account by DSMC method(appeared in section 1.3.3) or others.
• The sticking coeﬃcient of sputtered atoms on the chamber wall is assumed to be unity,
while it actually should be (slightly or considerably) less than 1. For the latter, newboundary condition must be introduced to the BEM treatment of the diﬀusion process of
atoms.
A part of these proposals has already been performed, which are mentioned in section 1.6.4.
It is desirable that these eﬀorts will be merged into the generic and comprehensive calculation
model of the sputter deposition process in the future.

Appendix A

Implementation
In this appendix, we describe the simulation program used in this study. At the ﬁrst section,some of the important algorisms and libraries are mentioned. In the second section, detailedof the program package is explained: availability and the usage of the program, our computing
environment, contents in the distribution archive, and the structure of the program.

A.1

Algorisms and libraries

A.1.1

Integration

As shown in section 2.3.2, the integration which gives the scattering angle has a singularity at
the lower limit of its range (eq. 2.40). In our program, we avoid this problem by using the
Double Exponential (DE) formula for numerical integration.
In DE formula, the integration
∞f (x)dx,(A.1)

x = exp (2 sinh t).

(A.2)

I=
0

is converted using
After this conversion and the discretization of the integration, we obtain the summation form:
Ih = h∞f (exp(2 sinh(ih))) 2 cosh(ih) exp(2 sinh(ih)).

(A.3)

i=−∞

The conversion given by eq. A.2 uses the fact that the trapezoidal rule gives extremely
accurate result for the integration of analytical function if the integral region is inﬁnite for both
positive and negative ends. In addition, the function in the summation of eq. A.3 approaches
to zero double exponentially for −∞ and ∞.
In the program, we have another integrals given by (eqs. 2.141 and 2.142). They are theintegrals of complete elliptic functions to construct the coeﬃcient matrix of the BEM method.
For these, we have used Newton Cotes (NC) method of 8th order, instead of the DE, becauseof the better convergence.
100APPENDIX A. IMPLEMENTATION

101

The Fortran source codes of DE and NC methods have been published in reference [87] withdetailed descriptions of theoretical backgrounds. In the courtesy of Professor Mori, they havebeen approved to use and to distribute with this program package.

A.1.2

Random number generation

For Monte Carlo simulations, it is critically important to use a pseudo random number generator
that has a long period as well as good “randomness” (for the measure of randomness, see [66]for example).
In this study, we have adopted the “Mersenne Twister (MT)” random number generator. It
has many preferable features, especially for MC calculation. Following is the quotation from the
Mersenne Twister Home Page at:
http://www.math.keio.ac.jp/~matumoto/emt.html
Mersenne Twister(MT) is a pseudo random number generator developed by Makoto
Matsumoto and Takuji Nishimura (alphabetical order) during 1996-1997. MT has
the following merits:• It is designed with consideration on the ﬂaws of various existing generators.
• The algorithm is coded into a C source downloadable below.
• Far longer period and far higher order of equidistribution than any other implemented generators. (It is proved that the period is 219937 − 1, and 623dimensional equidistribution property is assured.)• Fast generation. (Although it depends on the system, it is reported that MT issometimes faster than the standard ANSI-C library in a system with pipeline
and cache memory.) Eﬃcient use of the memory. (The implemented C-code
mt19937.c consumes only 624 words of working area.)
The C program is also available from the above url. We have modiﬁed it slightly to make itcallable from our Fortran routines. Reference [71] has given us a lot of information about theusage of Fortran program in Unix, as well as to build its interface to C.

A.1.3

Spline function

In our simulation, we have to utilize some functions that need large computational complexity, including some kinds of inverse functions which cannot be solved analytically. To reducethe computational load for these, we introduce the spline function [20] to approximate thesefunctions.
Spline function is a set of k-th polynomials, each of which has its own domain given by
splitting the domain of the original function. On each domain boundary, the 1st ... (k − 1)-thderivative of the polynomials of adjoining domains are requested to be continuous. A sequenceof points (“nodes”) are selected in the original domain, and values of the original function (which
may need complex calculation) are evaluated at each point. Then the factor of the polynomials
are determined so as to let the constructing “spline function” cross these points simultaneously.

APPENDIX A. IMPLEMENTATION102

In the following, we mention the applications of the spline function in this study.
• The function to represent the erosion track on the target. As described in section 2.2.1,
ejection positions of sputtered atoms are determined so as to reproduce the depth proﬁleof the erosion track. This depth proﬁle function is achieved by the 1D spline function
using the actually measured points shown in ﬁgure 3.2.
• The mean free path of the sputtered atom as a function of its speed (eq. 2.86). As shownin ﬁg. 2.4, the nodes are chosen at the same interval on the logarithmic axis of the kinetic
energy of the atom. The 1D spline function is constructed using these nodes, and is utilized
in the simulation.
• The speed of the colliding gas atom as a function of the sputtered particle speed vp andpseudo random number w (eq. 2.88). This is the 2D spline function, where the ﬁrst
parameter is taken to be the log of the sputtered particle energy log Ep , and the second
parameter is w. Nodes are taken at the same interval on both of these axes, and the spline
function is constructed.
However, since the spline function is the polynomial by its nature, it does not give agood approximation for eq. 2.88 at w ∼ 0 and at w ∼ 1. Therefore, we give the range[wmin , wmax ] only in which the spline function apply. At w < wmin , we deduce the vg usingw ∝ vg 3 , which is the asymptotic function of eq. 2.80 at vg << 1, while at w > wmax , wedirectly solve eq. 2.81 by the bisection method to obtain vg .
• The scattering angle χ0 as a function of collision parameter b and the relational kinetic
energy Er ≡ μv0 2 /2 in the potential scattering problem (eq. 2.40). The one parameter isb while the other is log Er . The nodes are taken on both axes at the same interval, andthe 2D spline function is deduced.
Since this is the even function for b, dχ0 /db is zero at b = 0. It is possible to incorporate
this restriction when constructing the spline function, but it has not been implemented
yet.
The 1D and 2D spline functions used in this study has been implemented from scratch basedon the reference [141].
When you construct the spline functions, you have to solve the linear simultaneous equations
with the rank equal to the number of nodes. Since the coeﬃcient matrix of this simultaneous
equations is sparse, it is solved using the SPARSE library, which has been developed for solving
large sparse linear systems. The SPARSE library can be obtained from:http://www.netlib.org/sparse/
Though this is written in C, the procedure to use it from FORTRAN is also documented in theirdistribution archive. In our program, it is utilized after applying this modiﬁcation.

A.1.4

Matrix calculation

In the boundary element method, you must solve the linear simultaneous equation given ineq. 2.145. It is done by the LAPACK library [4] which is a comprehensive FORTRAN library

APPENDIX A. IMPLEMENTATION

103

for versatile linear algebra operations. It has been used in various ﬁelds for scientiﬁc calculations,
including Finite Element Method, etc. The source code is in the public domain and is availableat:
http://www.netlib.org/lapack/

A.1.5

Special functions

In section 2.6, we have used the complete elliptic functions of 1st and 2nd kind on calculating thefundamental solutions in axisymmetrical Poisson’s problem (eqs. 2.134 and 2.135). The routines
to calculate these functions have been taken from the SLATEC Library, which can be found at:
http://www.netlib.org/slatec/
SLATEC is also a well known library written in FORTRAN, which includes over 1400 ofroutines for scientiﬁc calculations. The result of the complete elliptic functions are checked withtables in reference [88] before applying the calculation to our simulation.
Another special function used in our simulation is the Gaussian error function (eq. 2.83). The
g77 Fortran compiler we used (see below) has this routine in intrinsics as one of non-standard
extensions, and we use it. If your Fortran compiler does not have this function internally, you
can use the routine derf.f in the SLATEC.

A.2

Simulation Program

A.2.1

Availability and requirements

We named the simulation program used in this study as HISPUT. The archive including thewhole sources is available from our web site at
http://surf.ap.seikei.ac.jp/HISPUT/
in a tar+gzip’ed archive. It has been developed using the CVS (Concurrent Versions System 1 ).
The access via Anonymous CVS will be available in the near future.
The sources are written mainly in Fortran, and partly in C. CPP preprocessing language isused also in some Fortran sources to enable/disable debugging codes, to modify maximum sizesof various arrays, and to specify basic physical constants. Two routines are written in C. One
is the MT random number generator (see above), and the other is for the speciﬁer of the FPU
behaviour on error during the calculation.
We have used GCC (GNU Compiler Collection 2 ) to compile our program. It currently
contains C and Fortran compilers (which we have used), as well as ones for C++, Objective-C,
Java, etc. It has been proved to be matured and stable during the long period of its development.
The compilation/linking procedures are automated using “make” utility 3 . GNU version of make1

http://www.cvshome.org
http://gcc.gnu.org/
3
http://www.gnu.org/software/make/2

APPENDIX A. IMPLEMENTATION

104

has some extended syntaxes, and we use them. Therefore, automatic build may not work withtraditional make utility.
On developing and executing the programs in this study, we have used Debian GNU/Linux
operating system 4 running on IBM PC. Debian GNU/Linux is a freely available UNIX likesystem. It comes with more than 3000 software packages, which includes C/Fortran compilers,scientiﬁc libraries, graphic tools, etc. Version 2.2 (codename potato) is mainly used with someextra packages installed, e.g.:• cpp, g77, gcc• make• lapack, lapack-dev, blas1, blas-dev• slatec, slatec-dev
The last two items are packages of LAPACK and SLATEC libraries, respectively.

A.2.2

Directory structure of the distribution archive

In this subsection, we mention the contents of the distribution archive by its directory structure.
The ﬁles and routines included in each directory are described brieﬂy. The source codes in thisarchive can be redistributable under the terms of GNU General Public License 5 , except for theones under the contrib/ directory.
src/ Main routines of the HISPUT program. Files and routines in this directory are describedin the next subsection.
lib/ The routines which can be used for generic purpose. Currently, there is only one subdirectory.
nakano/ This subdirectory contains: 1D/2D spline functions, 1-argument equation solverusing bisection method, and the FPU control interface for FORTRAN.
contrib/ This directory is the collection of routines developed by other groups. If you run
make with the Makeﬁle in this directory, you can build the library ﬁle including these.
mt/ Mersenne Twister routines.
mori/ A part of routines included in reference [87], e.g., the DE formula and the NC
method of 8th order.
slatec/ A selection of SLATEC routines utilized in our simulation program.
sparse/ SPARSE library routines applied the patch to be callable from FORTRAN.
examples/ Originally, these have been developed to check the respective routines in the HISPUT program. Some of them also are used to prepare the ﬁgures presented in this thesis.
45

http://www.debian.org/
http://www.gnu.org/copyleft/gpl.html

APPENDIX A. IMPLEMENTATION105

mt/ Test code for Mersenne Twister.
spline/ Test code for spline function routines.
bem/ Test code for BEM routine. The BEM data in ﬁg. 2.6 are the result of this routine.
bemcmp/ Sample program to output the MC data in ﬁg. 2.6.
eject/ Test codes for routines treating the ejection of the sputtered atom from the target.
coldist/ Test code to calculate the distribution of colliding gas speed vg for given energyof the sputtered atom (eq. 2.88).
colgas/ Test codes to obtain several parameters for the sputtered atom ﬂying throughthe gas environment with ﬁnite temperature, e.g., mean free path of the sputteredatom, and velocity of colliding gas atom. Used to prepare ﬁg. 2.4.
vrel/ Sample code to calculate the relative speed between the sputtered atom and thegas atom. Used to prepare ﬁg. 2.5.
scatter/ Sample code to calculate the scattering angle from the collision parameter andthe kinetic energy of the sputtered atom when the gas atom is stationary. Used toprepare ﬁgs. 2.2 and 3.13.
trap/ Test code for the routine accounting the particle trap on the chamber wall.
inc/ Test code for the routine to calculate the residual time of the sputtered atom in anarbitrary cell located inside of the chamber.
spreadf/ Routines and execution rules to obtain the time evolution of high energy atoms
spread into the free space ﬁlled with gas. Used to prepare the ﬁgures shown insection 3.2.
spreadb/ Routines and execution rules to obtain the evolution of sputtered atoms in an
actual sputtering environment. Used to prepare ﬁg. 3.7.
thesis-ﬁgs/ Rules and tools to obtain a series of physical parameters using HISPUT for different sputtering conditions (pressures, ejection energies, etc.). They include the script toprepare the input data, the batch program to execute the HISPUT with these input dataﬁles, and the script to summarize the calculation results.
00common/ Some blocks for the input data for the HISPUT, which can be used commonly: chamber boundaries, depth proﬁle of the target erosion track, atomic parameters, etc.
copper-tt/ Rules to obtain the dependence of the calculation results on the thermalization threshold energy (section 2.6.7).
copper/ Rules to evaluate the thickness proﬁle of sputtered Cu ﬁlms (section 3.3).
LaB6/ Rules to evaluate the pressure dependence of the ﬁlm composition in LaB6 sputtered ﬁlms (section 3.4).
density/ Rules to evaluate the atomic density proﬁle of Cu sputtered atoms (section 3.5).

APPENDIX A. IMPLEMENTATION

A.2.3

106

Source ﬁle and routines

In this subsection, ﬁles in src/ directory are mentioned. Functions and subroutines containedin them are also explained brieﬂy.
It is ruled that the common blocks are hidden from outside of the ﬁle. In other words,
they are referenced only from the routines in the same ﬁle, and the special interface routine areprepared to set the values to the variables in them.
For input/output of the program, device numbers 0, 5, 6, 9, 10 are used. The purpose of
these and the relationships with the Unix (Linux) devices are summarized in table A.1.

Number0
5
6910DevicePurpose/dev/stderr/dev/stdin/dev/stdout(user deﬁned)(user deﬁned)

To output debug information and errors
To input calculation parameters
To output miscellaneous information
To input cell data (see inc.F)
To output calculation data

Table A.1: I/O Devices used in HISPUT

The ﬁles with suﬃx .F are the sources that need to be processed by the C preprocessor
before the compilation by the FORTRAN compiler, while those with suﬃx .f are the sources
directly passed to the FORTRAN compiler.
main.F Main program.
program hisput Main program routine.
readin.F The routine to read the parameters passed to the program.
subroutine readin The routine to read the parameters passed to the program. The
name of the output ﬁle which will be allocated to the device number 10 is also readfrom this input ﬁle.
input.dat.samp Sample ﬁle to be read by subroutine readin.
aeject.F Calculate the ejection angle of the sputtered atom (see section 2.2.3).
subroutine aeject Calculate the ejection angle of the sputtered atom from the target
using Monte Carlo method.
eeject.f Calculate the ejection energy of the sputtered atom (see section 2.2.2).
subroutine eeject The routine to obtain the ejection energy of the sputtered atom.

APPENDIX A. IMPLEMENTATION

107

function fe Thompson’s formula (eq. 2.18).
subroutine seteej Set the parameters in common block ceej to be used by eeject andfe.
function dfe Derivative of the fe function. Used to obtain the maximum of fe.
peject.F Calculate the ejection position of the sputtered atom (see section 2.2.1).
subroutine peject The routine to deduce the ejection position of the sputtered atom
using eq. 2.4. Set common block epej.
subroutine setpej Initialization routine for peject. Construct the spline function torepresent the depth proﬁle of the erosion track.
function ddpej Derivative of the depth proﬁle function of the erosion track. Used toobtain the maximum of the depth proﬁle.
bem.F Calculate the diﬀusive stage of the particle transport by boundary element method(BEM) (see section 2.6).
subroutine bem Main routine of BEM.
subroutine bmform Generate the matrix in eq. 2.145.
subroutine elmoﬀ Calculate the oﬀ-diagonal component of the matrix in eq. 2.145.
subroutine elmdia Calculate the diagonal component of the matrix in eq. 2.145.
subroutine pivttq After the matrix is generated, swap the known/unknown components
to be handled by the LAPACK routines.
subroutine boundr Resume the positions of components swapped by pivttq after theequation is solved.
bemsrc.F Calculate the source term from the thermalized position of the sputtered atom usingeq. 2.146.
subroutine setbemsrc Initialize the array.
subroutine bemsrc Calculate the value inside of the summation of eq. 2.146. This
subroutine will be called every time when the sputtered atom is thermalized during
the Monte Carlo step.
subroutine readbemsrc Read out the sum of each node after the Monte Carlo calculation is over.
itq.f Deduce the value of fundamental solution on each node (see sections 2.6.3–2.6.5).
function ts θ̄ ∗ in eq. 2.134.
function tsi Integrand θ¯∗ x(i) , x R (x ) in eq. 2.142.
function qsi Integrand q¯∗ x(i) , x R (x ) in eq. 2.141.

APPENDIX A. IMPLEMENTATION

108

subroutine preptq Initialization routine for tsi and qsi. Set common block citq.
ellipt.f Interface functions for the SLATEC routines of complete elliptic functions.
function dellf Interface function for the complete elliptic function of 1st kind.
function dells Interface function for the complete elliptic function of 2nd kind.
cdir.F Determine the velocity of the sputtered atom after the collision with gas atom (seesection 2.5.3).
subroutine cdir Determine the velocity of the sputtered atom after the collision withgas atom.
subroutine setcdir Preparation routine for cdir. Set values in common block ccdir.
subroutine roty Rotation around the y axis (eq. 2.91).
subroutine rotz Rotation around the z axis (eq. 2.92).
subroutine relv Deduce the relative speed between the sputtered atom and the gas atom
using eq. 2.89.
coldist.F Calculate the probability function of the gas speed colliding to the sputtered atom
ﬂying through the thermal gas environment (see section 2.5.2).
subroutine setdist Initialization routine for distf function. Set values in common block
ccoldf. It takes the speed of the sputtered particle as an argument. Note that you
have to call this routine before calling distf when the the sputtered particle speed
changes.
function distf Calculate the right hand side of eq. 2.88.
colgas.F Routines to determine the mean free path of the sputtered particle traveling throughthe gas environment, and the speed of the colliding gas atom to it.
subroutine colgas Deduce the speed of colliding gas atom using the inverse function ofeq. 2.88. Actually it uses a spline function set by setcol.
subroutine frpath Deduce the mean free path of the sputtered particle from eq. 2.86.
Actually it uses a spline function set by setcol.
subroutine setcol Prepare the spline functions for colgas and frpath. Set the variables
in common block ccol and cmfp.
function bfun The function called from setcol. Used on solving the eq. 2.88 for vg by
the bisection method.
subroutine setbfun Set the pseudo random number w to common block cgaswp to use
it as a internal parameter of bfun.
scatter.F Calculate the scattering angle in the potential scattering framework (see section 2.3.2).

APPENDIX A. IMPLEMENTATION

109

subroutine scat Return the scattering angle χ0 obtained by eq. 2.40. Actually it uses a
spline function set by setsca.
subroutine setsca Prepare the spline function for scat. Set the variables in common
block csct.
function scint1 Integrand in eq. 2.40.
function scint2 Denominator of the scint1. Also used to determine the lower limit of
the integration by eq. 2.41.
sp.f Born-Mayer (BM) type scattering potential (eq. 2.46).
function sp Returns the value of BM potential. It takes the distance r as a parameter.
subroutine setsp Set the parameters A and b of the BM potential to common blockcsp.
trap.F Deciding routine of the trap of sputtered atoms on chamber wall nodes (see section 2.4).
subroutine trap Decide the trap of sputtered atoms on chamber wall.
subroutine settrap Initialization routine for trap. Store the positions of boundary
nodes in common block ctrap.
inc.F Calculate the residual time of sputtered atoms (see section 3.5).
subroutine inc Calculate the residual time of sputtered atoms.
subroutine setinc Initialization routine for inc. Read the input from device 9 and set
values to common block cincp.
subroutine repinc Output routine.
subroutine pinner Decide if the sputtered particle is inside of the cell of the speciﬁed
index.
subroutine iswap Swap two integers.
subroutine dswap Swap two double precisions.
inc.dat.samp Sample ﬁle to be read by subroutine inc. The device number 9 will be allocated to this ﬁle.
ecross.f Routine to determine whether the linear path of the sputtered atom in 3D space crosses
the line node in RZ space of the axisymmetric coordinate (see section 2.4).
subroutine ecross Routine to determine if the linear path of the sputtered atom in 3D
space crosses the line node in RZ space of the axisymmetric coordinate.
subroutine scross Initialization routine for ecross.
ccross.
subroutine rcross Sort routine called by ecross.

Set variables in common blockAPPENDIX A. IMPLEMENTATION

110

axarea.F Calculate the area in 3D of the nodes obtained by splitting the axisymmetric chamberboundary in RZ plane.
subroutine axarea For axisymmetric space, calculate the area of the line node in RZ
space.
The call tree of these functions and subroutines are visualized in ﬁgure A.1. This only
includes the routines in the ﬁles at src/ directory, and not the ones in the ﬁles at lib/ nor atcontrib/.
Program hisput can be built automatically by running the make utility with the Makefilein this src/ directory.
With this program, you can execute the simulation by:% ./hisput < input.dat.samp > output.dat
For detail of this execution procedure, please refer the input data generated by executing the
mkinputs.pl script and the batch ﬁle do.sh contained in subdirectories under thesis-figs/.

APPENDIX A. IMPLEMENTATION

Figure A.1: Call tree of HISPUT subroutines.

111Appendix B

Atomic Parameters
This table gives the atomic/elemental parameters used in the MC simulation of this study, which
include mass M , Born-Mayer potential (BMP) parameters A, b, and surface binding energy Us .
Atomic mass is taken from [100], BMP parameters are from [1], and surface binding energy isfrom [177]. Atomic mass is shown both in atomic unit and in kg.
Table B.1: Atomic mass and potential parameters
Z
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
1

Sym
He
Li
Be
B
C
N
O
F
Ne
Na
Mg
Al
Si
P
Si
Cl
Ar
K
Ca
Sc
Ti

M /a.u.
4.002602
6.941
9.012182
10.81112.0107
14.00674
15.9994
18.9984032
20.1797
22.989770
24.3050
26.981538
28.0855
30.973761
32.06635.452739.94839.0983
40.07844.955910
47.867

M /10−27 kg
6.646482
11.52614.96509
17.95219.9443
23.25875
26.5676
31.547612
33.5092
38.175437
40.359444.803928
46.6371
51.433175
53.24758.8706
66.33564.9243
66.55174.651096
79.485A/ (e2 /a0 )8.6047
16.10924.59935.60648.36762.84078.77196.267114.72134.56140.72157.85186.39205.67222.99235.64255.82277.95298.57319.62343.72

A/eV
234.13438.33669.34 1
968.841316.1
1709.92143.42619.43121.53661.43829.04295.15071.75596.36067.6
6411.86960.97563.08124.18696.99352.6b/ (a0 −1 )2.20779
2.120812.05904
2.02771
2.01592
2.00900
2.00474
2.00229
1.99954
1.99899
1.956941.946811.95888
1.951371.94026
1.92449
1.91901
1.91634
1.910441.905441.90402b/Å−1
4.17217
4.00780
3.89107
3.83187
3.809593.796513.78846
3.78383
3.778633.777593.698133.678993.70180
3.68760
3.666613.636813.626453.621413.610263.600813.59812

It is 699.34 in ref. [1], but it must be an error since it is inconsistent with the value of e2 /a0 .

112

Us /eV

3.32
5.77
7.37

3.39
4.63

4.85

APPENDIX B. ATOMIC PARAMETERS

113

Table B.1: Atomic mass and potential parameters
Z
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63

Sym
V
Cr
Mn
Fe
Co
Ni
Cu
Zn
Ga
Ge
As
Se
Br
Kr
Rb
Sr
Y
Zr
Nb
Mo
Tc
Ru
Rh
Pd
Ag
Cd
In
Sn
Sb
Te
In
Xe
Cs
Ba
La
Ce
Pr
Nd
Pm
Sm
Eu

M /a.u.
50.9415
51.996154.938049
55.84558.933200
58.6934
63.54665.39
69.72372.61
74.92160
78.96
79.90483.80
85.467887.62
88.90585
91.22492.9063895.94
99
101.07102.90550
106.42107.8682
112.411114.818118.710121.760
127.60126.90447
131.29132.90545
137.327138.9055
140.116140.90765
144.24145
150.36151.964

M /10−27 kg
84.5904
86.341691.226839
92.73397.860948
97.4628105.52108.58115.778120.57124.41033
131.12132.684
139.15141.9227145.50147.63174
151.481154.27478
159.31164.4
167.83170.87872
176.71
179.1195186.663190.660
197.123202.187
211.88210.7300
218.01220.69484
228.037230.6582
232.668
233.98282
239.52240.7
249.68252.342A/ (e2 /a0 )366.70389.77414.53438.46463.35487.72511.53539.78564.68590.33618.26645.35672.43703.70731.90763.87792.10819.91850.89881.65911.85941.72972.661005.21040.71071.11102.61140.21171.41199.41231.21265.21298.91336.41370.21403.31436.01471.51508.41544.31580.5A/eV9977.9
10606
11279
11931
12608
13271
13919
14687
15365
16063
16823
17560
18297
19148
19915
20785
21553
22310
23153
23990
24811
25624
26466
27352
28318
29145
30002
31025
31874
32636
33501
34426
35343
36363
37283
38184
39074
40040
41044
42020
43005

b/ (a0 −1 )1.90077
1.89794
1.89607
1.89300
1.89116
1.88818
1.88457
1.884241.88099
1.87903
1.87833
1.87688
1.87495
1.87520
1.87373
1.87395
1.872491.870421.869901.86926
1.86489
1.86723
1.866431.86599
1.86657
1.865211.86416
1.865101.86382
1.861361.86020
1.85980
1.85901
1.85956
1.85891
1.85796
1.85679
1.856641.85651
1.856211.85606b/Å−1
3.59198
3.58663
3.58310
3.57730
3.57382
3.56819
3.561373.56074
3.55460
3.550903.54958
3.54684
3.54319
3.54366
3.54088
3.54130
3.538543.53463
3.53365
3.532443.52418
3.52860
3.52709
3.52626
3.52735
3.52478
3.52280
3.52458
3.52216
3.517513.51532
3.51456
3.513073.514113.51288
3.51108
3.50887
3.508593.50834
3.507783.50749

Us /eV
5.31
4.10
2.92
4.28
4.39
4.44
3.49

3.85

6.25
7.57
6.82
6.74
5.75
3.89
2.95

3.14

APPENDIX B. ATOMIC PARAMETERS

114

Table B.1: Atomic mass and potential parameters
Z
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92

Sym
Gd
Tb
Dy
Ho
Er
Tm
Yb
Lu
Hf
Ta
W
Re
Os
Ir
Pt
Au
Hg
Tl
Pb
Bi
Po
At
Rn
Fr
Ra
Ac
Th
Pa
U

M /a.u.
157.25158.92534
162.50164.93032
167.26168.93421
173.04174.967
178.49180.9479
183.84186.207190.23
192.217195.078
196.96655
200.59204.3833
207.2
208.98038
210
210
222
223
226
227
232.0381
231.03588
238.0289

M /10−27 kg
261.12263.90192
269.84273.87343
277.74280.52205
287.34290.540296.39300.4712
305.27309.204
315.88319.184323.934
327.07087
333.09339.3867344.1
347.02032
348.7
348.7
368.6
370.3
375.2
376.9
385.3086
383.64437
395.2566

A/ (e2 /a0 )1615.91651.41688.7
1724.91765.11801.41842.3
1879.81914.1
1952.7
1991.22028.12067.92104.62144.22185.72230.02270.72311.82349.22395.92434.62476.52517.42623.02665.62716.92761.12806.6

A/eV
43969
44935
45950
46935
48028
49016
50129
51149
52083
53133
54181
55185
56268
57266
58344
59473
60678
61786
62904
63922
65192
66246
67386
68499
71372
72531
73927
75130
76368

b/ (a0 −1 )1.85571
1.855121.855131.85465
1.85493
1.854491.85489
1.85459
1.853441.85346
1.85300
1.852311.85240
1.85190
1.85191
1.85204
1.85254
1.85246
1.85258
1.851641.85224
1.85173
1.85163
1.85137
1.86368
1.86360
1.86470
1.86480
1.86519b/Å−13.50683
3.50572
3.50573
3.50483
3.505363.50452
3.505283.50471
3.502543.50258
3.501713.500413.50058
3.499633.499653.49989
3.50084
3.500693.50092
3.499143.50027
3.499313.499123.498633.52189
3.52174
3.52382
3.52401
3.52475

Us /eV
4.05

2.42

6.44
8.1
8.9
8.03
8.17
6.94
5.84
3.81

6.2
5.55

Acknowledgement (in Japanese)
本論文のご指導を快諾していただき、多くの有益なご指摘を賜りました東京大学教授 吉田豊信 博
士に心より感謝致します。また本論文の審査を引き受けていただき、同じく多くのご意見を頂きました東京大学教授 堀池靖浩 博士、東京大学教授 鳥海 明 博士、東京大学助教授 渡邊 聡 博士、
東京大学助教授 近藤高志 博士に深く感謝します。
成蹊大学教授 馬場 茂 博士には、本論文に関連する研究の遂行にあたり常に快適な研究環境を与えていただくとともに、様々なご指導をいただきました。
金沢工業大学教授 金原 粲 博士は、筆者の学生時代からの延長でもある本研究に関して、継続して有益なご指摘をくださいました。
京都大学教授 森 正武 博士は、本研究で開発したシミュレーションプログラムにおいて、ご著書の FORTRAN ルーチンを利用すること、およびそのルーチンを含むソースコードを公開・配布することを快諾してくださいました。
京都大学助教授 松本 眞 博士には、乱数発生ルーチン ”Mersenne Twister” の利用に関してご教示をいただきました。
成蹊大学教授 小島紀徳 博士ならびに岐阜大学助教授 上宮成之 博士には、LaB6 膜の ICP-Mass
測定を行っていただきました。
国立科学博物館 洞口俊博 博士 には、京都大学大型計算機センター ATMLINE データベースの利用に関するご教示をいただきました。
バキュームプロダクツ株式会社の 大橋哲二 氏には、本研究で利用した製膜装置の設計にあたって有益なご意見をいただき、また実際に作製を行っていただきました。成蹊大学機械工作室 佐藤
道憲 氏には、実験装置の作製の一部をご指導いただきました。
最後に、成蹊大学物理情報工学科 (前 計測数理工学科) 表面物性研究室の学生諸君には、様々な形で本研究の遂行を助けていただきました。特に森 意佐央君にはガス熱運動を取り扱うルーチンの開発において、大貫徳之君には Cu プラズマの分光測定において、それぞれ協力していただ
きました。ここに記して感謝の意を表します。

115

Bibliography[1] A. A. Abrahamson. Born-Mayer-type interatomic potential for neutral ground-state atoms
with Z = 2 to Z = 105. Phys. Rev., 178:76 – 79, 1969.
[2] J. B. Almeida. Design of magnetrons for dc sputtering. Vacuum, 39:717 – 721, 1989.
[3] J. N. Avaritsiotis and C. D. Tsiogas. A reactive sputtering process model for symmetricalplanar diode system. Thin Solid Films, 209:17 – 25, 1992.
[4] Z Bai, C. Bischof, J. Demmel, and J. Dongarra. Lapack Users’ Guide. the Society ofIndustrial and Applied Mathematics, Philadelphia, 2nd package edition, 1995. (editor: E.
Anderson).
[5] L. T. Ball, I. S. Falconer, D. R. McKenzie, and J. M. Smelt. An interferometric investigation of the thermalization of copper atoms in a magnetron sputtering discharge. J. Appl.
Phys., 59(3):720 – 725, 1986.
[6] Albert-László Barabási and H. Eugene Stanley. Fractal concepts in surface growth. Cambridge Univ. Press, Cambridge, 1995.
[7] S. Berg and I. B. Katardjiev. Preferencital sputtering eﬀects in thin ﬁlm processing. J.
Vac. Sci. Technol. A, 17(4):1916 – 1925, 1999.
[8] G. A. Bird. Molecular gas dynamics and the direct simulation of gas flows. Oxford
University Press, New York, 1994.
[9] A. J. Blair, G. Metzger, and C. B. Fleddermann. Optical absorption spectroscopic diagnostics during sputter deposition of Y-Ba-Cu-oxide. J. Appl. Phys, 72(10):4792 – 4797,
1992.
[10] A. Bogaerts, R. Gijbels, and W. J. Goedheer. Hybrid Monter Carlo-ﬂuid model of a direct
current glow discharge. J. Appl. Phys., 78(4):2233 – 2241, 1995.
[11] A. Bogaerts, J. Naylor, M. Hatcher, W. J. Jones, and R. Mason. Inﬂuence of sticking
coeﬃcients on the behavior of sputtered atoms in an argon glow discharge: Modeling andcomparison with experiment. J. Vac. Sci. Technol. A, 16(4):2400 – 2410, 1998.
[12] Annemie Bogaerts, Mark van Straaten, and Renaat Gijbels. Description of the thermalization process of the sputtered atoms in a glow discharge using a three dimensional Monte
Carlo method. J. Appl. Phys., 77(5):1868 – 1874, 1995.
116

BIBLIOGRAPHY

117[13] R. Mark Bradley and James M. E. Harper. Theory of ripple topography induced by ion
bombardment. J. Vac. Sci. Technol. A, 6(4):2390 – 2395, 1988.
[14] L. C. Wrobel C. A. Brebbia, J. C. F. Telles. Boundary Element Techiniques Theory andApplications in Engineering. Springer-Verlag, Berlin, 1984.
[15] C. C. Chang. Energy dependence of angular distributions of sputtered particles from theAg{111} surface. Phys. Rev. B:, 48(20):12399 – 12405, 1993.
[16] J. W. Coburn and E. Kay. Positive-ion bombardment of substrates in rf diode glow
discharge sputtering. J. Appl. Phys., 43:4965 – 4971, 1972.
[17] Radlfo Cuerno and Albert-László Barabási. Dynamic scaling of ion-sputterd surfaces.
Phys. Rev. Lett., 74(23):4746 – 4749, 1995.
[18] Rodolfo Cuerno, Hernan A. Makse, Silvina Tomassone, Stephen T. Harrington, and H. Eugene Stanley. Stochastic model for surface erosion via ion sputtering: Dynamical evolution
from ripple morphology to rough morphology. Phys. Rev. Lett., 75(24):4464 – 4467, 1995.
[19] J. J. Cuomo and S. M. Rossnagel. Hollow-cathode-enhanced magnetron sputtering. J.
Vac. Sci. Technol. A, 4(3):393 – 396, 1986.
[20] Carl de Boor. A Practical Guide to Splines. Springer-Verlag, New York, 1978.
[21] M. P. Desjarlais. The eﬀect of charge exchange processes on ion diode impedance. J. Appl.
Phys., 66:4696 – 4701, 1989.
[22] S. K. Dew. Theoretical and practical aspects of collimated sputtering. J. Appl. Phys.,
76(8):4857 – 4862, 1994.
[23] S. K. Dew, T. Smy, R. N. Tait, and M. J. Brett. Modeling bias sputter planarization of
metal ﬁlms using a ballistic deposition simulation. J. Vac. Sci. Technol. A, 9(3):519 – 523,
1991.
[24] Steven Dew, Tom Smy, and Michael Brett. Step coverage, uniformity and composition
studies using intergrated vapour transport and ﬁlim deposition models. Jpn. J. Appl.
Phys., 32(2):1140 – 1145, 1994.
[25] M. F. Dony, J. P. Dauchot, M. Wautelet, M. Hecq, and A. Richard. Diagnostics by opticalabsorption of sputtered atom density in magnetron discharges. J. Vac. Sci. Technol. A,
18(3):809 – 813, 2000.
[26] M. F. Dony, A. Ricard, J. P. Dauchot, M. Hecq, and M. Wautelet. Optical diagnostics ofd.c. and r.f. argon magnetron discharges. Surf. Coatings Technol., 74-75:479 – 484, 1995.
[27] M. F. Dony, A. Ricard, M. Wautelet, J. P. Dauchot, and M. Hecq. Comparison of direct
current and radio frequency argon discharges by optical emission and absorption spectroscopy. J. Vac. Sci. Technol. A, 15:1890 – 1896, 1997.

BIBLIOGRAPHY

118[28] S. Dosa, S. Ghosal, R. L. Kosut, J. L. Ebert, T. E. Abrahamson, A. Kozak, D. W. Zou,
X. Zhou, J. F. Groves, and H. N. G. Wadley. Reactor-scale models for rf diode sputteringof metal thin ﬁlms. J. Vac. Sci. Technol. A, 17(4):1926 – 1933, 1999.
[29] Tilo P. Drüsedau, Mirko Löhmann, and Bernd Garke. Decay length of the pressure dependent deposition rate for magnatron sputtering. J. Vac. Sci. Technol. A, 16(4):2728 –
2732, 1998.
[30] C. Eisenmenger-Sittner, R. Beyerknecht, and A. Bergauer. Angular distribution of sputtered neutrals in a post magnetron geometry: Measurement and Monte Carlo simulation.
J. Vac. Sci. Technol. A, 13(5):2435 – 2443, 1995.
[31] R. Elsing. The practical use and application of Monte-Carlo studies in physical vapour
deposition technology. Surf. Coat. Technol., 49:132 – 138, 1991.
[32] C. B. Eom, J. Z. Sun, K. Yamamoto, A. F. Marshall, K. E. Luther, T. H. Geballe, andS. S. Landerman. In situ grown YBa2 Cu3 O7−δ thin ﬁlms from single-target magnetron
sputtering. Appl. Phys. Lett., 55:595 – 597, 1989.
[33] C. C. Fang, F. Jones, R. R. Kola, G. K. Celler, and V. Prasad. Stress and microstructureof sputter-deposited thin ﬁlms: Molecular dynamics simulations and experiment. J. Vac.
Sci. Technol. B, 11(6):2947 – 2952, 1993.
[34] C. C. Fang, V. Prasad, and F. Jones. Molecular dynamics modeling of microstructure andstresses in sputter-deposited thin ﬁlms. J. Vac. Sci. Technol. A, 11(5):2778 – 2789, 1993.
[35] L. J. Friedrich, S. K. Dew, M. Brett, and T. Smy. Thin ﬁlm microstructure modellingthrough line-segment simulation. Thin Solid Films, 266:83 – 88, 1995.
[36] L. J. Friedrich, S. K. Dew, M. J. Brett, and T. Smy. Integrating system and feature scalemodels to study copper reﬂow. J. Vac. Sci. Technol. B, 17(1):186 – 193, 1999.
[37] L. J. Friedrich, D. S. Gardner, S. K. Dew, M. J. Brett, and T. Smy. Study of the copper
reﬂow process using the GROFILMS simulator. J. Vac. Sci. Technol. B, 15(5):1780 –
1787, 1997.
[38] J. R. Gavaler, J. Talvacchio, T. T. Braggins, M. G. Forrester, and J. Greggi. Criticalparamters in the sinble-target sputtering of YBa2 Cu3 O7 . J. Appl. Phys., 70:4383 – 4391,
1991.
[39] Michael J. Grapperhaus, Zoran Krivokapic, and Mark J. Kushner. Design issues in ionized
metal physical vapor deposition of copper. J. Appl. Phys., 83(1):35 – 43, 1998.
[40] S. Hamaguchi and S. M. Rossnagel. Simulations of trench-ﬁlling proﬁles under ionizedmagnetron sputtermetal deposition. J. Vac. Sci. Technol. B, 13(2):183 – 191, 1995.
[41] G. L. Hansen and R. G. Ahonen. The eﬀect of gas scattering on the deposition proﬁle of
optical thin ﬁlms. J. Vac. Sci. Technol. A, 5(5):2898 – 2901, 1987.

BIBLIOGRAPHY

119[42] M. Hasegawa and K. Nanbu. Numerical analysis of ﬁlm growth and step coverage in thecollimation sputtering method. Vacuum, 48(10):825 – 831, 1997.
[43] T. Heberlein, G. Krautheim, and W. Wuttke. The sputter deposition process : a MonteCarlo study. Vacuum, 42:47 – 51, 1991.
[44] D. Henderson, M. H. Brodsky, and P. Chaudhari. Simulation of structural anisotropy andvoid formation in amorphous thin ﬁlms. Appl. Phys. Lett., 25(11):641 – 643, 1974.
[45] D. W. Hoﬀman. A sputteirng wind. J. Vac. Sci. Technol. A, 3(3):561 – 566, 1985.
[46] D. W. Hoﬀman. Perspective on stresses in magnetron-sputtered thin ﬁlms. J. Vac. Sci.
Technol. A, 12(4):953 – 961, 1994.
[47] Jeﬀrey A. Hopwood, editor. Ionized Physical Vapor Deposition. Academic Press, San
Diego, 2000.
[48] T. Horaguchi and T. Hirata. ATMLINE: relational database for atomic spectrum lines.
Kyoto University Data Processing Center Bulletin, 22:129 – 141, 1989. (in Japanese).
[49] Yoichi Hoshi and Eisuke Suzuki. Investigation of sputter deposotion process of iron ﬁlms
by computer simulation. IONICS, 21(1):37 – 44, 1995. (in Japanese).
[50] R. P. Howson, H. A. J’afer, and A. G. Spencer. Substrate eﬀects from an unbalanced
magnetron. Thin Solid Films, 193/194:127 – 137, 1990.
[51] T. C. Huang, G. Lim, F. Parmigiani, and E. Kay. Eﬀect of ion bombardment during deposition on the x-ray microstructure of thin silver ﬁlms. J. Vac. Sci. Technol. A, 3(6):2161– 2166, 1985.
[52] Shunji Ido, Mieko Kashiwagi, and Mikihiko Takahashi. Computational studies of plasma
generation and control in a magnetron sputtering system. Jpn. J. Appl. Phys., 38(7B):4450– 4454, 1999.
[53] Shunji Ido and Kosuke Nakamura. Computational studies on plasma generation anderosion in rectangular magnetron sputtering system. Jpn. J. Appl. Phys., 35(4A):2302 –
2305, 1996.
[54] Shunji Ido, Takeshi Suzuki, and Mieko Kashiwagi. Computational studies on the erosion
process in a magnetron sputtering system with a ferromagnetic target. Jpn. J. Appl. Phys.,
37(3A):965 – 969, 1998.
[55] Masanobu Ishida, Yukio Yamaguchi, and Yasunori Yamamura. Simulation of the MO ﬁlm
deposition in the magnetron sputtering discharge. Thin Solid Films,, 334:225 – 229, 1998.
[56] J. Jeans. The dynamical theory of gases. Dover publications, Dover, 4th edition, 1954.
[57] Robert W. Johnson and A. H. Daane. The langhanum-boron system. J. Phys. Chem.,
65:909 – 915, 1961.

BIBLIOGRAPHY

120[58] S. Kadlec and J. Musil. Optimized magnetic ﬁeld shape for low pressure magnetron
sputtering. J. Vac. Sci. Technol. A, 13(2):389 – 393, 1995.
[59] T. Kajiwara, T. Urakabe, K. Sano, K. Fukuyama, K. Watanabe, S. Baba, T. Nakano, andA. Kinbara. Mechanical and electrical properties of rf sputtered LaB6 thin ﬁlms on glass
substrates. Vacuum, 41:1224 – 1228, 1990.
[60] Manfred Kaminsky. Atomic and ionic impact phenomena on metal surfaces. SpringerVerlag, Berlin, 1965.
[61] Suk Tai Kang, Ryuichi Shimizu, and Tsuyoshi Okutani. Sputtering of Si with keV Ar+
ions. I. measurement and Monte Carlo calculations of sputtering yield. Jpn. J. Appl. Phys.,
18:1717 – 1725, 1979.
[62] E. Kay, F. Parmigiani, and W. Parrish. Microstructure of sputtered metal ﬁlms grown in
high- and low-pressure discharges. J. Vac. Sci. Technol. A, 6(6):3074 – 3081, 1988.
[63] J. H. Keller and R. G. Simmons. Sputtering process model of deposition rate. IBM J.
Res. Develop., 23(1):24 – 32, 1979.
[64] A. Kersch, W. Morokoﬀ, and Chr. Werner. Selfconsistent simulation of sputter deposition
with the Monte Carlo method. J. Appl. Phys., 75(4):2278 – 2285, 1994.
[65] A. Kinbara, T. Nakano, A. Kobayashi, S. Baba, and T. Kajiwara. LaBx thin ﬁlms preparedby magnetoron sputtering. Appl. Surf. Sci., 70/71:742 – 745, 1993.
[66] Donald E. Knuth. The Art of Computer Programming, Volume 2: Seminumerical Algorisms. Addisoin Wesley, Massachusetts, 3rd edition, 1997.
[67] P. Lecoeur, B. Mercey, and H. Murray. Inﬂuence of angular distribution on the deposition rate of species sputtered from a multicomponent target in diﬀernt conﬁgurations:Appliocations to mixed valene copper oxides. J. Vac. Sci. Technol. A, 13(4):2221 – 2227,
1995.
[68] Francisco A. Leon. Numerical modeling of glass ﬂow and spin-on planarization. IEEE
Trans. CAD, 7(2):168 – 173, 1988.
[69] Michael A. Lieberman and Allan J. Lichtenberg. Principles of plasma discharges andmaterials processing. Jon Wiley & Sons, Inc., New York, 1994.
[70] Z. Lin and T. S. Cale. Flux distributions and deposition proﬁles from hexagonal collimatorsdduring sputter deposition. J. Vac. Sci. Technol. A, 13(4):2183 – 2188, 1995.
[71] Mike Loukides. UNIX for FORTRAN Programmers. O’Reilly & Associates, Sebastopol,
1990.
[72] John E. Mahan and Andre Vantomme. A simpliﬁed collisional model of sputtering in thelinear cascade regime. J. Vac. Sci. Technol. A, 15(4):1976 – 1989, 1997.

BIBLIOGRAPHY

121[73] John E. Mahan and Andrè Vantomme. Trends in sputter yield data in the ﬁlm deposition
regime. Phys. Rev. B:, 61(12):8516 – 8525, 2000.
[74] Maxim A. Makeev and Albert-László Barabási. Ion-induced eﬀective surface diﬀusion in
ion sputtering. Appl. Phys. Lett., 71(19):2800 – 2802, 1997.
[75] A. Malaurie and A. Bessaudou. Numerical simulation of the charactersistics of the diﬀerent
metallic species falling on the growing ﬁlm in d.c. magnetron sputtering. Thin Solid Films,
286:305 – 316, 1996.
[76] A. Manenschijn and W. J. Goedheer. Angular ion and nertral energy distribution in acollisional rf sheath. J. Appl. Phys., 69(5):2923 – 2930, 1991.
[77] A. Manenschijn, G. C. A. M. Janssen, E. van der Drift, and S. Radelaar. Measurement ofion impact energy and ion ﬂux at the rf electrode of a parallel plate reactive ion etcher.
J. Appl. Phys., 69(3):1253 – 1262, 1991.
[78] Rod S. Mason, John C. Naylor, Michael G. Hatcher, and W. Jeremy Jones. Production,
thermalization, and transport of sputtered lithium atoms across a glow discharge plasma,studied by concentration modulated absorption spectroscopy. J. Appl. Phys., 84(3):1179– 1185, 1988.
[79] D. M. Mattox. Particle bombardment eﬀects on thin-ﬁlm depositoin: A review. J. Vac.
Sci. Technol. A, 7(3):1105 – 1114, 1989.
[80] T. M. Mayer, E. Chason, and A. J. Howard. Roughening instability and ion-induced
viscous relaxation of SiO2 surfaces. J. Appl. Phys., 76(3):1633 – 1643, 1994.
[81] A. A. Mayo, S. Hamaguchi, J. H. Joo, and S. M. Rossnagel. Across-wafer nonuniformityof long throw sputter deposition. J. Vac. Sci. Technol. B, 15(5):1788 – 1793, 1997.
[82] T. Mehdi, P. B. Legrand, J. P. Dauchot, M. Wautelet, and M. Hecq. Optical emission
diagnostics of an rf magnetron sputtering discharge. Spectrochim. Acta, 48B:1023 – 1033,
1993.
[83] Kevin Meyer, Ivan K. Schuller, and Charles M. Falco. Thermalization of sputtered atoms.
J. Appl. Phys., 52(9):5803 – 5805, 1981.
[84] S. Miyake, N. Shimura, and T. Makabe. Diagnostics of direct-current-magnetron discharges by the emission-selected computer-tomography technique. J. Vac. Sci. Technol.
A, 10:1135 – 1139, 1992.
[85] Dae Won Moon and Kyung Joong Kim. Surface topography development on ion-beamsputtered surface: Role of surface inhomogeneity induced by ion-beam bombardment. J.
Vac. Sci. Technol. A, 14:2744 – 2756, 1996.
[86] Charlotte E. Moore. Atomic Energy Levels, volume 2 of NSRBS-NDS 35. National Bureauof Standards, 1971.

BIBLIOGRAPHY

122[87] Masatake Mori. FORTRAN 77 Suuchi Keisan Programing. Iwanami Shoten, Tokyo, 1986.
(in Japanese).
[88] Shigeichi Moriguchi, Kanehisa Udagawa, and Shin Hitotsumatsu.
Koushiki. Iwanami Shoten, Tokyo, 1960. (in Japanese).

Iwanami Suugaku[89] Nobuhiro Motegi, Yuzou Kashimoto, Koji Nagatani, Seiichi Takahashi, Tomoyasu Kondo,
Yasushi Mizusawa, and Izumi Nakayama. Long-throw low-pressure sputtering technology
for very large-scale integrated devices. J. Vac. Sci. Technol. B, 13:1906 – 1909, 1995.
[90] T. Motohiro. Applications of Monte Carlo simulation in the analysis of sputter-deposition
process. J. Vac. Sci. Technol. A, 4(2):189 – 195, 1986.
[91] T. Motohiro and Y. Taga. Monte Carlo simulation of thermalization process of sputterd
particles. Surf. Sci., 134:L494 – L499, 1983.
[92] T. Motohiro and Y. Taga. Monte Carlo simmulation of the practice transport process in
sputter deposition. Thin Solid Films, 112:161 – 173, 1984.
[93] T. Motohiro and Y. Taga. On the elemental dependence of arrivval rate of the sputtered particles in the Monte Carlo simulation of the particle transport process in sputter
deposition. Thin Solid Films, 172:L105 – L107, 1989.
[94] S. J. Mroczkowski. Electron emissinon characteristics of sputtered lanthanum hexaboride.
J. Vac. Sci. Technol., A9:586 – 590, 1991.
[95] Karl-Heinz Müller. Stress and microstructure of sputter-deposited thin ﬁlms: Molecular
dynamics investigations. J. Appl. Phys., 62(5):1796 – 1799, 1987.
[96] G. K. Muralidhar, J. Musil, and S. Kadlec. An unbalanced magnetron sputtering device
for low and medium pressures. Rev. Sci. Instrum., 66(10):4961 – 4966, 1995.
[97] M. Muroi, T. Matsui, Y. Okamura, K. Tsuda, M. Nagano, and K. Mukae. Low temperaturesynthesis of YBa2 Cu3 O7−x thin ﬁlms. J. Mater. Res., 4:781 – 786, 1989.
[98] A. M. Myers, J. R. Doyle, J. R. Abelson, and D. N. Ruzic. Monte Carlo simulations ofmagnetron sputtering particle transport. J. Vac. Sci. Technol. A, 9(3):614 – 618, 1991.
[99] K. Nanbu, S. Uchida, and H. Yoshida. Computer simulation of growth of thin ﬁlms
fabricated by the sputtering method: comparison with experiment. Thin Solid Films,
228:330 – 332, 1993.
[100] National Astronomical Observatory, editor. Rika Nenpyo (Chronological Scientific Tables
1999). Maruzen, Tokyo, 1998.
[101] John F. O’Hanlon. A User’s Guide to Vacuum Technology. Jon Wiley & Sons, Inc., New
York, 1989.
[102] T. Ohta and H. Yamada. A sputter equipment simulation system for VLSI device. Vacuum,
51(4):479 – 484, 1998.

BIBLIOGRAPHY

123[103] T. Okutani, M. Shikata, S. Ichimura, and R. Shimizu. Auguiar distribution of Si atoms
sputtered by keV Ar+ ions. J. Appl. Phys., 51:2884 – 2887, 1980.
[104] W. Z. Park, T. Eguchi, C. Honda, K. Muraoka, Y. Yamagata, B. W. James, M. Maeda,
and M. Akazaki. Investigation of the thermalizaiton of sputtered atoms in a magnetron
discharge using laser-induced ﬂuorescence. Appl. Phys. Lett., 58(3):2564 – 2566, 1991.
[105] Robert Parsons. Sputter deposition processes. In John L. Vossen and Werner Kern, editors,
Thin Film Processes II, chapter II-4. Academic Press, San Diego, 1991.
[106] H. Patterson and D. H. Tomlin. Experiments by radioactive tracer methods on sputteringby rare-gas ions. Proc. Roy. Soc. London Ser.A, 265:474 – 488, 1962.
[107] A. M. C. Pérez-Martin and J. J. Jiménez-Rodrı́guez. Angular distribution and escapeprobability of sputtered atoms. Vacuum, 42:537 – 542, 1991.
[108] A. J. Perry, D. Vender, and R. W. Boswell. The application of the helicon source to plasma
sputtering. J. Vac. Sci. Technol. B, 9(2):310 – 317, 1991.
[109] I. Petrov, I. Ivanov, V. Orlinov, and J. E. Sundgren. Comparison of magnetron sputter
deposition conditions in neon, argon, krypton, and xenon discharges. J. Vac. Sci. Technol.
A, 11(5):2733–2741, 1993.
[110] Z. J. Radzimski, O. E. Hankins, W. P. Posadowski J. J. Cuomo, and S. Shingubara. Opticalemission spectroscopy of high density metal plasma formed during magnetron sputtering.
J. Vac. Sci. Techonol. B, 15:202 – 208, 1997.
[111] G. Mohan Rao and S. Mohan. Optimization studies on magnetic ﬁeld geometry for planar
magnetron sputtering targets. J. Vac. Sci. Technol. A, 9(6):3100 – 3104, 1991.
[112] R. L. Rhoades and S. M. Gorbatkin. Characterization of Ar/Cu electron-cycrotronresonance plasmas using optical emissio nspectroscopy. J. Appl. Phys., 80:2605 – 2613,
1996.
[113] R. S. Robinson. Energetic binary collisions in rare gas plasmas. J. Vac. Sci. Technol.,16(2):185 – 188, 1979.
[114] S. M. Rossnagael. Directional and preferential sputtering-based physical vapor deposition.
Thin Solid Films, 263:1 – 12, 1995.
[115] S. M. Rossnagel. Gas density reduction eﬀects in magnetrons. J. Vac. Sci. Technol. A,
6(1):19 – 24, 1988.
[116] S. M. Rossnagel. Energetic particle bombardment of ﬁlms during magnetron sputtering.
J. Vac. Sci. Technol. A, 7(3):1025 – 1029, 1989.
[117] S. M. Rossnagel. Directional and ionized physical vapor deposition for microelectronics
applications. J. Vac. Sci. Technol. B, 16:2585 – 2608, 1998.

BIBLIOGRAPHY

124[118] S. M. Rossnagel and J. J. Cuomo. Film modiﬁcation by low energy ion bombardment
during deposition. Thin Solid Films, 171:143 – 156, 1989.
[119] S. M. Rossnagel and J. Hopwood. Magnetron sputter deposition with high levels of metal
ionization. Appl. Phys. Lett., 63(24):3285 – 3287, 1993.
[120] S. M. Rossnagel and J. Hopwood. Metal ion deposition from ionized magnetron discgarge.
J. Vac. Sci. Technol. B, 12:449 – 453, 1994.
[121] S. M. Rossnagel, D. Mikalsen, H. Kinoshita, and J. J. Cuomo. Collimated magnetron
sputter deposition. J. Vac. Sci. Technol. A, 9(2):261 – 265, 1991.
[122] S. M. Rossnagel and K. L. Saeger. Optical emission in magnetrons: Nonlinear aspects. J.
Vac. Sci. Technol. A, 7(3):968 – 971, 1989.
[123] S. M. Rossnagel, I. Yang, and J. J. Cuomo. Compositional changes during magnetron
sputtering of alloys. Thin Solid Films, 199:59 – 69, 1991.
[124] D. N. Ruzic and H. K. Chiu. Modeling of particle-surface reﬂections including surface
roughness characterized by fractal geometry. J. Nucl. Mater., 162-164:904 – 909, 1989.
[125] David N. Ruzic. The eﬀects of surface roughness characterized by fractal geometry on
sputtering. Nucl. Instr. and Meth., B47(3):118 – 125, 1990.
[126] J. G. Ryan and S. Roberts. The formation and characterization of rare earth boride ﬁlms.
Thin Solid Films, 135:9 – 19, 1986.
[127] K. Sakuta, M. Iyori, Y. Katayama, and T. Kobayashi. Eﬀect of discharge gas pressure
on YBaCuO epitaxial ﬁlm formation by reactive RF magnetron sputtering. Jpn. J. Appl.
Phys., 29:L611 – L613, 1990.
[128] S. Scaglione, L. Caneve, and F. Sarto. Sputtering yield of optical materials: Sigmund’smodel and experimental results. J. Vac. Sci. Technol. A, 12(4):1523 – 1527, 1994.
[129] Ivan K. Schuller and Charles M. Falco. Structure and physical properties of sputteredmetallic superlattices. Surf. Sci., 113:443 – 453, 1982.
[130] Vladimir V. Serikov and Kenichi Nanbu. Monte Carlo numerical analysis of target erosion
and ﬁlm growth in a three-dimensional sputtering chamber. J. Vac. Sci. Technol. A,
14:3108 – 3123, 1996.
[131] Vladimir V. Serikov and Kenichi Nanbu. The analysis of background gas heating in direct
current sputtering discharges via particle simulation. J. Appl. Phys., 82(12):5948 – 5957,
1997.
[132] T. E. Sheridan, M. J. Goeckner, and J. Goree. Model of energetic electron transport inmagnetron discharges. J. Vac. Sci. Technol. A, 8:30 – 37, 1990.

BIBLIOGRAPHY

125[133] T. E. Sheridan, M. J. Goeckner, and J. Goree. Electron velocity distribution functions in a
sputtering magnetron discharge for the E×B direction. J. Vac. Sci. Technol. A, 16(4):2173– 2176, 1998.
[134] Eiji Shidoji, Hiroto Ohtake, Nobuhiko Nakano, and Toshiaki Makabe. Two-dimensional
self-consistent simulation of a dc magnetron discharge. Jpn. J. Appl. Phys., 38(4A):2131– 2136, 1999.
[135] Peter Sigmund. Theory of sputtering. I. sputtering yield of amorphous and polycrystalline
targets. Phys. Rev., 184:383 – 416, 1969.
[136] Peter Sigmund. A mechanism of surface micro-roughening by ion bombardment. J. Mater.
Sci., 8:1545 – 1553, 1973.
[137] T. Smy, L. Tan, S. S. Winterton, S. K. Dew, and M. J. Brett. Simulation of sputter
deposition at high pressures. J. Vac. Sci. Techinol. A, 15:2847 – 2853, 1997.
[138] Tom Smy, Kenneth L. Westra, and Michael J. Brett. Simulation of density variation andstep coverage for a variety of via/contact geometries using SIMBAD. IEEE Trans. Electron
Devices, 37(3):591 – 598, 1990.
[139] R. E. Somekh. The thermalization of energetic atoms during the sputtering process.
J. Vac. Sci. Technol. A, 2:1285 – 1291, 1984.
[140] J. Stache. Hybrid modeling of deposition proﬁles in magnetron sputtering systems. J.
Vac. Sci. Technol. A, 12(5):2867 – 2872, 1994.
[141] Keisuke Sugano, Kazuyoshi Yoshimura, and Fumio Takayama. C ni yoru spline kansuu.
Tokyo Denki University Press, Tokyo, 1993. (in Japanese).
[142] S. Swann. Film thickness distribution in magnetron sputtering. Vacuum, 38:791 – 794,
1988.
[143] R. N. Tait, S. K. Dew, T. Smy, and M. J. Brett. Density variation of tungsten ﬁlims
sputtered over topography. J. Appl. Phys., 70(8):4295 – 4300, 1991.
[144] R. N. Tait, T. Smy, and M. J. Brett. Simulation and measurement of density variation in
Mo ﬁlms sputter deposited over oxide steps. J. Vac. Sci. Technol. A, 8(3):1593 – 1596,
1990.
[145] N. Terada, H. Ihara, M. Jo, M. Hirabayashi, Y. Kimura, K. Matsutani, K. Hirata, E. Ohno,
R. Sugise, and F. Kawashima. Sputter synthesis of Ba2 YCu3 Oy as-deposited superconducting thin ﬁlms from stoichiometric target—A mechanism of compositional deviationand its control. Jpn. J. Appl. Phys., 27:L639 – L642, 1988.
[146] M. W. Thompson. II. the energy spectrum of ejected atoms during the high energy
sputteirng of gold. Philos. Mag., 18:377 – 414, 1968.

BIBLIOGRAPHY

126[147] J. A. Thornton. The microstructure of sputter-deposited coatings. J. Vac. Sci. Technol.
A, 4(6):3059 – 3065, 1986.
[148] John A. Thornton. Inﬂuence of apparatus geometry and deposition conditions on the
structure and topography of thick sputtered coatings. J. Vac. Sci. Technol., 11(4):666 –
670, 1974.
[149] H. Tsuge and S. Esho. Angular distribution of sputtered atoms from polycrystalline metaltargets. J. Appl. Phys., 52(7):4391 – 4395, 1981.
[150] R. T. C. Tsui. Calculation of ion bombarding energy and its distribution in rf sputtering.
Phys. Rev.,, 168(1):107 – 113, 1968.
[151] Mutumi Tuda, Kouichi Ono, and Akimasa Yuuki. Transport and deposition processesof sputtered particles in rf-microwave hybrid sputtering discharges. Jpn. J. Appl. Phys.,
33(7B):4473 – 4477, 1994.
[152] G. M. Turner. Monte Carlo calculations of gas rarefaction in a magnetron sputteringdischarge. J. Vac. Sci. Technol. A, 13(4):2161 – 2169, 1995.
[153] G. M. Turner, I. S. Falconer, B. W. James, and D. R. McKenzie. Monte Carlo calculation of the thermalization of atoms sputtered from the cathode of a sputtering discharge.
J. Appl. Phys., 65:3671 – 3679, 1989.
[154] G. M. Turner, I. S. Falconer, B. W. James, and D. R. McKenzie. Monte Carlo calculationsof the properties of sputtered atoms at a substrate surface in a magnetron discharge. J.
Vac. Sci. Technol. A, 10(3):455 – 461, 1992.
[155] G. M. Turner, A. Sikorski, D. R. McKenzie, G. B. Smith, K. Ng, and D. J. H. Cockayne. Spatial variations in the stoichiometry of sputtered YBaCuO thin ﬁlms: theory andexperiment. Physica C, 170:473 – 480, 1990.
[156] H. M. Urbassek and D. Sibold. Sputterd atom transport in high-current gas discharges:
A self-consistent computer simulation study. J. Vac. Sci. Technol. A, 11:676 – 681, 1993.
[157] M. A. Vidal and R. Asomoza. Monte Carlo simulation of the transport process in thegrowth of a-Si:H prepared by cathodic reactive sputtering. J. Appl. Phys., 67(1):477 –
482, 1990.
[158] Robert K. Waits. Planar magnetron sputtering. In John L. Vossen and Werner Kern,
editors, Thin Film Processes, chapter II-4. Academic Press, San Diego, 1978.
[159] W. Wang, J. Foster, A. E. Wendt, T̃. Onuoha J. H. Booske, P. W. Sandstorm, H. Liu,
S. S. Gearhart and N. Hershkowitz. Magnetic-ﬁeld-enhanced rf argon plasma for ionized
sputtering of copper. Appl. Phys. Lett., 71:1622 – 1624, 1997.
[160] X. Wang, H. Masumoto, Y. Someno, and T. Hirai. Optical properties of titania/silica multilayer ﬁlters prepared by helicon plasma sputtering. J. Vac. Sci. Technol. A, 16(5):2926– 2930, 1998.

BIBLIOGRAPHY

127[161] Kiyotaka Wasa and Shigeru Hayakawa. Handbook of sputter deposition technology. Noyes
Publications, New Jersey, 1992.
[162] Bernhard Welz and Michael Sperling. Atomic absorption spectrometry. Wiley-VCH, Weinheim, 3rd edition, 1999.
[163] W. D. Westwood. Calculation of deposition rates in diode sputtering systems. J. Vac.
Sci. Technol., 15(1):1 – 9, 1978.
[164] H. Windischmann. Intrinsic stress in sputtered thin ﬁlms. J. Vac. Sci. Technol. A,
9(4):2431 – 2436, 1991.
[165] Henry Windischmann. Intrinsic stress in sputter-deposited thin ﬁlms. Crit. Rev. Solid
State Mater. Sci., 17(6):547 – 596, 1992.
[166] B. Window. Removing the energetic neutral ploblem in sputtering. J. Vac. Sci. Technol.
A, 11(4):1522 – 1527, 1993.
[167] B. Window and N. Savvides. Unbalanced dc magnetrons as sources of high ion ﬂuxes. J.
Vac. Sci. Technol. A, 4(3):453 – 456, 1986.
[168] Harold F. Winters, H. J. Coufal, and W. Eckstein. Inﬂuence of energy reﬂected from thetarget on thin ﬁlm characteristics. J. Vac. Sci. Technol. A, 11(3):657 – 663, 1993.
[169] W.L.Patterson and G.A.Shirn. The sputtering of nickel-chromium alloys. J. Vac. Sci.
Technol., 4(6):343 – 346, 1967.
[170] L. C. Wrobel. Potential and viscous flow problems using the boundary element method.
PhD thesis, The University of Southampton, 1981.
[171] C. T. Wu, R. T. Kampwirth, and J. W. Hafstrom. High-rate magnetron sputtering of hightc Nb3 Sn ﬁlms. J. Vac. Sci. Technol., 14(1):134 – 137, 1977.
[172] S. C. Wu, H. T. Hsu, F. H. Chen, W. R. Chang, and T. Y. Tseng. Growth parameter eﬀect in superconducting YBa2 Cu4 O8 thin ﬁlms by d.c. magnetron sputtering.
J. Mater. Sci. Lett., 29:5593 – 5598, 1994.
[173] S. C. Wu, W. C. Tsai, C. K. Huang, H. T. Hsu, C. J. Huang, and T. Y. Tseng. Pressure
eﬀect on YBa2 Cu3 O7 thin ﬁlms growth in oﬀ-axis radio frequency magnetron sputtering.
J. Vac. Sci. Technol., A13:2412 – 2419, 1995.
[174] J. Xu, B. M. Moon, and K. V. Rao. Study of the resputtering eﬀect during rf-sputter
deposition of YBCO ﬁlms. J. Mater. Res., 10:798 – 802, 1995.
[175] Y. Yamamura and M. Ishida. Monte Carlo simulation of the thermalization of sputteredatoms and reﬂected atoms in the magnetron sputtering discharge. J. Vac. Sci. Technol. A,
13:101 – 112, 1995.
[176] Yasunori Yamamura. Computer studies of reemission and depth proﬁles for helium on
molybdenum. Nucl. Instr. and Meth., B28:17 – 26, 1987.

BIBLIOGRAPHY

128[177] Yasunori Yamamura and Hiro Tawara. Energy dependence of ion-induced sputtering
yields from monatomic solids at normal incidence. Atomic Data and Nuclear Data Tables,
62(2):149 – 253, 1996.
[178] YiFeng Yang, Ji Ung Lee, and James E. Nordman. Eﬀects of target presputtering on
stoichiometry of sputtered Bi-Sr-Ca-Cu-O thin ﬁlms. J. Vac. Sci. Technol. A, 10(5):3288– 3291, 1992.
[179] P. Ziemann and E. Kay. Correlation between the ion bombardment during ﬁlm growth of
Pd ﬁlms and their structural and electrical properties. J. Vac. Sci. Technol. A, 1(2):512– 516, 1983.