Title Page
 Table of Contents
 Modeling of low dose silicon implant...
 Modeling of the evolution of dislocation...
 Modeling of the stress effects...
 Conclusions and recommendation...
 Biographical sketch

Title: Point-defect-based two-dimensional modeling of dislocation loops and stress effects on dopant diffusion in silicon
Full Citation
Permanent Link: http://ufdc.ufl.edu/UF00082393/00001
 Material Information
Title: Point-defect-based two-dimensional modeling of dislocation loops and stress effects on dopant diffusion in silicon
Physical Description: viii, 145 leaves : ill. ; 29 cm.
Language: English
Creator: Park, Heemyong, 1965-
Publication Date: 1993
Subject: Electrical Engineering thesis Ph. D
Dissertations, Academic -- Electrical Engineering -- UF
Genre: bibliography   ( marcgt )
non-fiction   ( marcgt )
Thesis: Thesis (Ph. D.)--University of Florida, 1993.
Bibliography: Includes bibliographical references (leaves 137-144).
Statement of Responsibility: by Heemyong Park.
General Note: Typescript.
General Note: Vita.
 Record Information
Bibliographic ID: UF00082393
Volume ID: VID00001
Source Institution: University of Florida
Rights Management: All rights reserved by the source institution and holding location.
Resource Identifier: aleph - 001952299
oclc - 31274871
notis - AKC8864

Table of Contents
    Title Page
        Page i
        Page ii
        Page iii
        Page iv
    Table of Contents
        Page v
        Page vi
        Page vii
        Page viii
        Page 1
        Page 2
        Page 3
        Page 4
        Page 5
        Page 6
        Page 7
        Page 8
        Page 9
        Page 10
        Page 11
        Page 12
        Page 13
        Page 14
        Page 15
        Page 16
        Page 17
        Page 18
    Modeling of low dose silicon implant damage effects and oxidation enhanced diffusion
        Page 19
        Page 20
        Page 21
        Page 22
        Page 23
        Page 24
        Page 25
        Page 26
        Page 27
        Page 28
        Page 29
        Page 30
        Page 31
        Page 32
        Page 33
        Page 34
        Page 35
        Page 36
        Page 37
        Page 38
        Page 39
        Page 40
        Page 41
        Page 42
        Page 43
        Page 44
        Page 45
        Page 46
        Page 47
        Page 48
        Page 49
        Page 50
        Page 51
        Page 52
        Page 53
        Page 54
        Page 55
        Page 56
    Modeling of the evolution of dislocation loops and their effects on oxidation enhanced diffusion of boron in silicon
        Page 57
        Page 58
        Page 59
        Page 60
        Page 61
        Page 62
        Page 63
        Page 64
        Page 65
        Page 66
        Page 67
        Page 68
        Page 69
        Page 70
        Page 71
        Page 72
        Page 73
        Page 74
        Page 75
        Page 76
        Page 77
        Page 78
        Page 79
        Page 80
        Page 81
        Page 82
        Page 83
        Page 84
        Page 85
        Page 86
        Page 87
        Page 88
        Page 89
        Page 90
        Page 91
        Page 92
        Page 93
        Page 94
        Page 95
        Page 96
    Modeling of the stress effects on dopant diffusion in silicon
        Page 97
        Page 98
        Page 99
        Page 100
        Page 101
        Page 102
        Page 103
        Page 104
        Page 105
        Page 106
        Page 107
        Page 108
        Page 109
        Page 110
        Page 111
        Page 112
        Page 113
        Page 114
        Page 115
        Page 116
        Page 117
        Page 118
        Page 119
        Page 120
        Page 121
        Page 122
        Page 123
    Conclusions and recommendations
        Page 124
        Page 125
        Page 126
        Page 127
        Page 128
        Page 129
        Page 130
        Page 131
        Page 132
        Page 133
        Page 134
        Page 135
        Page 136
        Page 137
        Page 138
        Page 139
        Page 140
        Page 141
        Page 142
        Page 143
        Page 144
    Biographical sketch
        Page 145
        Page 146
        Page 147
Full Text







Copyright 1993

Heemyong Park


I would like to thank my advisor, Professor Mark E. Law, for many years
of guidance, encouragement, and trust. Not only the invaluable knowledge
and wisdom he shared, but also his way of looking at the world with the
sonorous laugh will never be forgotten and will lead me throughout my life.
I would also like to thank Professor Kevin S. Jones for providing his expertise
on dislocation loop characterization, his systematic experiments, and the
encouragement, without which this work would not have been possible.
I am grateful to Dr. Jim Slinkman in IBM for his substantial support for
the area of dislocation loops and for providing the experimental data on
boron segregation. I also wish to thank Michael Meng for performing the
experiments on loop evolution and sharing the data. Special thanks go to Jim
Listebarger, who helped me with TEM at the later stage of this work.
My life in Florida would have been quite different without the many

people in VLSI TCAD Group who have trusted me enough to listen. More
than anyone else, Minchang Liang deserves my heartfelt thanks for his help
in initiating me and keeping me up on the computers and for the timely dis-
cussions on FLOODS. Thanks also go to Stephen Cea for the help with stress
not only in substrate but also in my mind. I am also grateful to Serdar Kirli,
Xu Huang, and Omer Dokumaci, who came to me as my accidental neighbors
in the office and ended up three of my favorite people on the planet.
I thank my parents for the patience and for implanting in me their
values. Without them, the author of this dissertation would not even exist.
Finally, my most cordial accolade goes to my wife, Jung-Mee, for her love,
support, and everlasting inspiration.

I acknowledge the financial support from SEMATECH.


ACKNOWLEDGMENTS..................................................................................... iii

A BSTR A C T ........................................................................................................... vii


I INTRODUCTION............................................................................................. 1

1.1 The Dislocation Loop in Silicon............................ .............................. 3
1.1.1 Electrical Properties of Dislocations............ ......................... 6
1.1.2 Characterization of Dislocation Loops Created by Ion
Im plantation .......................................................... ........................ 7
1.2 The Effects of Dislocation Loops and Stress in Device Structures....... 10
1.3 Transient Dopant Diffusion Mediated by Point Defects....................... 14
1.4 O organization .............................................................................................. 16


2.1 Damage-Enhanced and Oxidation-Enhanced Diffusion................... 19
2.2 A Model for Dopant Diffusion Based on Pairing Theory...................... 20
2.2.1 Diffusion Equations Accounting for the Dopant-Defect Pairs..... 21
2.2.2 Significance of Binding Energies in Damage-Enhanced
D iffusion............................................................................................ 24
2.3 Defect Equations Incorporating Dual Reaction with Traps.............. 28
2.4 Simulation of Damage-Enhanced Diffusion and OED......................... 34
2.4.1 Initial Distribution of Point Defects................................ ........... 36
2.4.2 Modeling of the Enhanced Diffusion of Boron......................... 38
2.4.3 Modeling of the Enhanced Diffusion of Phosphorus............... 45
2.4.4 Modeling of the Oxidation Enhanced Diffusion....................... 51
2.5 Sum m ary.................................................................................................... 54

BORON IN SILICON..................................................................................... 57

3.1 Modeling of the Dislocation Loops as a Group...................................... 58
3.1.1 The Effective Pressure from the Dislocation Loops of Equal
S ize .......................................................................................................... 60

3.1.2 Asymmetric Triangular Density Distribution of Loop Radius.... 68
3.1.3 Interaction of Dislocation Loops and Point Defects.................. 72
3.1.4 Coalescence and Dissolution of the Dislocation Loops............ 75
3.2 Simulations of the Loop Evolution during Oxidation..................... 78
3.3 The Effect of Dislocation Loops on OED of Boron................................. 87
3.4 Sum m ary................................................................................................... 95

IN SILICO N ................................................................................................... 97

4.1 Importance of Stress in Moder Device Fabrication......................... 97
4.2 A Dopant Diffusion Model Including Stress Effects............................. 99
4.3 Boron Segregation around Dislocation Loop Layer.............................. 109
4.4 Phosphorus Diffusion Retarded by Nitride Film Stress...................... 114
4.5 Two-Dimensional Extension and the Effects on Threshold Voltage
of Short-Channel MOSFETs...................................................................... 117
4.6 Sum m ary................................................................................................... 123

V CONCLUSIONS AND RECOMMENDATIONS........................................... 124

5.1 Summary and Conclusions......................................................................... 125
5.2 Recommendations for Future Research................................................. 129

A PPEN D IX...................................................... .................................................... 132

REFEREN CES ........................................................................................................ 137

BIOGRAPHICAL SKETCH................................................. .................................. 145

Abstract of Dissertation Presented to the Graduate School
of the University of Florida in Partial Fulfillment of the
Requirements for the Degree of Doctor of Philosophy


December 1993

Chairman: Prof. Mark E. Law
Major Department: Electrical Engineering
Dopant diffusion in silicon is studied and modeled on the basis of point
defect kinetics associated with ion implantation damage. Point defect
parameters are extracted from the modeling of transient enhanced dopant

diffusion due to oxidation and low dose implant damage without extended

defects. The theory of dopant-defect pairing is found to be crucial in modeling
the implantation damage effects, and the effective binding energies for boron-
defect and phosphorus-defect pairs are experimentally determined. The
extracted parameters provide an important reference for further modeling of
diffusion under high dose implantation conditions involving extended
Evolution of dislocation loops through their interaction with point
defects is modeled in two dimensions by accounting for the pressure around

the ensemble of loops as well as loop coalescence and dissolution as observed
in transmission electron microscopy (TEM) measurements. Assuming an

asymmetric triangular density distribution of periodically oriented circular
dislocation loops leads to estimation of the effective pressure and an efficient

model for the statistical loop-to-loop interaction. Simulation with the model
correctly predicts variation of the number of captured silicon atoms and the
radii and densities of the dislocation loops during oxidation in agreement
with the TEM data. It also shows significant reduction in oxidation enhanced
diffusion of boron in a buried layer in agreement with measured profiles,
confirming the role of dislocation loops as an efficient sink for interstitials.
A point-defect-based atomistic model for the stress effects on dopant
diffusion is developed by accounting for variation in formation enthalpy of
dopant-defect pairs due to the hydrostatic pressure. Binding energies and
diffusivities of dopant-defect pairs under the pressure are modeled and
incorporated into diffusion equations. Boron segregation around dislocation
loops in silicon is explained by the pressure effects, and the simulation agrees
with the measured profiles. The model also shows that retarded diffusion of
phosphorus under oxide-padded nitride film of various widths is caused by
the stress at the film edge. Two-dimensional simulation of diffusion in the
pressure field leads to better prediction of threshold voltage shift in short
channel MOS transistors.



In the past several years, process simulation has been recognized as an
efficient way to reduce cost and time in development of semiconductor
device manufacturing technologies. The pragmatic goal of simulation is
achieved not only when the process models in the simulators correctly
represent the phenomena observed at each processing step, but also when the
models are capable of predicting unseen effects. Development of the process
models is an ongoing process, more or less following up today's rapid
evolution of technologies. When process simulators go beyond the stage of
mere emulation of reality, they become a practically essential guide to new
technologies. To fulfill this extensive purpose, it is crucial to use accurate and
predictive models founded on what is actually going on inside the material of

The physical models provide process simulators with a potential to
foresee an unprecedented phenomenon in a wide range. If the physics of a
certain processing step is not understood well, phenomenological and
empirical modeling is unavoidable and often produces satisfactory results in a
limited range. When we understand the physical mechanism under good
assumptions, it is the correct choice to develop a physics-based model as far as

computation time due to model complexity is not a concern. Performance of
physical models heavily depends on the parameters that are meaningful in
physical theory. Numerical solution for quantities such as dopant concentra-
tion is usually required with correct values of the parameters. Some
parameters have established values, others can have a certain valid range

expected from theory. In determining the complexity of models to develop,
we need to consider the number of unknown parameters, so that the model
can be parametrized reasonably on the basis of experimental data. Too
complicated a model leads to impractical, arbitrary simulations with a lot of
random parameters, even if it may be based on physical theory in its form.
Parametrization of a physical model is done on the basis of accumulated
data from experiments in a wide range of conditions. Systematic experiments
with advanced measurement reveal the physics of phenomena in the device
materials. Synergetic efforts of material scientists, theorists, and modelers are
more and more required to establish the physical parameters. Validity of
physical models is attested at the same time as the properties of materials
under process conditions are explored. Extensive usage of the physical
models helps to understand the experimental observation and to design
further experiments. Also, limitations of measurement such as two-
dimensional doping profiles can be overcome by the extended use of models.
The models with consistent parameters extracted from a large body of
experimental data eventually lead to accurate and predictive simulations
necessitated for the future technology development.
Dopant diffusion in silicon is one of the fields in process simulation
where judicious, physics-based modeling is strongly demanded. Doping
profile and junction depths in multi-dimensions are the crucial process
design factors that should be modeled most accurately. Diffusion of impurity
atoms, which occurs in the crystal at high temperatures, results in various
doping profiles depending upon the crystallographic, mechanical, thermo-
dynamic, and electrical properties of substrate and interface associated
chemically with the impurities. The diffusion process during the fabrication
of integrated circuits becomes more complicated due to the non-equilibrium

environments initiated by the other processing steps, typically by ion
Modern device fabrication processing includes ion implantation as an
essential technique for doping the semiconductor. The high energy bombard-
ment of incident ions inevitably creates damage in the crystalline structure of

substrate. The ion implantation damage governs the subsequent dopant
diffusion process during thermal annealing cycle, which is necessary for
substrate recrystallization and dopant activation. Transient diffusion of
dopant associated with the ion implantation damage produces relatively
more pronounced impacts in device structures of smaller dimensions. As
devices are scaled down, more accurate prediction of dopant redistribution
due to the transient diffusion is required for designing shallow junctions.

1.1 The Dislocation Loop in Silicon

High-dose ion implantation is crucial in obtaining highly-doped regions
in silicon such as the source and drain of MOSFETs and DRAM cells. The

high-dose implants of common dopants such as arsenic are known to

amorphize the surface region in silicon substrate, simultaneously producing a
large amount of point defects, i.e., interstitials and vacancies. The subsequent
annealing leads to solid phase epitaxial regrowth of the amorphous region,
and the extended defects such as dislocation loops are formed below the
amorphous/crystalline interface. The dislocation loops are inherently
accompanied by a stress field in the crystalline silicon, interacting with the

point defects and the dopant atoms. A simplified picture of a cross-section of
an extrinsic dislocation loop is shown in Figure 1.1. An extra layer of silicon
atoms form an approximately circular shape to attain the lowest self-energy in

the crystalline structure. By intuition, we can visualize the stress distribution


I- Loop Radius R


3-D view

Figure 1.1. A cross-section of an extrinsic dislocation loop composed of an
extra layer of silicon atoms with dangling bonds at the core boundary.

around the loop in a qualitative way. Inside the dislocation loop region, the
atoms tend to be more compactly spaced, thereby causing a compressive
pressure. On the other hand, the region just outside the dislocation loop
boundary is generally under tensile stress, since the atoms are more sparse
due to the misfit. The loop grows by absorption of an interstitial or by
emission of a vacancy, and it shrinks in the reverse way, changing the point
defect distribution around its boundary. Conversely, the non-equilibrium
distribution of point defects around the loop primarily determines the growth
and shrinkage, i.e., the evolution of the dislocation loop.
As the device size is more scaled down, the performance of the shallow
junction, short channel devices becomes more contingent on the process-
induced defects and the related material properties such as stress. The defects
and the stress that arise during the integrated manufacturing processes affect
the characteristics of the scaled-down devices, often in an anomalous way

unobserved or insignificant in the large scale devices of relatively outdated
technology. There are three main aspects in the effects of dislocation loops on
the device characteristics. First, significant leakage currents can be caused by
the dislocations when they are decorated with metallic impurities and located
across the junctions. This is definitely an adverse effect on devices. It can be
suppressed either by avoiding ion implantation for the junction formation or
by removing or relocating the formed dislocations deliberately through a

subsequent process. As far as the dislocations are confined to areas outside
the space-charge regions, they do not usually have detrimental effects on
device operation.
Second, the dislocation loops indirectly affect dopant redistribution by
capturing and emitting point defects during subsequent anneals. Since
dopant atoms diffuse by pairing with point defects, the dopant redistribution
and the resultant doping profiles and final junction depths are largely affected
by the dislocation loop evolution in association with the point defects. The
role of dislocation loops as capturing source for interstitials can be utilized to

reduce the junction depth during oxidation and annealing.

Third, the stress field from the locally disturbed crystalline structure
directly influences diffusion of dopant atoms paired with the point defects.
The stress arises not only from the extended defects, but also from thin films
and device isolation structures. Especially the dislocation loops can getter
dopant atoms as well as silicon interstitials and metallic impurities in the
pressure field around their boundary, causing dopant segregation. The
superimposed interactions of dopants, point defects, and extended defects
make the diffusion process more complicated. In this work, the second and
the third aspects of the effects of dislocation loops are investigated by building

up physics-based atomistic models of the dislocation loop evolution and the
stress effects on dopant diffusion in general.

1.1.1 Electrical Properties of Dislocations

The electrical properties of dislocations in semiconductors have been
studied mainly on the theoretical basis. Shockley [1] remarked that the
dangling bonds in the core of an edge dislocation may form a one-
dimensional partially-filled band of edge-states analogous to two-dimensional
surface state bands. The broken bonds in the extra half-plane of atoms are

considered to be occupied by only one electron in a neutral crystal. Read [2, 3]
developed a model postulating that the dislocations act as acceptor centers,
which explains the observed reduction of the carrier density in n-type Ge with
plastic deformation [4, 5]. He assumed that all conduction electrons are
expelled from the vicinity of a negatively charged dislocation in n-type
material and a cylinder of positively ionized donor atoms is created around
the dislocation core. Labusch and Schr6ter [6, 7, 8] presented a different model
based on the assumption of half-filled dislocation band in the neutral state,
which accounts for the donor-like behavior of dislocations in p-type Ge and Si
at very low temperatures observed through Hall measurements of holes.
However, all these models were only about the charge states of dislocations at
temperatures much lower than common annealing temperatures used in
Recently, Ross et al. [9] measured the reverse leakage current of GeSi/Si
p-n diode during the in situ formation of misfit dislocations at about 7000C,
simultaneously characterizing the dislocation density and length by
transmission electron microscopy (TEM). They observed the proportionality

between the leakage current and the dislocation density and length, and

found that a generation-recombination process at the dislocation cores is not
sufficient to explain the substantially large generation current. Consequently
it was suggested that device degradation due to the dislocations should be
related to the point defects or the metals diffusing rapidly along the
dislocations. The same argument may also be applied to the implantation-
induced dislocation loops, for they are essentially a sort of edge dislocations
with a morphological change. Cerva and Bergholz [10] showed that the half-
loop dislocations at the oxide mask edges act as effective nuclei for the
formation of Ni precipitates, which are originated from the reactive ion
etching process step. In that way, the dislocation loops can cause fairly high
leakage current, unless their formation is suppressed by such means as
implanting carbon [11] which has been conjectured to be gettering centers for
Si self-interstitials [12].

1.1.2 Characterization of Dislocation Loops Created by Ion Implantation

Unlike the point defects, the extended defects such as dislocation loops
are usually large enough to be observed as they are formed inside the silicon
substrate. The transmission electron microscopy (TEM) produces a manifest
view of the extended defects with principal sizes ranging around several
hundred angstroms on the average. From the modeler's point of view, this is
a critically aiding factor in building up a predictive model based on
experimental observation. Jones et al. [13] systematically analyzed the
morphological characteristics of the extended defects due to ion implantation
in silicon by using the TEM. Over 2500 different plan-view and cross-
sectional TEM specimens have been examined to investigate the effects of
implant species, dose, energy, annealing time and temperature, wafer
temperature and orientation, and pre- and post-amorphization on the defect

formation. The classification scheme groups all secondary defects into five
categories based on the origin of the damage. Of those five categories, the

categories I and II refer to the dislocation loops induced by typical ion
implantation conditions used in silicon device fabrication. The major criteria
that distinguish the formation of the two categories of extended defects were
found to be implant dose and implanted ion species, i.e., ion mass [13].
Figure 1.2 shows the criteria of dislocation loop formation due to implants of
common dopants in silicon, extracted by Jones [13]. Category I damage is
subthresholdd" damage that results when the implant dose exceeds the
critical dose (~ 2x1014 cm-2) and simultaneously no amorphous layer is
formed. These defects are typically extrinsic dislocation loops that are
observed at a depth corresponding approximately to the peak of the
implanted dopant distribution. For room temperature implantation, this
damage is usually induced by implants of light ions such as boron, and its
density is found to be proportional to the dose.
If the dose is increased sufficiently to cause amorphization of the surface
region of silicon substrate, then defects classified as category II ("end of range")
damage are formed beneath the amorphous-crystalline interface in the
heavily damaged but still crystalline material. The category II defects evolve
into extrinsic dislocation loops after the regrowth of the amorphous layer
with the subsequent annealing. This damage arises whenever an amorphous
layer is formed during implantation. High dose implants of heavy ions such
as arsenic are bound to induce the category II dislocation loops in silicon
devices. In addition, half loop dislocations are observed to form upon
annealing after arsenic implants with very high dose larger than 5x1015 cm-2,

evolving from damage classified as category V. The category V defects,
usually precipitates or dislocation loops, occur when the solid solubility of the

Criteria of Extended Defect Generation

U Critical dose for Category I defects
"i amorphization
E 1o16
0o .. Category H defects

I Category I threshold dose

S- Category H threshold dose
No extended defects

0 20 40 60 80 100 120 140
Ion mass 4amu)

50 ~190 keV room temperature implants
~5jA/cm2, (100) wafers [Data by Kevin S. Jones]

Figure 1.2. Ion mass-dependent critical implant dose of common dopants in
silicon for generation of category I and II extended defects. Data by K. S. Jones.

implanted species at the annealing temperature is exceeded. In this thesis, the

modeling is focused on the category II dislocation loops, based on their

observation through TEM.

Figure 1.2 shows that the dislocation loops are not formed as far as the

implant dose is below the critical doses (~ 2x1014 cm-2 for B, P, and Si; 5x1013

cm-2 for As and Sb). In this regime of implant conditions, an excessive

amount of point defects are created after the implantation, and they cause

transient diffusion of dopants without nucleating extended defects during the

subsequent thermal processes. It is this regime that should be used to

investigate the non-equilibrium diffusion mechanism of dopant atoms and

point defects in the absence of their complicated interaction with extended


1.2 The Effects of Dislocation Loops and Stress in Device Structures

The effects of the dislocation loops on device performance were
previously observed in shallow junction bipolar technology [14, 15]. The
dislocation loops emanate from the emitter-base junctions due to the emitter-

formation arsenic implants and also from the heavily implanted buried layer.
The dislocations decorated with metallic impurities can be conductive
enough to cause emitter-collector leakage. When they extend through the
narrow base, they become emitter-collector pipes which allow significant
emitter-to-collector leakage current even when the base terminal is open.
This phenomenon ultimately gives rise to severe yield problems [15]. Since
the growth of dislocation loops is critically determined by the point defect
distribution around them, a point defect based model will provide the most
fundamental estimates on the tolerable size and density of the dislocations
under a certain processing condition for a given device structure.
It is well known that the stresses from traditional LOCOS (local
oxidation of silicon) isolation structures are apt to induce dislocations in the
silicon substrate. In scaling down the device size, the schemes for reducing
the extent of laterally merging field oxide, so called "bird's beak," tend to
make the LOCOS isolation more susceptible to the stress-induced dislocation
generation. Fahey et al. [16] investigated the effects of dislocations and stress
observed in the recent process technologies including the trench isolation and

trench capacitor structures suitable for a bipolar IC and 4-Mbit and 16-Mbit
DRAM (dynamic random access memory) fabrication. The trench isolation
method, composed of reactive ion etching (RIE) of the substrate and chemical
vapor deposition of SiO2 for filling the trench, was developed to avoid the
scaling limitation of the LOCOS-based isolation techniques. Although it

achieves good planarity and scalability, substantial stresses can be generated
during thermal oxidation steps usually required following the planarization.
The thermal oxidation of the trench sidewalls, accompanied by a substantial
volume expansion, is analogous to driving a wedge into the trench,
overfilling it, creating stress and defects such as dislocations in the adjacent
silicon [16]. The shape of the trench and the temperature of oxidation are
critical factors in the reduction of stress-induced defects. At the corner
locations of the trenches, larger stresses can be generated. Stress build-up
becomes worse as the oxidation temperature is decreased.
Figure 1.3 shows a simplified cross-section of a trench capacitor DRAM
cell with the substrate-plate-trench structure [17]. Without proper control
over the processes, the leakage current may be a very significant problem in

Bit line

/ Word line

Dislocation loop




po+ poly

in oxide





Figure 1.3. A simplified cross-section of a trench capacitor DRAM cell; ref.


this structure due to superimposed effects of the implantation-induced
dislocation loops near the p+n junction and the stress at the corer of the

trench capacitor. The stresses originate from the oxidation process, from
thermal expansion mismatch of the silicon and the oxide, and from the
intrinsic compressive stress of the polysilicon that fills the trench [16]. Under

the stress from the trench, the dislocations can move great distances outside
the original implanted region by the gliding process, causing a detrimental
effect on device performance. The excessive leakage currents increase power
dissipation of the circuits and also degrade refresh times in the DRAMs.
Detailed stress analysis combined with TEM observations is required to
examine the defect generation and propagation mechanism.
The MOSFET device characteristics are also influenced by the
dislocation loops due to high-dose ion implantation as well as the stress-
induced edge dislocations during the oxidation processes. Figure 1.4 shows a
simplified cross-section of a typical n-channel lightly-doped drain (LDD) MOS

structure, which features the self-aligned phosphorus doped n- regions
between the channel and the n+ source and drain regions doped by high-dose
arsenic implantation. During the annealing cycle for activating the
source/drain region, a layer of dislocation loops is produced near the
source/drain-to-substrate junctions. The dislocation loops, which can cross
the junctions, may trap the fast-diffusing impurities and provide generation-
recombination centers at their sites. This can lead to large leakage currents
through the reverse-biased source/drain-body junctions, as in the case of
misfit dislocations observed by Ross et al. [9] The excessive leakage currents
increase power dissipation of the integrated circuits as a whole. Furthermore,
reduction in the junction breakdown voltage is expected. Minority carrier

Polysilicon gate Si02

Al ........H ..1 ... Al

.n -; n\
nn n+

Dislocation loop layer I I Antipunch implant
Sp-type substrate
LDD region

Figure 1.4. A simplified, unsealed cross-section view of an n-channel LDD
MOSFET structure.

lifetime is also decreased by the extended defects through the introduction of

localized energy levels within the silicon bandgap.

As the channel length of the MOS transistor scales down below 2 Pm, a

series of effects arise that can not be predicted by the conventional long

channel device models. One of the short channel effects is the decrease of

threshold voltage VT due to the contribution of the depletion regions of the

source and drain junctions on the channel depletion region charge. The

primary factor that determines the threshold voltage is the background

doping concentration of the substrate or the wells of CMOS. Recently, Sadana

et al. [18] found that a pronounced segregation of boron from the bulk silicon

into the arsenic-implanted region occurs during the subsequent anneal. The

local redistribution of dopant in the substrate can make a significant

difference in the VT roll-off in the short channel MOSFETs [18]. The observed

boron pile-up was exactly centered around the arsenic implant-induced

dislocation loop layer, which strongly suggests the direct interaction of dopant
atoms with the dislocation loops. A similar local redistribution was also
reported of the implanted boron profiles in the post-amorphized silicon
substrate [19]. The stress field around the dislocation loops is supposed to be
the major source of change in the dopant diffusion mechanisms. It is also
possible that the capture of dopant impurities takes place in a way similar to
the gettering of fast-diffusing impurities through the extended defects. In this
thesis, it is attempted to model this phenomenon on the basis of the stress
field and an atomistic dopant diffusion theory.

1.3 Transient Dopant Diffusion Mediated by Point Defects

According to the generally accepted theory [20], dopant diffusion in
silicon occurs through pairing with point defects, both vacancies and
interstitials. The point defects and dopant diffusion in the absence of the
extended defects should be clarified first in order to develop an extensive
model for the dislocation loop effects. It is difficult to extract point defect
parameter values from experiment, since they cannot be measured directly.
Only by examining the increase or decrease in dopant diffusion can point
defects be investigated. Oxidation-enhanced diffusion (OED) is the most
widely measured and understood non-equilibrium diffusion phenomenon,
which makes it a suitable starting point for understanding point defect
transport and recombination. It is believed in general that OED is caused by
injection of interstitials from the growing oxide.
One of the OED effects on the short channel device characteristics can be
seen in the reverse short-channel effect [21, 22]. It was observed that the
threshold voltage initially increases as the channel length decreases until the
final VT fall-off begins. The anomalous VT behavior was attributed to the

OED effects on redistribution of the channel dopant (boron in n-channel
devices) during polysilicon gate sidewall reoxidation. A typical channel
doping profile has a positive concentration gradient towards the SiO2/Si
interface, since an antipunch implantation is usually used to suppress short
channel effects (Figure 1.4). The interstitials injected laterally during the
oxidation causes enhancement of the channel dopant diffusion, which results
in a local increase of the surface concentration near the source and drain
regions. Additionally, it is suggested that the channel dopant redistribution
around the dislocation loops near the source/drain junctions can also
complicate the situation regarding VT variation in the short-channel
MOSFETs, as mentioned in the previous section.
It has been reported [23, 24, 25] that low dose silicon implantation also
provides non-equilibrium diffusion conditions by creating point defects or
defect busters. If the silicon implant dose is below 2x1014 cm-2, no extended
defects are formed [13]. The effects of point defects on dopant diffusion in
intrinsically doped silicon are crucial especially in the processing of shallow
junctions and short-channel devices. The lightly-doped drain regions shown
in Figure 1.4 are made typically through phosphorus implantation with a
dose range of 1013 cm-2, which produces the junction depth and length
dimensions less than 0.3 gm. The breakdown voltage and the threshold
voltage in the LDD device are strongly dependent on the dopant
redistribution of the lightly doped regions [26]. It has been found that the
formation of arsenic/phosphorus junctions in the LDD is influenced by the
point defects created by the low-dose phosphorus implants [27]. The damage-
enhanced transient diffusion of phosphorus can not be explained by the
traditional diffusion theory of gaussian distribution and Fick's law. To

correctly model the diffusion, it is necessary to assess the pairs of dopant and
point defects according to their electrochemical reaction.

1.4 Organization

The primary goal of this work is to understand the dopant diffusion
phenomena related with ion implantation damage. Based on the knowledge
on the diffusion mechanisms, physics-based models are developed for the
interactions among dopant atoms, point defects, and extended defects.
Material-related effects such as traps and mechanical stress are also dealt with
in this thesis. This work unveils several crucial points in diffusion studies,
and uses them in simulating the experimental observations correctly.

Predictive and quantitative simulations are achieved by implementing the
physics of diffusion.
Figure 1.5 shows a schematic diagram for the interrelation of the three
main elements of diffusion, i.e., dopant atoms, point defects, and extended
defects. In this work, the extended defects refer only to dislocation loops of
category II origin. The interaction between each of the elements, signified as

A, B, and C, corresponds to the topic of each chapter in this thesis. All of the
three elements can be created after ion implantation through the process
denoted as a, P, and y, and they go through nucleation, evolution, or
diffusion processes during the subsequent thermal annealing with different
ambient conditions, e.g., oxidation, nitridation or inert annealing. By
deliberately controlling the conditions of implantation and annealing, we can
selectively create the elements of diffusion to achieve an uncomplicated

experimental environment suitable for probing each interaction between two
of the elements. It is also required to avoid any extrinsic diffusion of dopant
at high concentration above solid solubility limit (8 in Figure 1.5) involving

Figure 1.5. The interrelations of the three major elements of diffusion
conditioned by ion implantation.

precipitation and/or electric field effects, which is not the subject of this
Chapter II deals with transient enhanced diffusion of boron and
phosphorus through pairing with the point defects induced by silicon ion
implantation and oxidation in the absence of extended defects (interaction A).
The point defect based diffusion model in the process simulator SUPREM-IV
is re-examined by implementing the pairing of dopant and point defects
under excessive supersaturation of point defects due to ion implantation.
The pairing model clarifies the thermal nature of the diffusion enhancement,
and leads to a physically-meaningful quantity that governs the transient
diffusion. A set of point defect parameters including effective binding

energies of dopant-defect pairs is extracted by simulating the experimental

data, and it is consistently used for modeling diffusion phenomena under
extended conditions in Chapter II and IV.

In Chapter III, the implantation-induced dislocation loops are modeled
based on TEM experiments, and their evolution associated with point defects
during oxidation (interaction B) is discussed and correctly simulated in
FLOOPS (Florida Object-Oriented Process Simulator). Furthermore, the
combination of the dislocation model and the pair diffusion model in
Chapter II provides a consistent simulation of the indirect influence of the
dislocation loops on boron diffusion through the redistribution of point
defects around the layer of loops (interaction A).

Chapter IV describes an atomistic model for stress effects on dopant
diffusion in general. A pressure-dependent dopant diffusion equation is
derived by accounting for the variation in enthalpy, binding energy, and
diffusivity of dopant-defect pairs. The model leads to quantitatively
consistent simulation of the boron segregation around the dislocation loop
layer (interaction C) and the nitride film stress effects on phosphorus
diffusion. The effects of the stress and the boron segregation around the

dislocation loops on device characteristics are also discussed.
Finally, this work is summarized in Chapter V. Conclusions and
recommendations for future extension of this work are also made.


This chapter describes point defect models and compares them with
experimental results for intrinsically doped material. Transient dopant
diffusion due to low dose silicon implant damage can be modeled with the
same parameters as oxidation enhanced diffusion, and therefore provides an
additional technique to probe point defect behavior. Parameters are extracted
consistently for both experimental conditions and fit to Arrhenius
relationships. The theory of dopant-defect pairing is found to be crucial in
modeling the implantation damage effects, and the effective binding energies
for boron-defect and phosphorus-defect pairs are estimated from the

2.1 Damage-Enhanced and Oxidation-Enhanced Diffusion

The low dose implant damage enhanced diffusion and the OED should
be directly comparable. There is generally a large difference in defect amount
and distribution between the two cases. The concentration of point defects
just after implantation usually reaches several orders of magnitude larger
than the equilibrium concentrations. Interstitials and vacancies created

locally at the surface region diminish rapidly through recombination during
the initial annealing period. It is this short time period when the transient

diffusion of dopant atoms occurs, and normal diffusion is resumed after this
transient period. On the contrary, oxidation is believed to induce only
interstitials from the Si/SiO2 interface. The interstitials are injected during

dry oxidation to a much smaller amount than the implantation-induced
defects. However, the injection takes place continuously as long as the
oxidation goes on, and the supersaturation of interstitials persists and usually
extends to the deep bulk silicon region through uninterrupted interstitial
diffusion. As a result, the OED proceeds almost constantly during oxidation.

The modeling of oxidation effects on intrinsically doped layers has been
previously described by Law [28]. It was possible to fit the OED experimental
results with two different sets of parameters obtained by using two different
measurements of the interstitial diffusivity under gettering conditions [29,
30]. Park and Law [31] developed a point defect based model for the enhanced
diffusion of boron and phosphorus due to the low dose silicon implantation
damage. This chapter describes the model and simulation results in ref. [31]
extensively. A new set of model parameters was extracted to consistently fit
both the experiments of implant damage and of OED. The effect of both of
these conditions creates non-equilibrium diffusion, and a good test of a
dopant diffusion model is whether it can be used for both conditions.

2.2 A Model for Dopant Diffusion Based on Pairing Theory

This modeling work is based on the assumption that the point defect-
dopant pairs are at local equilibrium, as has been used previously [20, 32 ~ 39].
Other researchers [40, 41, 42] have developed models that do not depend on

this assumption. However, these models involve much more computation

time since they have to solve more sets of partial differential equations. They
also have to estimate both the pairing and depairing reactions rates. Given

the present state of knowledge on the pairing mechanism of charged defects
and dopant atoms in non-equilibrium, it seems premature to determine the

forward and reverse rates of pairing for different charge states through
modeling. For example, it can be very arbitrary to extract the values of
binding energies of dopant-defect pairs with different charge states, and the
model can be prone to many assumptions. It is therefore useful to see if the
local equilibrium assumption allows modeling of implant damage effects.
While this work was in progress, other researchers [43, 44] were
developing a diffusion model under another different assumption. In
addition to the dual mechanism of diffusion with interstitialcy and vacancy,
they assumed other possibilities of dopant-defect reaction such as Frank-
Turnbull mechanism recombinationn of vacancy with dopant-interstitial pair)
and dissociation mechanism recombinationn of interstitial with dopant-
vacancy pair). Those mechanisms have been used in explaining diffusion of
fast-diffusing metallic impurities. By applying the two additional
mechanisms to diffusion of common dopants in silicon, they modeled the
diffusion of arsenic and phosphorus under nitridation condition [43]. In the
case of implant damage enhanced diffusion with the excessively large
amount of point defects, however, it is not certain whether Frank-Turnbull
mechanism is applicable to the common dopants. Estimation of binding
energies is required to model the transient diffusion of dopants in silicon
under excessive supersaturation of point defects as created by ion
implantation. In this chapter, both the damage enhanced diffusion and the
oxidation enhanced diffusion of boron and phosphorus are modeled by the

dual mechanism of diffusion with pairing theory.

2.2.1 Diffusion Equations Accounting for the Dopant-Defect Pairs

When it is assumed that local equilibrium is attained between the
dopant, the defects, and the dopant-defect pairs at each point in the diffusing

dopant profile, the following equations can be applied to modeling of
damage-induced transient diffusion. The dopant equations proposed by
Mathiot and Pfister [33] are modified to limit the maximum number of
dopant-defect pairs as described in Law et al. [45]:

-JAX dAcKAxGx-Ccxc))Cl C Vlog CA+ C 2-1
v c' / x \e Cx n

where JAX is the diffusive flux of donor dopant A paired with defect X which
is I (interstitial) or V (vacancy), Gx' is the charge state term (n/ni)1-c, where c
denotes charge state of the pair (0 for neutral pairs, -1 for negative pairs, and
+1 for positive pairs), Cxc- is the concentration of interstitials in a charge state
one electron more negative than the pair state c under inert intrinsic
condition, CA is the concentration of unpaired dopant atoms in the
substitutional sites, and dA is the diffusivity of AXc pair. The consideration
of charge states leads to the expressions for concentration of defect X under
non-inert condition in terms of C~x- as follows [46]:

Cx-1 = Gxc x- 2-2

where Cx = GxC'xc, and Cx = C Cxe 2-3, 2-4

The above equations were derived for a donor dopant of singly positive
state. Similar equations are applied to an acceptor dopant case when charge
states of defect and dopant are considered. When we assume that electronic
processes are faster than any other and that the dopant atoms and defects are
both dilute compared to the silicon lattice sites, the number of dopant-defect
pairs can be expressed by a simple mass action relationship with a pairing
coefficient KAXC:

CAX = KAX CX-1 CA 2-5

The relationship of the pair concentration to the total concentration of

dopant CA can be expressed as:

CA= CA + Y KAIcCIC CA + 1 KAV CV- CX = CA + CAI + CAV 2-6

where the summation is over all pair charge states and CA is used for the

unpaired substitutional donor atom concentration.

According to statistical thermodynamics based on the dilute

concentration approximation, the pairing coefficient KAXc in Eq. 2-5 is
expressed as [20, 32, 47]:
0 EbAx(\
AXc bAXc
KAXC = Cs exp WkTJ 2-7

where Cs is the concentration of lattice sites, EbAXc is the binding energy of the
AXc pair, and OAXc is the coordination number which represents the number
of equivalent ways of forming the AXC state at a particular site, which is
assumed to be 4 for both dopant-vacancy and dopant-interstitial pairs in

silicon. In fact, Eq. 2-7 is based on diffusion theory via a vacancy mechanism

without considering the charge states of the defects. However, it has been
extended self-consistently to an interstitial(cy) mechanism [20], and it is also

reasonable that the binding energy EbAXc depends on the charge states of the

defects. In intrinsic material, it is impossible to distinguish the charge states,
and an effective pairing coefficient KAX is used:

= KAX CX-1
C ex

The effective pairing coefficient KAX is determined by temperature and the
intrinsic charge state distribution of defects in intrinsic conditions. Therefore,
we can define effective binding energy EbAX of dopant-defect pairs AX in the
same way as Eq. 2-7:

OAAX (EbAX 2-9
KAX = S exp kT 2-9

Thus, the effective binding energy EbAx incorporates the charge distribution of
defect X, which is constant under intrinsic doping conditions. When doping
is low and moderate, the EbAX is a physical parameter that describes the total
average energy required to separate the dopant atom and the defect from a
pair state. The value of EbAx can be estimated by simulating the enhanced
diffusion in the lightly-doped layer, as will be demonstrated in the remaining
part of this chapter. When the doping is extrinsic, the actual binding energies
of the pair AXc may vary significantly depending on the charge state of the

2.2.2 Significance of Binding Energies in Damage-Enhanced Diffusion

Since ion implantation creates such a large number of excess defects, it
becomes critical to account for the dopant-defect pairs. During the initial
short period of transient diffusion, complete pairing of all the dopant atoms

may occur, so that CA = CAI + CAV. The enhancement of dopant diffusivity
under these conditions can be investigated by deriving an expression for

DA/DX in terms of effective pairing coefficients KAI, KAV, and defect
concentrations. With the assumptions that defects are dilute compared to
lattice sites of silicon and that the defect gradient is negligible, Eqs. 2-1, 2-6, and
2-8 can be used to derive an equation for the instantaneous diffusivity
enhancement of a dopant A under intrinsic conditions:

fAI + ( 1 AI )

A I C + +n t IKAIC +KAVC] 2-10
/AI + -fAOc 1 + KAICI + KAV Cv2

where fAI is the fractional interstitial component of diffusivity in
equilibrium. To simplify further, it can be assumed that the dopant diffuses
by interacting only with interstitials, i.e., fAI = 1 and KAV = 0, and Eq. 2-10

DA KAICI +(CI/C 2-11

This formulation qualitatively explains the transient enhanced diffusion of
phosphorus and boron, since they are known to diffuse mainly via an
interstitialcy mechanism. There are two limiting conditions to be considered:

(a) When the supersaturation of interstitials, CI / Cl, is moderate and
below a threshold (KAIC ) 1, i. e., when KAICI < 1, the diffusivity
enhancement DA / DI is proportional to CI / C*. This is the situation during
(b) When CI / Cl exceeds the threshold (KAI C ) 1, i. e., when
KAI CI > 1, Eq. 2-11 is approximately
S= 1+ (KAIC)-1 2-12

Eq. 2-12 characterizes the damage-enhanced diffusion during the initial
transient period. Under intrinsic doping conditions, Eq. 2-12 indicates that
the diffusivity enhancement is no longer dependent on the interstitial
supersaturation, and is determined only by temperature when the
supersaturation is higher than the threshold (KAI C1) 1 that is also

thermally determined. Therefore, this model predicts that the short time
diffusion enhancement is determined solely by the KAI C\ product, and it is
this product that is critical in modeling the ion implantation damage effects.
Figure 2.1 shows the diffusivity enhancement DA / D; as a function of
interstitial supersaturation CI / Ct, calculated from Eq. 2-11 for three different
values of KAI. The instantaneous diffusion enhancement DA / D* is limited
by the value of threshold ( KAI C ) -1 when CI / C\ is excessive, in contrast to
the case of conventional modeling that neglects the amount of dopant-defect
pairs CAI. When the pairing theory is ignored in modeling the transient
diffusion under supersaturation of point defects, the diffusion is erroneously
overestimated in the calculation. The result would be obviously wrong as the
straight line of proportion in Figure 2.1 suggests. Simulations without

1 0 8 I I " 11111 I I 1 1 i , ,I , ,,,,,,I , .,.,I , ,, ...,I ,I .il I , ,,,,,,I

-- KAI=5x10-17 c3
- KAI=5x10-16 cm3
10 -- _- KAI=5x10-15 cm3
-- -with CAI neglected

0< 104

C/ -I* ---

1 I i I II I I II I I llli I t IIII- |l I I I I I I I I 1 m I I "Ii I| I mIlI I
1 102 104 106 108

Figure 2.1. Enhancement of instantaneous diffusivity of dopant dominated
by interstitialcy mechanism (fAI = 1) as a function of supersaturation of
unpaired interstitials, with and without considering the dopant-defect pairs.

considering the pairing model actually showed excessively enhanced
diffusion of phosphorus and boron with implant damage.
We can draw from Eq. 2-12 a very important fact about the equilibrium
concentration of interstitials and the binding energy of dopant-interstitial
pair. Experiments show that the enhancement of phosphorus and boron
diffusion is always larger at lower temperatures [23, 24, 25]. The enhanced
diffusivities measured from the experiments are time-averaged effective
diffusivities, which represents only the net amount of diffusion shown final
profiles [48]. However, it is fairly reasonable that the total diffusion amount is
determined by the enhancement of the instantaneous diffusivity during the
initial time period with high supersaturation of defects, as shown in Eq. 2-12.
From Eq. 2-12, therefore, it is inferred that the value KAI CI should be smaller
at lower temperatures. The KAI is determined by the binding energy EbAI in
the relationship of Eq. 2-9. A theoretical expression of the equilibrium
interstitial concentration C\ can be obtained from a general model of
equilibrium concentrations of point defects in statistical thermodynamics [20]:

C 0 =- I Cs exp exp - 2-13

where Cs is the number of available lattice sites in the crystal, HI is the
formation enthalpy of silicon self-interstitial, Sf is the formation entropy that
is usually attributed to lattice vibrations, and O8 is the number of degrees of
internal freedom of the defect on a lattice site. As mentioned in Fahey et al.
[20], no experiment has definitively measured the equilibrium concentrations
of vacancies or interstitials in silicon, or even the enthalpies of formation.
C*j in Eq. 2-13 does not include the effects of multiple charge states of the
interstitial and the Fermi level shift, which is insignificant under intrinsic

doping conditions. Combining Eq. 2-9 and Eq. 2-13, we obtain the relationship
of the factor KAI C}I with temperature and energies as follows:

KAI CI = OAI 1 exp ) exp (.- .EbAI k 2-14

To be consistent with the observed temperature dependence of enhanced
diffusion, therefore, the binding energy EbAl should be smaller than the
formation enthalpy of interstitial HI, i.e., the activation energy of the
equilibrium concentration of interstitials. This is expectable since the energy
difference H{ EbAI is equal to the formation enthalpy of the dopant-
interstitial pair HAI, which should be positive in order to be physically
meaningful. Thus, the product KAI CL, a critical factor determining diffusion
enhancement, should be considered as the proper measure of the activation
of dopant-interstitial pair. With larger HAI, the KAI Ci is smaller, and the
concentration of the pairs in thermal equilibrium is also smaller. It brings
about more enhanced diffusion in non-equilibrium with high super-
saturation of interstitials. In addition, the large value of HAI results in
stronger dependence of the enhancement on temperature. The same
argument based on the energetic can be applied to the vacancy mechanism,
and the dopant-vacancy binding energy EbAv should be smaller than the
vacancy formation enthalpy Hy. These relationships provide a crucial
restriction on physics-based estimation of the point defect parameters.

2.3 Defect Equations Incorporating Dual Reaction with Traps

It has been observed that diffusion of point defects depends on the
silicon material used in the diffusion studies. Interstitial diffusivities in the
literature vary by several orders of magnitude from experiment to

experiment. It was attributed to the difference in material and silicon crystal
growth method used in each experiment. Impurities such as carbon and
oxygen in silicon are known to interact with point defects and thereby affect

their diffusion. Griffin and Plummer [49] explained the material dependence

of interstitial diffusivity by introducing the concept of bulk traps capturing

interstitials. They considered the reaction of an empty bulk trap T (i.e.,
unpaired oxygen or carbon atom) and an interstitial I:

I + T IT 2-15

The above reaction accounts only for the reverse reaction as the dissociation

of the filled trap IT (i.e., interstitial-trap pair). However, vacancies are also
expected to contribute to the trapping of interstitials. In this study, an

additional interaction of vacancy and filled trap is modeled to explain the

trap-mediated recombination of point defects:

V + IT T 2-16

In Eq. 2-16, the forward reaction represents the recombination of a vacancy

and a trapped interstitial, while the reverse reaction is the trapping of a

silicon atom out of its substitutional site, leaving a vacancy behind.
The empty trap concentration CET is computed by the following

equation that incorporates both interstitial and vacancy trap interaction


S= RTI + RTV 2-17

where RTI and RTV are the rate of empty trap concentration change due to the
interstitial reaction in Eq. 2-15 and the vacancy reaction in Eq. 2-16,

respectively. They can be derived by considering equilibrium states for the
traps and each kind of point defects:

RTI = Ktrap [CET T -C C (CT CET)] 2-18

RTV = KtrapV (CT CET) CV CT ET C* CET] 2-19

where CT is the total trap concentration which depends on the silicon
material, and CET is the empty trap concentration in equilibrium. KtrapI is the
rate of forward reaction between interstitials and empty traps, while KtrapV is
that of vacancies and filled traps. Eq. 2-18 is similar to that given in Law [28],
but the trap equation (Eq. 2-17) now includes the vacancy reaction term.
The trap reaction rates KtrapI and KtrapV are limited by the diffusivity of
each kind of point defects, provided that the traps are immobile and
unlocalized. Moreover, there should be a certain amount of energy barrier in
the trap interaction, similarly to the case of direct recombination of interstitial
and vacancy as described by Fahey et al. [20]:

KR a=va (DV + Dv) exp AEI 2-20
KR Cs kT

where AEIv is the recombination energy barrier of free interstitial and
vacancy, aiv is the capture radius for bulk recombination, 0 is the volume of
the silicon unit cell, CS is the density of lattice sites. DI and DV are the
diffusivities of interstitial and vacancy, respectively. In analogy to Eq. 2-20,
the trap reaction rates can be formulated as follows:

47C aEx-I _AEET- I 2-21
KtrapI = CT1 D exp AEET 2-21

4xc asrv AEFT-V
Ktrap =4 a DV exp A-kT) 2-22
Kap TCs

where AEET-I and AEFr-V are energy barriers of recombination of empty trap-
interstitial and filled trap-vacancy, respectively. The corresponding capture
radii aET-I and aFT-V can be assumed to be equal to aiv for direct I-V
The defect continuity equations used in the model are based on those
given in Law and Pfiester [27] which include terms for bulk recombination
and reaction of empty traps and interstitials. In this work, they are modified
to incorporate the above model for the dual interaction of traps with both
interstitials and vacancies. The equations for interstitial and vacancy are:
a(CI + CAI)
= V(JI + JAI) KR (CI CV C CV) + RTI 2-23

(Cv + CA) = V(Jv + JAV) KR (CI Cv CCv) RTV 2-24

where JI and JAI are the flux of free unpaired interstitials and that of
interstitials paired with dopant A, respectively. RTI is the trap mediated
recombination affecting the interstitial amount as given in Eq. 2-18. In
Eq. 2-24, Jv, JAV, and RTV are similarly defined. Assumption of detailed
balance shown in Eqs. 2-18 and 2-19 is represented in the above defect
continuity equations by separating the trap interaction with interstitial and
vacancy as RTI and RTV.
The boundary conditions for the interstitials are [50]:

DI V CI + KI(CI-C) = g 2-25

where KI is the surface recombination velocity, and g, is the injection flux. A
similar equation holds for vacancies, but gv is zero for oxidizing conditions.
The interstitial injection flux at an oxidizing interface, gl, is assumed to be

proportional to the number of silicon atoms consumed by the oxidizing

interface, and the proportionality constant 0 is the fraction of consumed
silicon lattice atoms that are re-injected into the crystal as interstitials. The

surface recombination velocity, KI, can also depend on the surface growth

rate. The surface recombination is expressed as:

KI = KImax 0ox +1 Kmin

where vOx is the initial oxide growth velocity for a bare silicon wafer, KImax is

the surface recombination velocity maximum at a growing interface, Kimin is
the velocity found at the inert interface, and ai is the decay dependence.
The defect diffusion model shown in this section is tested by simulating

a recent experiment on dependence of boron diffusion on silicon material.

Van Oostrum et al. [51] characterized epitaxial silicon layers grown in

different ways by monitoring the differences in boron diffusion during

surface oxidation. After oxidation-enhanced diffusion of boron spikes, they
observed remarkable differences in boron diffusion with increasing depth

position in the layers grown by three different techniques: molecular beam

epitaxy (MBE), fast gas switching chemical vapor deposition (FGCVD), and

low-temperature CVD (LTCVD). As a tentative means of testing the trap-

mediated diffusion model in this section, two of their results are simulated
with SUPREM-IV in this study. Figure 2.2 shows the simulation of the delta-

doped boron spikes as grown by MBE and after 20 minute annealing at 9000C.

The as-grown profile is a rough approximation to the SIMS data shown in

[51]. The total trap concentration CT used in the simulation is 1x1017 cm-3, a
typical value known for epitaxial layers. Figure 2.3 shows another simulation
for the SIMS boron profiles in the layer grown by FGCVD. In this case, the

simulation assumes larger trap concentration 1x1019 cm-3. Both Figure 2.2








1012 i 1 1 1 1 1
0.0 0.2 0.4 0.6 0.8 1.0 1.2
Depth (gpm)
Figure 2.2. SUPREM-IV simulation of the boron diffusion at 9000C with total
trap concentration equal to 1x1017 cm-2, emulating the data [51] for the case of
MBE-grown epi silicon.

- Boron at 0 min. - -Interstitial at 0 min.

1 019 Boron at 20 min. - Interstitial at 20 min.


S10L A -

0 15_ L L J -
a 1017
o 1014

1013 --_
1012 I

0.0 0.2 0.4 0.6 0.8 1.0 1.2
Depth (glm)
Figure 2.3. SUPREM-IV simulation of the boron diffusion at 9000C with total
trap concentration equal to 1x1019 cm-2, emulating the data [51] for the case of
FGCVD-grown epi silicon.

and Figure 2.3 are in good agreement with the SIMS data shown in [51]. The
simulated interstitial distributions at 20 minutes explain the difference in

depth-dependent boron diffusion for the two cases. In the MBE-grown
material, the interstitials injected at the surface diffuse almost uninterrupted
by the small number of traps in the bulk, thereby causing the enhanced
diffusion even in 1.0 pm depth position. In contrast, the decrease in boron
diffusion enhancement at deeper region in Figure 2.3 is due to the reinforced
capturing of interstitials by the bulk traps in large amount, which is probably
the characteristics of FGCVD-grown epilayers. Although the assumed value

of CT (1x1019 cm-3) may be considered too large, the quenching effect of the
OED in FGCVD material can be attributed to possible existence of impurities
from the reactor contamination, as Van Oostrum et al. [51] suggested.

Another possible reason for the reduced boron diffusion in the bulk is the
oversaturation of vacancies acting as interstitial traps [51]. Non-equilibrium
states of vacancies work as the trapping source of interstitials by
recombination, in the similar way to the bulk traps. A simulation with initial
vacancy concentration larger than the equilibrium concentration showed the
same result as in Figure 2.3. The tentative results consistent with the data
initially attest the validity of the model, which will be fully confirmed in the
next section.

2.4 Simulation of Damage-Enhanced Diffusion and OED

The above model of dopant-defect pairing was implemented into
SUPREM-IV, and the optimal estimates of the point defect parameters were

found by simulations. The optimization was performed to obtain the global
fit between the simulations and the data from both experiments on silicon
implant damage-enhanced diffusion of boron and phosphorus and on OED.

Table 2.1 shows the parameter values extracted from the modeling in this

work. Many of the parameters were newly determined by fitting the

simulations to the experimental data, which will be described in this section.

Table 2.1 Parameters extracted from the simulations of implantation damage
and oxidation effects on phosphorus and boron diffusion.

Parameters Pre-exponential Activation

CT (Float-zone)
CT (Czochralski)
DAI (w/ neutral I)

600.0 cm2/sec
5.0 x 1022 cm-3
13.0 cm2/sec
2.31 x 1021 cm-3
1.95 x 102 cm/sec
7.60 x 109 cm/sec
1.11 x 103 cm/sec
2.92 x 1015 cm/sec
2.10 x 106
1.79 x 10-7
8.16 x 10-4 cm3/sec
1.0 x 1017 cm-3
1.0 x 1018 cm-3
2.39 x 105
8.16 x 10-4 cm3/sec
1.77 x 10-5 cm3/sec
1.7 cm2/sec
1.52 eV
0.9 eV
3.85 cm2/sec
1.49 eV
0.6 eV

- 1.91

1.57 eV
2.94 eV
3.92 eV
3.56 eV

3.66 eV
0.05 eV

__________________________ ________________________________ I __________________________

All the damage-enhanced diffusion data used in this work are obtained
from experiments that satisfy the critical assumptions required for the pairing
theory; (i) the silicon damage implant dose is low enough that neither Si

lattice structure is amorphized nor any dislocations are created; also the
concentration of defects is much less than that of silicon lattice sites Cs so that
we can apply dilute concentration approximation; (ii) the silicon is
intrinsically doped to avoid high concentration diffusion effects, such as
precipitation and electric field effects on diffusivity. Only under the two
conditions is it possible to obtain physically meaningful values for the
binding energy EbAx.

2.4.1 Initial Distribution of Point Defects

Profiles of the point defects created during implantation are necessary as

an initial condition for SUPREM-IV, which numerically solves the dopant
and defect diffusion equations (Eqs. 2-1, 2-23, and 2-24) in the subsequent
diffusion steps. The program used for the initial defect calculation in this
work is a Monte Carlo model [52] implemented in SUPREM-III [53]. To

smooth the noisy defect profiles produced by the Monte Carlo method, the
Pearson IV distribution function was used with parameters determined by an
optimizer. Figure 2.4 shows one of the as-implanted interstitial profiles used
as an initial condition for the simulation of phosphorus diffusion in this
work. For the interstitial distribution, it was assumed that all the incident
silicon atoms serve as self-interstitials in addition to recoils. The Monte Carlo
simulations show that vacancies are created in the same amount as
interstitials, with its distribution located nearer to the surface. The separation

between the peak positions of CI and CV is larger with higher implantation

-- From Monte Carlo simulation
- -- Pearson IV fitting

1 022 I I
MC simulation input
E 0Si implant dose= 1x1014 cm2
S10 20 energy =- 60 keV
S- # of trajectories 150,000
| io18 Tilt angle -7
Temperature = 35C

8 116 1 Pearson IV fitting
16 -2
C dose 1.4132x10 cm
I\ R = 0.061 micron
(D 14 '
S10 St.dev = 0.039 micron
Skewness = 0.919
SKurtosis 4.008
0 0.2 0.4 0.6 0.8
Depth (micron)

Figure 2.4. As-implanted initial distribution of interstitials after the 28Si
implantation. Monte Carlo simulation from SUPREM-III is compared with
Pearson IV fitting.

Throughout the simulation in this work, it was found that the

distributions of point defects used as initial conditions make a significant

difference in the absolute amount of resultant diffusion. In all the

simulations, it was assumed that in the background below the implanted

region, the point defects exist in constant equilibrium concentrations,

although it has not been confirmed experimentally. When the initial

background concentration of vacancies was assumed to be larger than the

equilibrium concentration by an order of magnitude, the simulations led to

more than 20 % difference in final junction depths of boron, depending on

the depth location of boron profiles. Since the absolute amount of diffusion

enhancement is sensitively determined by the initial defect distributions, the

defect parameters extracted in this work have almost the same error range as
the initial defect simulation. However, the relative values of the parameters
are not much dependent on the simulations of initial point defects. In
general, the damage-enhanced diffusion predicted by any point defect-based
model is contingent on the simulation of initial as-implanted point defect
distribution. More accurate diffusion simulation requires conclusively
quantitative knowledge on the implantation damage based on possibly direct
observations through experiments.

2.4.2 Modeling of the Enhanced Diffusion of Boron

Packan and Plummer [23, 24] performed experiments which investigate

the enhanced diffusion of boron due to the low-dose 29Si implants and its
dependence on the damage implantation dose, energy, annealing time, and
temperature. They observed that the diffusion enhancement is the largest at
the lowest temperature and that its dependence on implant dose and energy
is nonlinear. Four sets of their representative data are simulated and fitted to
verify the pairing theory model in this study. The inert diffusivity of boron
DB and its fraction due to boron-interstitial mechanism fBI used for the
simulation are from the surface oxidation experiments of Packan and
Plummer [54] (see Table 2.1). Those values were also used for the
representation of their damage enhancement data [24]. Other important
defect parameters in Table 2.1 are consistently determined through this
modeling work.
Figure 2.5 shows the data and the simulation for dependence of the
boron profile movement 4 Deff t on anneal temperature and 29Si implant
energy. The time conditions are 60 minutes at 8000C, 5 minutes at 9000C, and
2 minutes at 10000C, within which the overall transient diffusion is supposed

Boron, 29Si dose=lxl 014 -2
0 50 100 150 200 250
100 -I, I, 1 .

75- -
I .-'.- --
50 -

a 25
800C Data --- 8000C Simulation

0 I- 900C Data - -- 9000C Simulation

[ 1000C Data 1000C Simulation
-25 .. '.. . .. ... .
0 50 100 150 200 250
Energy (keV)

Figure 2.5. Total boron profile motion 5 Deff t due to 29Si implants of dose
1x1014 cm-2 as a function of implant energy and anneal temperature. The
data [24] and the SUPREM-IV simulations are compared. The anneal times
are 60 minutes at 8000C, 5 minutes at 9000C, and 2 minutes at 10000C.

to occur and finish at each temperature. The results of the simulation with

the finally extracted parameters in Table 2.1 match the data very well. It is

shown that the damage enhancement of diffusion is larger at lower

temperatures. Moreover, for the implants greater than 40 keV, the energy

dependence is stronger at lower temperatures. The location of the boron

profile is deeper in the substrate than as-implanted defect profiles of 200 keV

energy condition in this experiment [24]. Therefore, it implies that the rapidly

diffusing interstitials play the critical role in determining the transient boron

diffusion. Also, the equilibrium interstitial concentration C*, which is lower

at lower temperature, will affect the duration of the transient diffusion of the

dopant-interstitial pairs. DI and C* are the most important parameters to be

decided along with the binding energy. In this boron experiment, float-zone
substrate was used [23, 24], and the total trap concentration in that material is
assumed to be so small (= 1.0x1017 cm-3) that the trap effect on resultant
dopant movement is negligible.
There are two sources for DI and C* in the literature. One is estimated
from the studies on gettering of gold by Bronner and Plummer [29], and the
other from gold diffusion investigation through rapid optical annealing by
Boit et al. [30]. SUPREM-IV can be used to simulate the major data of OED
effects successfully with both values of DI [28]. For the damage modeling in
this study, a large amount of fitting was performed by using Boit's values of
DI and C\ first, but it was not possible to obtain an acceptable global fitting
that matches both the OED and the damage data. With Bronner's DI and C\
shown in Table 2.1, however, correct simulation results were obtained,
especially with regard to the temperature dependence of damage
enhancement. Since the maximum enhanced diffusivity is limited by
(KAI C*") -1 as shown in Eq. 2-12, the difference between the effective binding
energy EbAI and the activation energy of Cl will determine the temperature

dependence shown in the boron data. The time-averaged effective diffusivity
and junction depth 4 Deff t will be roughly proportional to the maximum
instantaneous diffusivity enhancement for the chosen time conditions of
these data. The difference between the two energies should be negative and
its magnitude decides the degree of the temperature dependence, as Eq. 2-14
demonstrated. The EbBI for boron-interstitial pair is found to be around 1.52
eV from general fitting of all the Packan's data. With the Boit's value 1.58 eV
for the Ct activation energy [30], the simulation results in reversed
temperature dependence to the data in Figure 2.5, due to such a small energy
difference. Rather, correct dependence was obtained from the Bronner's

larger activation energy for C1 (2.36 eV in Table 2.1). Since the product DI Ci is
experimentally established to follow an Arrhenius relationship [55], DI and

CI compensate for the effects of each other on damage-induced diffusion.
The surface recombination velocity of interstitial Kimin is also critical in
simulating correct temperature dependence. The vacancy surface recombina-
tion velocity Kvmin makes a larger difference in the OED simulation and
therefore has been determined from the OED fitting. The fit in Figure 2.5 is

possible only by increasing the activation energy of Kimin so that larger
concentration of interstitials can be reduced near the inert surface at high
temperature. Surprisingly, it was observed that Kimin considerably affects the
diffusion enhancement not only around the surface but also in deep regions
when anneal time is large. This is probably because the boundary condition
given by Eqs. 2-25 and 2-26 becomes more sensitively dependent on KImin as
defect distribution approaches the equilibrium, and even a small amount of
remaining excess CI near equilibrium in all the region may cause additional
diffusion of dopant.

The dependence of boron diffusion on silicon implant dose is shown in
Figure 2.6 at various temperatures. From the data it was observed that
doubling the implant dose did not double the effective diffusivity [23]. As the
maximum diffusivity is limited by temperature, the only factor that makes

difference in Deff t for each dose condition is the duration of defect
supersaturation above the threshold shown in deriving Eq. 2-12. In this
model, the saturation time is determined mainly by bulk recombination, and
it can not be linearly proportional to defect amount; it can be shown that the

total defect concentration during the initial transient time period is a

cotangent hyperbolic function of time, by solving the defect equations


Boron, 29Si implant at 180 keV
1 5 0 , ,, ,, , ,,,1 , ,, ,i . .
O 8000C Data --- 8000C Simulation
125 A 9000C Data - A- -900C Simulation

] 10000C Data 1 -1000C Simulation
O 100-

C 75


2 5 . .. I . I I .
1011 1012 1013 1014 1015
Dose (cm-2)

Figure 2.6. Total boron profile motion q Deff t due to 29Si implants at 180
keV as a function of damage dose and anneal temperature. The data [23, 24]
and the SUPREM-IV simulations are compared. The anneal times are 30
minutes at 8000C, 10 minutes at 9000C, and 4 minutes at 10000C.

(Eqs. 2-23 and 2-24) under the assumption that the bulk recombination is the
only dominant process during the short time [48].
The simulation for 1x1012 cm-2 dose condition at 10000C in Figure 2.6
show slightly smaller enhancement than the data, where the annealing time
is 4 minutes rather than 2 minutes for the data of Figure 2.5. It is probable

that the transient diffusion has already finished within less than 2 minutes at

10000C. The junction depth difference in simulation is within the error range
of inert diffusivity of boron DB; calculation with the Di used in the

simulation shows that the junction depth increase due to intrinsic diffusion
for the 2 minutes is approximately 13 nm. The experiments are less sensitive

to the implantation damage amount at higher temperature and with less


Boron, 8000C, 29Si implant at 200 keV
1 7 5 , . .. 1 . . 1 'I -,
0 5x1013 cm-2 Data -- 5x1013 cm2 Simu.

150 1x1014 cm2 Data --I-- x1014cm2 Simu.

S ~ 2x1014cm2 Data - 2x1014cm2 Simu.

C100 ----- -------


5 0 I I I I I
0 25 50 75 100 125 150
Time (min)

Figure 2.7. Anneal time dependence of total boron profile motion / Def t
due to 29Si implants at 200 keV with various doses. The data [23, 24] and the
SUPREM-IV simulations are compared.

damage amount, because the inert diffusion tends to dominate the junction


The time dependence of the damage-enhanced diffusion at 8000C is

shown in Figures 2.7 and 2.8 under different dose and energy conditions,

respectively. The simulation superbly matches the data in Figure 2.7. In

Figure 2.8, the junction depth motion in simulation is 15 20 nm smaller

than data for anneal times less than 20 minutes. However, the important

point in the data shown in Figures 2.7 and 2.8 is that the enhanced diffusivity

is almost independent of both dose and energy for times shorter than 20 25

minutes. When anneal time is below the time constant of transient

diffusion, the boron diffusivity is determined thermally with the dopant-

defect binding energy. With larger dose or energy, it takes more time for the

Boron, 8000C, 29Si dose=lx104 cm2
1 5 0 I . . 1 . . I . . . .
0 40keV Data ---40 keV Simulaton

125 80 keV Data - --80 keV Simulation
E 200 keV data -200 keV Simulation



25 I
0 10 20 30 40 50 60 70 80
Time (min)

Figure 2.8. Anneal time dependence of total boron profile motion 1 Deff t
due to 29Si implants of dose 1x1014 cm-2 with different energy conditions.
The data [24] and the SUPREM-IV simulations are compared.

defect concentrations to approach their equilibrium value. Therefore, the
time constant becomes larger but not linearly, and the enhancement lasts
longer as shown by the data. The SUPREM-IV simulation based on the
pairing theory correctly models this behavior of transient diffusion in both
Figures 2.7 and 2.8.
The diffusion during the initial short time period seems to be effectively
modeled by the bulk recombination rate KR. However, its value is limited by
DI and Dv [56], and is controllable only by introducing the energy barrier for I-
V recombination AEiv. There exist discrepancies among different researchers
with regard to the temperature dependence of the I-V recombination [20, 57,
58]. Moreover, dopant-mediated recombination and the possible effects of
charged defects may complicate meaningful extraction of the parameter [28].

Due to such uncertainty of this parameter, the extraction of all the other
parameters was performed first, so that they do not rely as strongly on KR. It
was found that the assumed value 0.75 eV for AEIv gives correct dose and
time dependence of Figure 2.7. For both OED and implant damage modeling,
that value of AEIv generated generally correct results, even though it could

not be precisely determined. Accurate assessment of KR will be possible only
by monitoring the defect recombination for an extremely short initial period
of diffusion, which the present data do not directly represent.

With consistent parameters, the effective binding energies for boron-
interstitial and boron-vacancy pairs were extracted to be 1.52 eV and 0.9 eV,
respectively. Since fBI is about 0.8, the simulation is much more sensitive to
interstitial parameters than to vacancy parameters for most anneal time
conditions. Accordingly, the error range for EbBI is less than that for EbBv. It
was observed that the dopant distribution becomes obviously non-Gaussian
with values of EbBv larger than 1.0 eV. According to the physics-based
argument in section 2.2.2, the EbBV should be smaller than the activation
energy of the equilibrium concentration of vacancies C which is extracted to

be 1.08 eV in this work. Thus the simulation definitely shows that the model
parameters are extracted consistently with the prediction from energetic.
The C, and Dv were mainly determined from the OED fitting, since they were
found to be the most crucial parameters for the OED modeling in two

2.4.3 Modeling of the Enhanced Diffusion of Phosphorus

The data for silicon implant damage effects on phosphorus were
obtained recently by Park and Law [25]. The experimental procedure is
consistent with the assumptions of pairing model, and searches for the

Phosphorus, 28Si dose=1xl014 -2 at 60 keV
1 0 0 , 1, , 1 1 1 , I 1, 1 l
O 8000C Data -- -8000C Simu.
80 A 900C Data -A--9000C Simu.
S. [] 1000C Data 1 -1000C Simu. /
o 60 /
1100C Data 1100C Simu.


20 .

0.1 1 10 100
Time (min)

Figure 2.9. Phosphorus profile motion /Defft due to 28Si implants of dose
1x1014 cm-2 at 60 keV as a function of anneal time and temperature. The data
[25] and the SUPREM-IV simulations are compared.

dependence of transient diffusion on anneal temperature and time, including

RTA conditions. The data and the simulation which led to the extracted

parameters of Table 2.1 are shown in Figure 2.9. Since the anneal times span

15 seconds to 60 minutes, careful interpretation of the results is needed. It

should be regarded that the time constant for the transient diffusion steeply

decreases as temperature increases. The junction depth from the damage

enhancement at 8000C is larger than at 9000C for all time conditions,

consistent with the previous boron diffusion results. Comparing the

junction depths between the damage-enhanced diffusion and undamaged

inert diffusion shows that the time constant at 10000C is less than 15 seconds

[25]. The time dependence of phosphorus diffusion is not necessarily the

same as for boron for several reasons. First, the boron experiments of Packan

and Plummer [23, 24] had a spatial displacement between the damage profile
and the boron profile, and there was no such displacement in the phosphorus
experiment of Park and Law [25]. Second, phosphorus has a larger interstitial
fraction than boron. Because of these two reasons, the vacancy transient
tends to be longer in the phosphorus case; since larger amounts of
interstitials are paired with phosphorus atoms, more vacancies remain in the
shallow region and are governed more by surface recombination. Third,
Packan and Plummer [23, 24] used float-zone material and Park and Law [25]

used Czochralski. The different number of traps contributes to different time
The three data at 10000C for the times larger than 15 minutes obviously
show inert intrinsic diffusion. There is a spread in the reported values of
intrinsic diffusivity of phosphorus [54, 59 ~ 65], and therefore the deviation in
the results could be due to the uncertainty in phosphorus intrinsic diffusivity
at high temperatures. However, Figure 2.9 shows that the pairing model
simulates the overall time and temperature dependence reasonably well.
Without considering the dopant-defect pair contribution, the diffusivity
enhancement would be proportional to the excessive defect supersaturation
from the beginning, which would cause extremely large junction depths at all
the temperatures.
The optimal phosphorus-interstitial binding energy Ebpj was found to be
1.49 eV. This compares well with the theoretical lower limit of Ebpi of 1.4 eV

[20]. Since phosphorus diffuses mainly by an interstitial mechanism, the Ebpi
critically determines the simulation results. The phosphorus-vacancy
binding energy Ebpv was extracted to be approximately 0.6 eV. The Ebpv has an
upper limit of 0.8 eV above which most simulated profiles take a non-
Gaussian shape. An example of this is shown in Figure 2.10 where Ebpv is

Distribution of phosphorus and defects
(Eb pl1.49eV; b P=1.2eV) (at 8000C)








I I I I I I. .
0 0.1 0.2 0.3 0.4 0.5
Depth (micron)
Figure 2.10. An example of non-Gaussian dopant diffusion shown in
SUPREM-IV simulation; for the sake of illustration, Ebpv is increased up to
1.2eV while Ebpi is kept as the extracted value, 1.49 eV.

increased up to 1.2 eV with other parameters unchanged. As can be seen in

the solid line, there is a substantial increase in the tail diffusion of

phosphorus. It can be attributed to the steep gradient of vacancy distribution

and the higher binding energy, which locally affects the diffusive flux of

phosphorus-vacancy pairs during the transient diffusion period.

Since Czochralski material was used for the substrate in this

phosphorus experiment, a larger total trap concentration of 1.0x1018 cm-3 was

assumed in contrast to the previous boron case. Therefore, the reaction of

defects with traps has more impact on the simulation, particularly on the

time dependence characteristics. Two parameters in the trap implementation

were found to be important in modeling the phosphorus enhanced diffusion.

First, the total trap concentration CT is the measure of total influence of traps

-- -Phosphorus t=O min
\ Phosphorus t=1 min
----- Interstitial t=1 min
\ - -Vacancy t=1 min



-- -
- -- ,

on interstitial transport. At the initial stage of transient diffusion, the

excessive amount of interstitials quickly fill most of the empty traps
determined by CT Increasing CT will reduce the amount of free interstitials
to be paired with dopant atoms, and will decrease the transient dopant
diffusion period. After the initial stage, the trapped interstitials will be
dissociated from the traps, which approach equilibrium, and the delay of the
diffusion enhancement at the initial stage will be compensated for. In
modeling 9000C data in Figure 2.9, it was observed that increasing CT up to
3.0x1018 cm-3 reduces the junction depth at 15 and 30 seconds, while at 60
minutes the results are the same as when CT is 1.0x1018 cm-3.
The second important factor in trap implementation is the vacancy-
filled trap reaction, which was necessarily introduced to improve the time-
dependence modeling. As shown previously in this chapter, the dual
reaction of traps with both interstitial and vacancy is physically more

reasonable. The reaction rate KtrapV is diffusion-limited with an energy
barrier AEFT-V which in general will be also dependent upon charge states.
For the phosphorus simulation, an optimal value of 1.0 eV was used for the
AEFT-v. The effects of the vacancy reaction with filled traps can be observed
usually in large time annealing data. The forward reaction of vacancy and
filled trap, i. e., trap-mediated recombination will increase the empty trap
amount and also annihilate the interstitials that were paired with the traps.
At the initial period of excessive defect supersaturation, its effect on CI is
insignificant because the interstitial-trap reaction rate KtrapI is larger than

KtrapV. However, when transient diffusion is almost finished and CI goes
near equilibrium, the accumulated empty traps will start to effectively interact
with the small amount of interstitials, decreasing the CI down to C* much
more rapidly. This is shown in Figure 2.11, where the time variation of the

Time variation of C / C* at z=0.1 pm; 9000C; CT=1xl 018 c3
106 I

5 A EFT-V = 1.9eV
- --- EFTV = 1.0 eV
104 - Without traps





10-3 10-2 10-1 100 101 102
Time (min)

Figure 2.11. Interstitial supersaturation at 0.1 im depth position for different
AEFT-V values in the simulation.

interstitial supersaturation C / C* is simulated at one depth position.

Compared to the case without traps, the trap-mediated recombination results

in strong effects on the defect diffusion near equilibrium for a long time.

Hence the vacancy reaction with traps leads to reduced phosphorus diffusion

at large times. When the AEFr-v is reduced below 1.0 eV, the junction depth

of the 9000C, 60 minute simulation profile decreases, whereas the simulation

at 15 seconds does not change.

To refine the results at 9000C, temperature-dependent inert diffusivity

fraction fpi was used for the phosphorus modeling. Its values from the

Arrhenius relationship in Table 2.1 are 0.95 at 11000C and 0.884 at 9000C. This

range does not contradict the minimum value 0.94 at 11000C that was

extracted from nitridation studies [66]. At lower temperatures, the diffusivity

of phosphorus-interstitial pair DP*I was kept as the default from Fair [65], and
vacancy contribution Dp* was added to match the short-time RTA data at

9000C. However, the fpI adopted in this study can not be considered to be
conclusive, in so far as such minute variation of the fpI does not critically
improve the simulated results. In the previous section, the boron simulation
used the value of fBI (= 0.8) extracted by Packan and Plummer [54] based on
their value of fpI from oxidation experiment. With regard to significance of
the difference, however, it was observed that with fBI of about 0.7, the boron
simulation results do not change beyond the experimental error range, which
might be compensated for by changing the boron-interstitial binding energy

within its error range (about 0.02 eV). Therefore, all the extracted
parameters are self-consistent in modeling both phosphorus and boron

2.4.4 Modeling of the Oxidation Enhanced Diffusion

Although fits for most of the intrinsic OED data have been presented

[28], simultaneous fitting to the damage data required a substantially larger
value of Kimin. It is easy to obtain fits to the OED of phosphorus by adjusting
the fraction of consumed silicon injected as interstitials, 0. The larger value
of Kimin makes it difficult to fit the measured interstitial lateral decay length,
which required changes in the values of the vacancy equilibrium
concentration, C*, and the bulk recombination constant KR. The vacancy
diffusivity is also affected since the vacancy equilibrium concentration-
diffusivity product is held fixed to the value estimated by Tan and G6sele [55].

These parameters were the keys to simultaneously fitting both the OED and
silicon damage effects. Since a substantially smaller value of KR is used, there
is less interstitial recombination with vacancies in the near surface region.

This compensates for larger surface recombination by decreasing the amount

of bulk recombination flux for the interstitials. This can be further localized

to a surface effect by adjusting the value of Kvmin which controls the rate at

which vacancies are replenished at the surface of the wafer. The higher value

of C* and the smaller value of Dv result in a shorter time to achieve

interstitial-vacancy equilibrium, which affects the time dependence of the

near-surface bulk recombination. Table 2.1 shows the final parameters which

are used to model both the silicon implantation damage and the oxidation

effects on diffusion.

The experimental result of Griffin and Plummer [49] is shown in Figure

2.12 along with the simulation fit. The plot shows the lateral decay length of

the interstitial supersaturation under an inert pad oxide as a function of time

I I . .i

11000C 0

0* 10000C


* Data



* ^^0

S I TI 1 I
Time in Minutes


Sr 1 04

Figure 2.12. Lateral decay length of interstitial diffusion as a function of time
and temperature from the measured data [49] and the SUPREM-IV







i-1 r

I I I I .


Interstitial Supersaturation @1100C
5 1 Packan
20A Initial oxide M Ahn
4.5 -
4 o Ahn 201pm membrane
A Ahn 55pm membrane
3.5 -- Simulation

S400A initial oxide


1.5 backside enhancement/ 0
1 -
1 10 102 103 104
Time (minutes)
Figure 2.13. Supersaturation of interstitial in depth direction under a growing
oxide layer as a function of oxidation time. The simulation is compared with
the data extracted from the lightly doped phosphorus diffusion experiments

and temperature. This experiment is especially sensitive to the recombina-

tion along the surface of the wafer. The effects of Kvmin and the increased CV

are evident here as they provide the same lateral decay length as obtained

from earlier fits with smaller values of Kimin [28]. Figure 2.13 shows the

supersaturation of interstitial in depth direction under a growing oxide layer

as a function of oxidation time. The measurements were all from the lightly

doped phosphorus diffusion enhancement. The frontside enhancement

results are for two cases of different initial oxide thickness, matching the data

from several experimental investigations [54, 67, 68]. The backside enhance-

ment data are from the OED experiment with silicon membranes of two

different widths at the backside of wafer by Ahn et al. [68] The simulation

shows reasonable match with the data, including the onset time of

enhancement. The results shown in Figures 2.12 and 2.13 indicate that the
high interstitial diffusivity coupled with bulk traps and vacancies is effective

in modeling OED results in both the vertical and lateral directions [28, 49].
Griffin's extracted diffusivity is several orders of magnitude below that used
in these simulations, but the time dependence of the interstitial decay length
in Figure 2.12 is still correct, which demonstrates the effectiveness of the trap
Although the activation energies of the parameters in Table 2.1 and in
Law [28] are different, the values of the parameters were adjusted in
complimentary ways to maintain a good fit with the OED data. The bulk

recombination was decreased from the value in Law [28] and to compensate
for this, the equilibrium vacancy concentration was increased. Since a larger
value of KImin was required to fit the damage data, this was compensated by
less interstitial-vacancy recombination in the near surface region by
decreasing Kvmin. These adjustments maintained a good fit to the interstitial
decay length as shown in Figure 2.12, but substantially changed the amount of
vacancy undersaturation which is shown as a function of time at 11000C in
Figure 2.14. The vacancy surface recombination replenishes the surface
concentration of vacancies, and the resulting profile has a fairly strong depth

dependence. The smaller value of Kvmin and increased C* do not perturb the
simulation far from the data. In fact, the data of Packan are upper bounds, so
the fit is not unacceptable at short times, and agrees well at longer times

considering the wide spread in the data.

2.5 Summary

Transient dopant diffusion in intrinsic doping conditions is simulated
correctly by point defect models based on the theory of dopant-defect pairing.

Vacancy Undersaturation @1100C Oxidation
-0 .5 I I.' ' . ..- ' 1 . -11
-- Simulation + Mizuo
SA Packan A Antoniadis

-1.5 A

-2 A
0* m


-3 .5 I I.. .... I . . I . ..... I . ..
1 10 102 103 104
Time in Minutes

Figure 2.14. Vacancy undersaturation near a growing oxide layer as a function
of time at 11000C from several researchers and SUPREM-IV simulations. The
data are from ref. [57, 58, 69]; P. Packan's data are unpublished.

The model parameters are extracted to consistently fit both the experiments of

low dose silicon implant damage effects and of oxidation enhanced diffusion.

This supports the assumption of local equilibrium in SUPREM-IV diffusion

model. The implantation damage effects on boron and phosphorus diffusion

can be characterized through the pairing coefficients determined by the

effective binding energies. Interpretation of the pairing theory in terms of

diffusivity enhancement clarifies the thermal nature of the transient

diffusion under excessive point defect supersaturation, accounting for

observations from the measured data. The thermally-determined quantity

KAI CI is found to quantify the enhanced diffusion of intrinsically-doped

boron and phosphorus due to implant damage. The intrinsic doping


condition validates the simplification of charge state effects in the models
shown in this work. The extracted parameters provide a substantive basis for
further modeling of diffusion under extended conditions, such as the

extended defects model presented in the next chapter.


High dose ion implantation is an essential technique for obtaining
heavily-doped regions in silicon such as the source and drain of MOSFETs
and DRAM cells as well as the emitter of bipolar devices. The high energy
bombardment of incident ions naturally disrupts the crystalline structure of
silicon and introduces an excessive number of point defects which work as

nucleation source for the extended defects. As we saw in the previous

chapter, ion implantation damage almost predominates the subsequent
dopant diffusion during the thermal annealing cycle required for substrate
recrystallization and dopant activation. The damage effects become much

more complicated when the extended defects coexist with the point defects,
interacting with each other through a spatial and temporal variation.

The implants of common dopants at a dose above a certain ion mass-
dependent threshold amorphize the surface region in silicon substrate [13],
simultaneously producing a large amount of point defects. The subsequent
annealing leads to solid phase epitaxial regrowth of the amorphous region,
and extended defects such as dislocation loops are formed below the
amorphous-crystalline interface. The dislocation loops are inherently

accompanied by a stress field in the crystal, interacting with the point defects.

It is generally accepted [19, 70 74] that the end-of-range dislocation loops

affect the distribution of point defects by absorption of interstitials or by

emission of vacancies at their core boundary during growth, and by the
reverse processes during shrinkage. There has been a lot of effort to model

the dopant diffusion by investigating the interaction of dopant and point
defects without any extended defects under low dose implant damage and
oxidation/nitridation conditions. This chapter is focused on modeling the
evolution of dislocation loops and its effects on the point defect diffusion, and
its resultant impact on the dopant redistribution during oxidation.

3.1 Modeling of the Dislocation Loops as a Group

Previous work [75 ~ 78] established theoretical models for a single
circular dislocation loop and its interaction with point defects. Bullough et al.
[76] studied the migration of an interstitial impurity atom around a single
dislocation loop on the basis of the stress field from the loop solved by
Bastecka and Kroupa [75]. Borucki [78] proposed a model for the growth and
shrinkage of a single dislocation loop due to the capture and emission of

point defects, and simulated the point defect variation from an assumed
initial high supersaturation around a periodic array of the loops in a three
dimensional numerical solver of diffusion equations. He built his model for
the single dislocation loop growth based on the pressure field solved by
Bastecka and Kroupa. However, it is necessary to model the effects from the
group of dislocation loops formed in the substrate as observed through TEM
pictures. TEM measurements [79, 80, 81] show that the variations in
distribution and size of the actual dislocation loops during oxidation or
annealing are generally not so homogeneous and simple as in the case of one
single loop. The dislocation loops usually form a network by merging with
each other during oxidation. Coalescence and dissolution of dislocation loops
are a statistically complicated process, which also depends on the implanted
ion species [79, 80, 81]. Since the TEM measurements can only give statistical
data on density and size of the dislocation loops, a useful and correct model

for the ensemble of dislocation loops should be made in a way to reflect the
statistical data from the TEM pictures.
This work is to build up a point defect based model for the evolution of
the group of dislocation loops with both accuracy and efficiency, so it can be
implemented in two dimensional process simulators. It involves two
dimensional conversion of Borucki's model for the interaction of a single
loop and point defects in three dimensions, and it is extended to the ensemble
of the loops. In addition, we developed an efficient model for the network
formation of the dislocation loops via a statistical density distribution
function of loop radius. The evolution of the loops and the redistribution of
point defects during oxidation are simulated where the surface injection,
recombination, and diffusion of point defects have been implemented
through a point defect parameter set (Table 2.1) consistent with the current
available data on oxidation enhanced diffusion (OED). Finally, the interstitial
capture efficiency of the dislocation loops is attested by monitoring the OED of
boron in a buried layer located below the dislocation loops.
The model for the group of dislocation loops is developed on the basis
of thermodynamics and linear elasticity which govern a single dislocation
loop. Major assumptions in the model are: (i) dislocation loops are all
circular and evenly distributed on a plane interconnecting their centers; (ii)
orientation of the loops is periodic in two perpendicular directions,
approximating the morphology shown in cross-section TEM and plan-view
TEM pictures; (iii) radius and density of the loops follow an asymmetric
triangular distribution function; (iv) formulation of pressure field from a
single dislocation loop in linear elastic material can be used as a basis in
calculating the effective pressure from the ensemble of dislocation loops; (v)


ro (1+e)



Figure 3.1. Coordinate system for calculating the pressure from a single
circular dislocation loop; b denotes the direction of Burgers vector of the loop,
and e is the dilatation factor of the impurity atom. Excerpted from ref. [76].

local equilibrium is attained around the dislocation loops through their
reaction with point defects.

3.1.1 The Effective Pressure from the Dislocation Loops of Equal Size

In order to accomplish a physics-based model with both accuracy and
efficiency in two dimensions, it is very important to account for the pressure
field around the loops. The pressure from a single circular prismatic
dislocation loop has been formulated by Bastecka and Kroupa [75] from the
diagonal components of the stress tensor in linear elasticity theory. Figure 3.1
shows the geometry used in calculating the pressure from a single loop in
terms of location of an impurity atom in cylindrical coordinates. The loop
has the Burgers vector b in the axial direction. The dilatation factor e of an

1.5x10 '
Radius R=200A

"E 5.0x109

2. o

)-5.0xl 0

-1.0x10 _

-1.5xl 010
-800 -400 0 400 800
Radial distance r (A)
Figure 3.2. The pressure from a single dislocation loop calculated by Eq. 3-1 at
different relative positions from the loop center.

atom near the loop measures the space for an elastic inclusion of the atom.
The pressure is expressed in terms of the loop radius R, the radial distance r,
and the height h from the loop [75]:

b R [R2-_ r2- h2 1
p [ y ( R r ) + h2 E() +K3-1

where a = { 4 r R / [ (r + R )2 + h2] }0.5, E(a) and K(a) are the complete elliptic
integrals of the first and the second kind, b is the magnitude of the Burgers
vector, and p is the shear modulus, and y = 3 ( 1 v) / (1 + v ), where v is the
Poisson's ratio.
Figure 3.2 shows an example of the calculated pressure at different
positions of the atom as a function of the radial distance and the height. In
general, the pressure is positive inside the cylindrical region of the loop

(up to surface)

(into bulk substrate)

Figure 3.3. The cross-section viewed in a two-dimensional process simulator
(shaded plane) is perpendicular to the layer of dislocation loops (linkage of
dotted lines).

radius R, and negative outside. Its magnitude decreases approximately in the
inverse proportion of the cube of the distance in the outside region, if the
distance from the center exceeds the loop diameter [76].
In two dimensional process simulators, the cross-section where dopant

and point defect concentrations are calculated is perpendicular to the layer of
dislocation loops formed inside the substrate in parallel with the surface. A
simplified diagram in Figure 3.3 helps us visualize the dislocation loop layer
center plane, which is assumed to link the centers of all the circular loops.
The loop layer formed by ion implantation is usually parallel to the wafer
surface, except at the mask edge. Therefore, it is necessary to obtain an
effective pressure from the group of dislocation loops at one depth position by

cross-section of a dislocation loop
L ------------------


t L

, /

R i

( t
,-------- - -
(0,0): x [D]

---------- ----------- L





i I


SPloop(i) dx dy

Figure 3.4. The assumed configuration of dislocation loop ensemble viewed
from the top of substrate; the shaded rectangle is a symmetrical integration
area for calculating the average pressure from the six nearest loops.

considering a certain configuration of the loop ensemble. The configuration

of the ensemble of equal-sized loops assumed in this model is shown in

Figure 3.4 which is viewed from the top of the substrate. We assume that the

loops of one radius are evenly distributed on one plane interconnecting their

centers. The orientation of the loops are assumed to be periodic in two

perpendicular directions, with their Burgers vectors lying on the center plane.

This is a good approximation of the loop morphology as observed in plan-

view TEM and cross-section TEM pictures, where most of the category II end-

of-range dislocation loops in {100} Si substrate have their Burgers vectors in

either <110> or <111> directions [13]. The effective average pressure at one

depth position due to the ensemble of dislocation loops of identical size is

- -


by numerical integration
5.0x109 I -Pressure along (x=0,y=0) [A]
- -Pressure along (x=L,y=L) [B]
--- Pressure along (x=L,y=50A) [C]
'2.5x109- - -Pressure along (x=L-R,y=0) [D]

0- 0* /

-2.5x109 /
i I I
R = 200A
/ L = 500A
-5.0x1 0 1 1 , 1
0 200 400 600 800 1000
Distance Z from the center plane (A)
Figure 3.5. The average pressure

from the equal-sized loops compared
with the individual pressures at four different vertical locations marked as
[A], [B], [C], and [D] in Figure 3.4. The loop layer center plane corresponds to
the plane viewed in Figure 3.4.

calculated by integrating the pressure components from the six nearest

individual loops on the shaded square region. The formulation of the

pressure from a single circular loop shown in Eq. 3-1 was used for the

calculation. The pressure components from the loops farther than the six

neighbors are negligible, for the pressure from a single loop diminishes

rapidly outside.

The calculated average pressure

from the uniform circular loops is

shown in Figure 3.5 as a function of the depth distance Z from the loop layer

center plane, in comparison with the actual pressure at four different

locations [A] to [D] in Figure 3.4. The pressure in Figure 3.5 is only for the

half-side with respect to the layer center plane. It is symmetrically extended

below and above the dislocation loop layer. The radius R and the unit

interloop distance L used in the calculation are 200A and 500A, respectively.
The average pressure

is found to be compressive inside the loop layer of
thickness 2R, and to be tensile outside. This is a meaningful approximation
consistent with physical reasoning on the configuration of the dislocation
loops as observed in cross-section TEM (XTEM) pictures. The XTEM pictures
show that the loops are oriented such that they form a layer with definite
thickness roughly equal to the average diameter of the loops. The magnitude

of positive pressure inside the individual loops is larger than that of negative
pressure outside, as shown in Figure 3.2. Therefore, the average pressure

shown in Figure 3.5, calculated by superposition based on linear elasticity of
silicon substrate, approximates meaningfully and effectively the pressure
around the loop layer in two dimensions with regard to its shape and
Analytic functions for the average pressure are extracted in this work to
replace the time-consuming numerical integration of the pressure
components. The average pressure from one-sized dislocation loops can be
expressed as a function of three variables, i.e., radius R of the loops, uniform

loop density D or unit interloop distance L, and normal distance Z from the
loop layer center. The detailed subroutine including the optimized analytic
functions for the average pressure is shown in Appendix. Figure 3.6 shows
the numerical calculation of average pressure with three different pairs of R
and L. For a given L ( or density ), the magnitudes of the maximum average
pressure at Z = 0 and the minimum (negative maximum ) average pressure
at Z R increase with larger R, but it is not linearly proportional. The
pressure decreases with larger L, i.e., with lower density, as expected. As a
rough approximation, the calculated average pressure at the loop layer center
and at the layer boundary is found to be almost linearly proportional to the

3.0x109 ,- , I I I I I 1
R = 100 A; L = 500 A

S- -R = 200 A; L = 500 A
- -- R = 200 A; L = 1000A

>- 1.0x10

0 200 400 600 800 1000
Distance Z from the center plane (A)
Figure 3.6. The average pressure

from the group of uniform dislocation
loops for different radius R and unit interloop distance L.

density of loops, especially when the density is large. The analytic function

for average pressure

was found by deriving the best-fit coefficients of the

function through an optimizer program. The maximum error between the

analytic fit and the numerical solution was less than 1% in the most

interested range of R, L, and Z. Figure 3.7 shows the functional variation of

the maximum pressure at the loop layer center and the minimum pressure at

the layer boundary, as calculated from the derived analytic function. It

should be noted that the largest possible magnitude of

max or

min for a

given radius, which is the value at the end point of each line, is actually

larger with smaller radius. It implies that the pressure

max or


from the dislocation loop layer can be larger when the loop size is smaller, as

far as the density is high enough. As it will be discussed in the next sections,

the density of loops decreases as the loops grow during oxidation. The

relative values of growth and coarsening rates determine whether the







" n 4 n


max ; R=100A

max; R=200A
- - <>max R=400A


max; R=800A
:- -

min ; R=10OOA

min ; R-200A
- --

min ; R400A

min; R=800A

Z :-7 - - - -

-O.UA IV 'I I I l i I
0 400 800 1200 1600 2000
Unit interloop distance L (A)

Figure 3.7. The maximum and minimum values of the average pressure as a
function of unit interloop distance L for different loop radius R.

maximum pressure decreases or increases during a thermal process. It is

reasonable that the strain from the extended defects tends to be relaxed during

annealing without any external sources of point defects. In case of the loop

growth process such as oxidation, however, the time dependence of pressure

variation should be investigated through simulation or be confirmed directly

by measurement. In Figure 3.7, the pressure variation can be visualized by

following a certain line spanning from the left to the right as oxidation goes

on, and the slope of the line is determined by comparing growth and

coarsening rates of the dislocation loops. The simulation result on time

dependence of the pressure variation during dry oxidation will be discussed

later in this chapter.



. I t t

3.1.2 Asymmetric Triangular Density Distribution of Loop Radius

A preliminary model for the equal-sized loops was reported previously

by Park and Law [82]. However, since TEM data show non-uniform radius

and density of the dislocation loops, an advanced model should encapsulate

the distribution of the loop size via a statistical function. An asymmetric

triangular distribution function is used in this work for its simplicity and

approximate match to the TEM data. An example is shown in Figure 3.8,

where the triangular density distributions are compared with the loop radius

distribution statistically extracted from the TEM pictures [81]. The evolution

of loops was monitored through TEM measurements during dry oxidation at

9000C for 4 hours [81]. The scatter plot of loop density was produced by

Si Implant 2x1015 ac-2; 9000C Oxidation
1x1010- I II I
o 0.5 hour
0 o 1 hour
8 8x109 a hours
SA 4 hours
S6x109 -- triangular distribution

.2 4x109 -

S2x109 /A CD A
Sa AA&

0 1 / II' I I I

0 200 400 600 800 1000 1200
Loop radius (A)

Figure 3.8. The asymmetric triangular density distribution function is applied
to the statistical distribution of loop radius extracted from the TEM measure-
ments (ref. [81]) under 9000C dry oxidation condition.

counting the number of loops within each range of radius in randomly

sampled areas of the TEM micrographs, as explained in reference [81].

Extracting the histogram of loop radius distribution necessarily involves

discretization with a certain sampling period of radius. Application of the

mathematical triangular distribution is achieved first by formulating the

probability of finding dislocation loops with a certain radius without

discretization. The unnormalized probability density function fd(R) for the

distribution of loop radius R can be expressed as:

2 Dall ( R Rmin)
fd(R) = RmaxRmin (Rp Rmin) when Rmin < R < Rp 3-2
( Rmax Rmin ) ( Rp Rmin )

d(R) = 2 Dall ( Rmax R)3-3
fd(R (Rmax Rmin ) ( RmaxRp wen Rp < < max 3-3

where Rmax, Rmin, and Rp represent the maximum, the minimum, and the

majority loop radii, respectively. The total density Dali of all the loops is the

normalizing factor of fd(R), satisfying the following relationship:

Dall = fd(R) dR 3-4

In order to obtain the loop size distribution data, it is necessary to interpret

the TEM pictures by counting the number of the loops whose radius is within

a certain range AR. Thus, the discretized density distribution function D(R)

for the loops of radius R depends on the sampling period AR:

D(R) = 2 fd(R') dR' = fd(R) AR 3-5

The discretized triangular density distribution is shown in Figure 3.9. The

discretized density of loops with each range of radius is denoted as an array


o AR


Rmin Rp Rmax Radius

,4 time = tl
"o Dp
I te

D1, D2, etc. Temporal change of the distribution represents coalescence and
dissolution of the dislocation loops as observed during oxidation and
Rmin Rp Rmax Radius

Figure 3.9. The discretized asymmetric triangular density distribution func-
tion for radius of dislocation loops.

D1, D2, etc. Temporal change of the distribution represents coalescence and

dissolution of the dislocation loops as observed during oxidation and

annealing. The density Dp of the majority loops with radius Rp and its

relation with the total density Dali can be expressed by using Eqs. 3-2, 3-3, and


2 Dall AR 2 Dali
Dp = fd(Rp) AR = Rmax- Rmin m 3-6
Rmax Rmin m

where m is defined to be ( Rmax Rmin ) / AR 2 Dail / Dp. The factor m, or

the sampling period AR of the radius in the interpretation of TEM data

correlates the total density Dali with the majority loop density Dp. Thus, the

discretized loop size distribution is represented by parameters Rmin, Rmax, Rp,

Dp, and m. The Dall and m are directly determined from TEM measure-

ments. The average radius of the distributed loops is derived as:

Rmax fR-ax
Rave = C RiDi / Di = R-/d(R) dR fd(R) dR
i i Rmtn JRmin

= 3 (Rmin + Rp + Rmax ) 3-7

Therefore, the quantities such as Rave and Dall, which are measurable from

the TEM pictures by pattern recognition without discretization, can be given

by closed-form functions in terms of those distribution parameters.

Based on the linear elasticity, the total effective pressure P from the

entire distribution of dislocation loops can be calculated by superposing the

average pressure

from the loops of various size, under a simplifying but

reasonable assumption that the centers of all the loops lie on one layer center

plane. In this model, it is calculated by integrating

with respect to radius

R over the triangular distribution at each time step:

P= =--

i AR J
= 2 DRal f

dR 3-8
Dp (Rmax Rmin) R

In Figure 3.10, the pressure P from an entire distribution of dislocation loops

is shown along with the pressure

from the loops of major size ( R = Rp,

D = Dp), and also with the pressure calculated simply by assuming an

equivalent distribution of one-sized loops of the average radius, that is, by

simply substituting R = Rave and D = Dall in the analytic function

. It is

noticeable that the total effective pressure P changes from compression to



$ 2.0x109-
2 1.0x109-

2 0-


P from the entire distribution

from loops of radius Rp

calculated with R=R ave; L=L all
Rp= 200A Rmin= 50A Rmax =400A
Lp=800A LalI=500A Rave=217A

I -

-2.0x109 -1 1I 1 I 1 I I I I .
0 200 400 600 800 1000
Distance Z from the center plane (A)
Figure 3.10. The total effective pressure from the entire distribution in
comparison with the pressure component from majority of the loops
, and the pressure from an assumed distribution of all
average-sized loops , where Dp = 0.5/Lp2 and Dall =

tension around the loop layer boundary less abruptly than

( R = Rave; D =

Dall). It implies that the fuzziness in the transition region due to the various

size of the loops is effectively modeled through the superposition of pressure

following the distribution function.

3.1.3 Interaction of Dislocation Loops and Point Defects

The interaction energies between the ensemble of dislocation loops and

the point defects are expressed in terms of the total effective pressure P and

the elastic volume expansion susceptible to the external pressure effect on

interstitial and vacancy, AVI and AVV, respectively. As in the case of single

loop model [76, 78, 83], those volumes are given by assuming the sphericity of

the unrelaxed point defects with the values of radii, ro = 1.11A for the
interstitial, rs = 2.47A for the vacancy. The effective interstitial volume
expansion AVI is dependent on the dilatation factor e shown in Figure 3.1,
which determines the elastic inclusion of an interstitial [76]. The pressure
effect from the three dimensional morphology of the dislocation loops is
effectively modeled in two dimensions by using the unestablished physical
factor as a model parameter. The value of e used in the simulation is 0.5,
which is reasonably within its meaningful range ( 0 < e < 1).

The interaction between the loops and the point defects is primarily
reflected on the effective equilibrium concentrations of point defects around
the dislocation loop layer [83]:

CI*(p) = C'*(p=o) exp -kT 3-9

Cv*(p) = Cv*(P=O) exp (I ) 3-10

where P is the total effective pressure given by Eq. 3-8. The above expressions
for the pressure-dependent equilibrium point defect concentrations are

physically valid, for the formation enthalpy of an interstitial increases by the
amount of the interaction energy PAV, whereas that of a vacancy decreases by
-PAVy. Under a compressive medium as inside the dislocation loop layer,
therefore, CI*(p) decreases whereas CV*(p) increases with respect to those
nominal values in the absence of the external pressure, CI*(P=0) and CV*(P=0),
respectively. The gradient of pressure is the driving force for the point defect

movement around the dislocation loop layer.
Growth and shrinkage of the dislocation loops are modeled in terms of
their reaction with point defects at the loop layer boundaries. The model
simultaneously accounts for emission and absorption of point defects at the

dislocation loop layer boundaries in two dimensions. The layer boundaries
or the locations of reaction are assumed to be at the distance Rave from the
center of the layer. The boundary conditions are given by the reaction rates of
dislocation loops and the point defect formation energy change due to the
loop growth or shrinkage. The interstitial continuity equation in the
presence of dislocation loops can be derived by considering an additional flux
due to the local variation of the free energy as in previous work [78, 84, 85, 86],
and by incorporating the pressure-dependent equilibrium interstitial
concentration as shown in Eq. 3-9:

S= V [DI CI(p) V )] KR ( CI CV CI*(p) CV*(P))

KIL ( CI CIb) at loop layer boundaries 3-11

where KR is the bulk recombination rate, and KIL is the constant of reaction
between the interstitials and the dislocation loop ensemble. It should be
noted that the flux term now includes the pressure effects from the
dislocation loops in terms of the varied equilibrium concentration of
interstitials CI*(p) as well as the driving force for the point defect movement
around the dislocation loop layer due to the gradient of pressure. CIb is the
effective local equilibrium concentration of interstitials at the loop layer
boundary, which is modified for this two dimensional model from the
formulation of the boundary condition for a single loop [78, 87]:

CIb = gbc CI*(P) exp( -AEf / kT) 3-12

where gbc is a geometry factor (= 0.3 ) which approximately converts the three
dimensional boundary condition into two dimensions. AEf is the change in

defect formation energy due to the self-force of a dislocation loop during
emission and absorption processes at its edge [77]:

AEf=- gb n 8Rave 2V 1 3-13
47x(l-v)Rave I rc 4v 4

where p is the shear modulus, b is the magnitude of Burgers vector of the
loop, Q is the atomic volume of silicon, rc is the core (torus) radius of the
loop, v is Poisson's ratio, and Rave is the average radius of the dislocation
loop ensemble as given in Eq. 3-7. Vacancy continuity equation can also be
derived in an expression similar to Eq. 3-11, with a different reaction rate KVL
and a boundary condition as:

CVb = gbc-1 CV*(P) exp( AEf / kT) 3-14

which in conjunction with Eq. 3-12 precludes Frenkel pair generation at the
loop layer boundary, consistently with the situation at the core boundary of a
single dislocation loop [78].

3.1.4 Coalescence and Dissolution of the Dislocation Loops

The TEM experiments show that the density of dislocation loops created
by Si ion implantation generally decreases during oxidation, whereas the size
of them increases [81]. The model based on the triangular distribution
represents the loop distribution change in agreement with the observations.
The radii and the density of the loops ( or the unit distance L between the
loops as shown in Figure 3.4 ) can be correlated with the number of Si atoms
bound by the dislocation loops per unit area (N), considering that the density
of the majority size loops Dp is equal to 0.5 Lp-2 :

N = na a R7Di = 7 na R2fd(R) dR
i JRai

= m L2 [(Rp + Rmax)(Rp + Rmin) + Rax + Rmin] 3-15

where na is the atomic density of Si atoms on (111) plane ( = 1.5x1015 cm-2 ),
m is the ratio of Dall and Dp as defined in Eq. 3-6. Time-derivative of Eq. 3-15
relates the growth rates of the loop radii with the rate of their coarsening in
terms of dLp/dt. Furthermore, the net capture rate dN/dt of the dislocation
loop layer should be proportional to the rates of emission and absorption of
point defects at the layer boundaries:

dN= aN aN dL aN dRmin aN dRmax aN dm 1(1R)
dt -- = Rp '+ LpdRp + BRmin dRp + Rmax dRp + mdRp J t )

= aKIL( CI CIb) aKVL( CV CVb) Iat loop layer edges 3-16

where a is an effective cross-section of the loop layer in the unit of linear
length, KIL, KVL, CIb and Cvb are the same parameters of the loop reactions
with interstitials and vacancies as given in Eqs. 3-11, 3-12, and 3-14. It is
noticeable in Eq. 3-16 that the capture and emission rate constants of
dislocation loops are the product of a and KIL or KVL for interstitial and
vacancy, respectively. The partial derivatives of N, i.e., NN/ARp, NN/ALp,
aN/ARmin, aN/fRmax, and NN/am can be easily calculated from Eq. 3-15 as
closed-form expressions in terms of those distribution parameters. It is
assumed in Eq. 3-16 that the linear measures of the loop distribution (Lp,
Rmin, Rmax) have experimentally observable simple relationships with respect
to the majority loop radius Rp during the coalescence and dissolution
processes. For the loop growth, the relation between dRp/dt and dLp/dt, i.e.,
dLp/dRp can be extracted statistically from various TEM data as an analytic

function of Rp and Lp. It is determined by the more easily obtainable rate

dLaii/dRave, which is empirically extracted to be approximately 0.9 according
to the oxidation data [81]. In the current model, a constant value for

dRmax/dRp (= 2.3) is also extracted from the statistical interpretation of TEM
pictures [81], as shown previously in Figure 3.8. The dRmin/dRp is set to be
equal to zero, and dm/dRp is analytically calculable from the definition of m
in terms of Rmin and Rmax, as far as the loop distribution from TEM data is
generated with constant sampling period AR at each time step. Even if AR is
not maintained to be constant in time in the extraction of the loop radius
distribution from the TEM data, it is possible to calibrate the raw TEM data

before the asymmetric triangular distribution is applied. When the loops
shrink, the loop density is assumed to be constant, so dLp/dRp is equal to
Thus, Eq. 3-16 makes it possible to model the loop coalescence during
oxidation and annealing in a phenomenological way. Figure 3.11 shows the
scheme for modeling the evolution of the distributed dislocation loops,
which was implemented for simulation. The point defect distribution
around the layer of dislocation loops determines the capturing process of Si
atoms by the loops, and the majority loop growth rate aRp/at can be solved at
each time. The maximum radius and the density of the loops are obtained
eventually from their empirical relationship with Rp. Although this model

does not physically represent the actual loop-to-loop interactions or Ostwald
ripening, it encapsulates the statistical phenomenon through the change of
loop size and density restricted by the empirical triangular distribution func-
tion. The loop coalescence during oxidation is efficiently modeled through
the ratio of interloop distance variation with respect to loop radius as ob-
served in the statistical interpretation of TEM data. In the following sections,

Rmin, Rmax, Rp, Dp, AR Pressure from a single loop

Density distribution function D(R) Average pressure from one-
isized loops

Average radius
Total effective pressure P
Total density of the loops
from all the loops
p ~ ~----------------------"
Number of captured silicon atoms

Other Perturbation
Point defect concentrations around the loop layer *"" e.g., Oxidation

Figure 3.11. The scheme used for modeling the evolution of dislocation loops
based on pressure and point defects.

the loop reaction constants KIL and KVL will be more accurately estimated by

simulating and quantifying the loop evolution and the reduced OED of boron

in the buried layer region deeper than the dislocation loop layer.

3.2 Simulations of the Loop Evolution during Oxidation

The evolution of dislocation loops was simulated in a two dimensional

process simulator FLOOPS in which the above model was implemented. The

point defect parameters used in the simulation are shown in Table 2.1. As

described in Chapter II, those parameters were extracted by a dopant-defect

pairing model for both experiments on the low dose Si implant damage-

enhanced diffusion and the oxidation enhanced diffusion of boron and

phosphorus. Since the dislocation loop evolution is determined by the point

defect behavior, it is crucial to use the point defect parameters consistently for
all the diffusion simulations.
The data of loop growth are from the TEM experiments by Meng et al.
[81] on the dry oxidation-induced evolution of the end-of-range dislocation
loops introduced by Si ion implantation at 50 keV into (100) Si wafers. The
two Si implant dose conditions ( 2x1015 and 5x1015 cm-2 ) determine the
initial dislocation loop geometry. Since the oxidation was preceded by a

preannealing step at 5500C for 16 hours in the experiment, the initial point
defect concentrations were reasonably assumed to be at equilibrium in the
simulation. The data provided the loop distribution characteristics at each
oxidation time in the form of histogram. There are two sets of loop
evolution data corresponding to the Si implantation dose conditions.
Figure 3.12 shows the variation of pressure around the dislocation loop
layer during a four-hour oxidation at 9000C as obtained from the simulation
for the 5x1015 cm-2 dose condition. The maximum pressure at the layer
center is about the same order of magnitude as typical average values beneath
nitride films. As the loops grow, the pressure gradient around the loop layer

decreases, while the maximum and the minimum pressures remain in the
same range. This observation from simulation suggests that the strain is not
relaxed on the average during the growth period of the dislocation loops.
This simulation result is based on the theoretical calculation of the pressure
from individual loops, and it agrees with intuition qualitatively. However, it
will be necessary to directly verify the pressure in its quantity and distribution
by measuring the magnitude of the strain introduced by the dislocation loops
and by monitoring how the strain distribution changes upon annealing.
Zaumseil et al. [88] showed that it is possible to measure the strain from the
dislocation loops by using triple crystal x-ray rocking curve analysis. The x-ray

5.0x109 - -Pat 2 hours
Si implant dose = 5x1015 cm-2
- -P at 4 hours

E 3.0x1 09-

| 2.0x109 -

1.0x10 -

9 -1

-2.0x1 09 "
-2 .0x109 I I I I I I. I. I I I I I I I-
0 0.1 0.2 0.3
dislocation loop layer center Depth (gpm)

Figure 3.12. Variation of the pressure around the dislocation loop layer
during dry oxidation at 9000C for 4 hours, simulated in FLOOPS.

diffraction studies will lead to evaluation of the strain field around the

dislocation loop layer, and the calculated average pressure in Figure 3.12 can

be confirmed by the measurements.

In Figure 3.13, the simulation shows variation of the free unpaired

interstitial distribution around the dislocation loop layer during the dry

oxidation at 900'C for 4 hours where interstitials are injected from the

surface. The compressive pressure inside the loop layer as shown in Figure

3.12 causes the local dip in the interstitial distribution inside the loop layer

region. The model correctly represents the interstitial movement towards the

loop layer caused by the local variation of the pressure and the boundary

conditions restricted by the loop self-force. Growth of the loops can be

visualized from the change in the interstitial distribution near the layer

? 5x1013



_ 5x1012


- -0.0 hr. with loops

--0.5 hr. with loops

- -1.0 hrs. with loops

----- 2.0 hrs. with loops

- -- 4.0 hrs. with loops

--- --4.0 hrs. without loops

Si implant dose = 5x1015 cm-2

'- - -2 7


- S
- - - - - -

vA,' 'I I .
0 0.1 0.2 0.3 0.4 0.5
dislocation loop layer center Depth (pm)

Figure 3.13. Simulation of the variation in free interstitial distribution
around the dislocation loop layer with its center at 0.15 pm during the dry
oxidation, compared with the case without the dislocation loops.

boundaries from 0 to 4 hours. The interstitial supersaturation is already

limited at short times and gradually diminishes as time goes on, since more

interstitials are captured by the growing dislocation loops. The reduction of

interstitials due to the capturing action of the loops is well compared with the

case without the loops, which is the normal OED simulation as described in

the previous chapter. As seen in this log-scaled plot, the dislocation loops

capture most of the injected interstitials, which suggests the substantially

strong interaction of the dislocation loops and the interstitials. Figure 3.13 is

the result from the FLOOPS simulation with one parameter set that

simultaneously generated Figure 3.12 and other pictures in this chapter,

including the next section.


900C Oxidation
5 0 .0 , , ,1 , ,,,, ..L
- .- .--without dislocation loops
---- with dislocation loops
C ( Si implant dose = 5x1015 cm-2)



0 1 .0 . . .. .. .

0.01 0.1 1 2 4 10
Oxidation time (hour)
Figure 3.14. Simulation of interstitial supersaturation CI/CI* at 0.5 umn depth
position during the four-hour oxidation, with and without the dislocation

In Figure 3.14, the model predicts that the interstitial supersaturation at

0.5 pm depth position is significantly lower in the presence of the dislocation

loops than in the case without any extended defects during most time period

of oxidation. The simulation result shown in Figures 3.13 and 3.14 agrees

with several experiments on transient dopant diffusion due to high dose

implant damage, which suggest that the dislocation loops in silicon work as

an efficient sink for interstitials [19, 71, 74, 89]. The degree of reduction in the

CI supersaturation mostly depends on the loop-interstitial reaction constant

KIL. The value for KIL is determined by fitting the simulation with the TEM

data and SIMS data of boron OED that will be discussed in the next section.

The reaction constant KVL between loops and vacancies also affects the loop

growth during oxidation, but its value is much smaller than KIL at 9000C.

9000C Oxidation
1.8x1015 . . ..... ,
IE Data; 2x1015 cm-
o 1.6x1015
S x --8--Simulation; 2xlo15cm-2

0 1.4x1015- Data; 5x1015 cm-2
S- O Simulation; 5x1015 cm-

5 1.0x1015

S6.0x1014- --

4.0x1014 '
0 1 2 3 4 5
Oxidation time (hour)
Figure 3.15. Density of the Si atoms bound by dislocation loops as a function
of oxidation time from data (ref. [81]) and FLOOPS simulation in the two
different cases of Si implant dose, 2x1015 and 5x1015 cm-2.

More accurate evaluation of KVL may be possible through more experiments

in vacancy injection environment such as nitridation or silicidation.

Figure 3.15 shows a good agreement between the simulation and the

data [81] on the temporal change in the number N of Si atoms bound by the

dislocation loops per unit area during the oxidation at 9000C. There can be

two ways to obtain the data for N from the plan-view TEM pictures. First,

direct pattern recognition of the shaded area in the PTEM picture from the top

will lead to a good measure of the number of atoms confined in the non-

circular loops as shown in the picture. However, this method should

consider the random orientation of the loops in <110> or <111> directions,

and a proper projection of the loop area will be required. Second, N can be

calculated from the measured radii and density of the loops, assuming the

loops are all circular. This was the method adopted for the extraction of the
data [81] fit in this work. This method will involve a certain range of error by
the assumption of circular loops, and the error range increases at larger times
when the loops are more non-circular due to the network formation. In
Figure 3.15, the constant 5% error bar at 4 hour condition may be under-
estimating the probably larger error in the TEM measurement. Considering
the overall fit, the simulation result implies that the two-dimensional
simplification can lead to an effective model for the interaction of extended
defects and point defects, which is an essentially three-dimensional
Figure 3.16 shows the variation in total density of dislocation loops
during oxidation from the simulation and the data [81]. Again, non-
uniformity of loop size and shape is seen in the TEM pictures typically at
larger times; when the loops are forming a network, there can be a sizable
difference in the density of the actual non-circular loops of various sizes and
that of the loops modeled to be circles. In this case, the meaningful loop
radius and density depend more on statistical interpretation of the TEM
pictures. The missing data point at four hours for 5x1015 cm-2 dose condition
in Figure 3.16 corresponds to this case. However, the simulation at large
times in general has been improved by modeling the loop coalescence using
the triangular distribution function for loop radius. The preliminary model
with equal-sized loops led to a larger discrepancy at the four hour condition
The average radius changes during the oxidation as shown in Figure
3.17. The data show little difference in the initial Rave between the two
silicon implantation dose conditions, and the simulation predicts that the
loop size will vary at almost the same rate during oxidation, consistently with


0 4x1010


Z 2x1010

C x

, , I , I ,

I , I ,

, I I 1 I I 1 I I I I 1 1 1 1 I I I
0 1 2 3 4 5
Oxidation time (hour)

Figure 3.16. Total density of the dislocation loops from data
simulation in the two different Si implant dose conditions.

9000C Oxidation
. . I . . I I


(ref. [81]) and

2 3 4 5
Oxidation time (hour)

Figure 3.17. Variation in the average radius of the dislocation loops from data
(ref. [81]) and simulation in the two Si implantation dose conditions.

Data; 2x1015 cm-2

-E Simulation 2x1015 cm-2

Data; 5x1015 cm-2

- O Simulation 5x1015 cm2

Data; 2x1015 cm-

---Simulation; 2x1015 cm2

Data; 5x1015 cm-2

- Simulation; 5x1015 cm"2



500 -

400 -




: : : : , ,


I ,

the data. Figure 3.18 shows that the variation of the loop distribution

parameters in the simulation agrees roughly with the data [81]. The error

range of Lp and Rmax especially at large times depends on interpretation and

calibration of the actual non-triangular distribution of the loops. Modeling

the loop coalescence with the distribution function has allowed us to

simulate the maximum radius of the dislocation loops, which is an

important factor to consider in designing short channel MOSFETs and

shallow junctions formed by high dose ion implantation. The model can be

further developed to predict the tolerable implantation and annealing

conditions that prevents the dislocation loops of maximum size from

reaching the p-n junction region and from causing detrimental leakage

currents. More systematic data extraction from TEM pictures will lead to

9000C Oxidation

,3000 -
-P- :
S2500 -

- 2000-
: -
2 1500



Lp data

--- Lp simulation
Rmax data

- O Rmax simulation

-0 Rp simulation
- A-- Rmin simulation

Si implant dose = 2x1015 cm-2

--A A --- -----A -- --- -- ---- -A
^-'- -0-----0~---
_-A.- _-- _-- -- ._-. .......- _A

0 1 2 3 4 5
Oxidation time (hour)
Figure 3.18. Variation in the radii of the dislocation loops and the unit
interloop distance Lp for the loops of radius Rp, which are obtained from data
(ref. [81]) and predicted by FLOOPS simulation.

more accurate estimation of the rates of growth and coarsening of the
dislocation loops.

3.3 The Effect of Dislocation Loops on OED of Boron

It is necessary to extract accurate values of the point defect capture rates
of dislocation loops in order to predict the distribution of point defects around
the loop layer more quantitatively. In this work, it was achieved by
monitoring the reduction in oxidation enhanced diffusion (OED) of boron in
the region deeper than the dislocation loop layer. This modeling work uses
the SIMS data from an experiment by Meng et al. [90] where boron buried
layer was used to probe the reaction between point defects and dislocation
loops during oxidation. In that experiment, a thin (< 500 A) boron marker
layer was grown epitaxially on <100> Czochralski silicon wafers, followed by
epitaxial growth of an overlayer of 6000 A of undoped silicon. As-grown
boron profiles were obtained by SIMS. The end-of-range dislocation loops
were introduced into some of the samples by Si implantation at 50 keV to a
dose of 2x1015 cm-2. All the samples were then annealed at 9000C from 10
minutes to 4 hours in either nitrogen or dry oxygen ambient. Thus the
samples were categorized into four different groups: (1) without loops and
annealed in N2; (2) without loops and annealed in dry 02; (3) with loops and
annealed in N2; (4) with loops and annealed in dry 02. For the case (4), the
samples were preannealed at 9000C for 10 minutes in nitrogen ambient prior
to oxidation so as to form the dislocation loops and anneal out excess point
defects [90]. The depth locations of the dislocation loop layer and of the boron
buried layer are at about 0.15 pm and 0.6 pm from the surface, respectively.
Boron redistribution during oxidation at 9000C was simulated with the
model described in the previous sections. The initial loop distribution

observed in the TEM experiments [81] were used in the simulation, since the

implant conditions are the same for the two experiments. The dislocation

loop parameters used in the simulation are shown in Table 3.1. Those are the

same parameters used for simulation of loop evolution in the previous

section. To complete this table, it will be necessary to estimate the

temperature dependence of some parameters such as KIL and KVL by

monitoring the loop evolution at temperatures other than 9000C.

Table 3.1. Parameters and constants of the interaction between dislocation
loops and point defects used in the simulation.


KIL (loop-interstitial reaction constant) @9000C

KVL (loop-vacancy reaction constant) @9000C

a (capture and emission cross-section)

ro (radius of a silicon interstitial)

rs (radius of a vacancy)

e dilatationn factor of a silicon interstitial)

dRmax/dRp @9000C

dRmin/dRp @9000C

dLail/dRave @9000C


4.3x103 sec-1

1.0x10-4 sec-1

7.5x10-7 cm

1.11 A

2.47 A






In Figure 3.19, the simulation correctly shows the reduction in OED due
to the dislocation loops capturing the interstitials injected from the surface

during 2 hour annealing in dry oxygen. For the case with the loops, the

simulation started from the Gaussian fit to the B profile obtained after the

preannealing step required for loop formation. Figure 3.20 also shows a very




0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1
Depth (glm)
Figure 3.19. SIMS data [90] and FLOOPS simulation of boron diffusion with
and without dislocation loops during the dry oxidation at 9000C for 2 hours.

1 019

1 018

1017- b

1016 -
0.2 0.3 0.4 0.5 0.6 0.7
Depth (pLm)

0.8 0.9 1

Figure 3.20. SIMS data [90] and FLOOPS simulation of boron diffusion with
and without dislocation loops during the dry oxidation at 9000C for 4 hours.


---- Data: No lop; oxidation - Simulation: No loop; oxidation

Data: With loops; oxidation - Simulation: With loops; oxidation

S -A- Data: No loop; in nitrogen -A Simulation: No loop; in nitrogen

= At 9000C
o 11
S 1

I 7

"n 5
!E 0---. j^ ^_ ,-- -

o 3
m -i A-

0 1 2 3 4 5
Time (hours)
Figure 3.21. The reduction in OED of boron due to the dislocation loops
shown in terms of diffusivity enhancement, also compared with the case of
nitrogen ambient without the loops.

good match between simulation and measurement of boron profiles after 4

hour annealing in dry oxygen.

The enhancement in boron diffusivity was quantified through profile

matching with a boron diffusivity DB* under inert intrinsic condition [24].

Figure 3.21 shows the extracted diffusivity enhancement with respect to the

DB* for the simulation and the SIMS profiles for the three different groups of

samples (1), (2), and (4). The data show approximately 50 % reduction in OED

of boron when the dislocation loops exist. The simulation agrees with the

data fairly well for all three time conditions. The seemingly larger error for

the half an hour condition is due to the larger error range in calculating

effective diffusivity at shorter times. Actual junction depth movement,

Data: No loop; oxidation

-p Simulation: No loop; oxidation 900C
9000C Oxidation
120 ---Data: With loops; oxidation .L
m - Simulation: With loops; oxidation
o- 100

g 80


0 40

2 0 . i.i .
0 1 2 3 4 5
Time (hours)
Figure 3.22. Re-presentation of OED data and simulation shown in Figure
3.20 in terms of the linear measure of junction depth movement ] Dff t
with 10 % error bars applied to all the time conditions.

which is proportional to square root of product of effective diffusivity and

anneal time, does not show such strong time dependence in error range.

Figure 3.22 re-presents the displacement of boron profiles during oxidation in

terms of Deff t rather than Deff/DB*. Since the unit is linear measure of

length, the constant 10 % error bars in the Figure 3.22 are now meaningful.

Comparison of Figure 3.21 with Figure 3.22 tells that the representation of

enhanced diffusion in terms of effective diffusivity can be misleading when it

is used for matching simulation and measurement. More detailed discussion

on error range interpretation has been provided previously by the author [48].

Figure 3.21 also shows that the diffusivity is enhanced by factor of 3 or 4

during annealing in nitrogen. The unusual enhancement of boron diffusion

in the nitrogen ambient suggests that thin native oxide or oxygen precipitates
at the epi/substrate interface have worked as an interstitial injection source in
the bulk. The SIMS profile of oxygen peak at the epi/substrate interface
shown in Meng et al. [90] provides an experimental basis for the possibility.
The possibility was considered in the simulations by implementing an
interfacial oxide layer injecting a reasonable amount of interstitials in the
bulk. For the simulation of the profile in the absence of dislocation loops, the
surface injection level was chosen to fit the measured profile. To obtain the
gaussian profiles, the location of interfacial oxide was assumed to be at 1.0 gm

from the surface. The purpose of this imaginary injection source is to
approximately mimic the effect of oxygen precipitates. The same bulk
injection level was used for the case with the dislocation loops, so that we can
estimate the effectiveness of interstitial sink action of dislocation loops.
Figure 3.23 shows the simulation and the data for boron profiles after 4
hour annealing in nitrogen ambient with and without dislocation loops. The
simulated profile movement is less than the measured one, and the

interstitial sink action of dislocation loop layer is slightly overestimated in
the simulation with the reaction rates in Table 3.1. However, it should be
noted that the same reaction rates and the bulk injection source were used for
all the best-fit simulations in the previous figures, which are more critical
than the nitrogen ambient case in Figure 3.23. In addition, the simulation at
other time conditions is within the error range of SIMS measurement as
shown in Figure 3.24. Since the boron buried layer is located between the
interstitial injection source at 1.0 gm and the dislocation loop layer near the
surface, the interstitial sink action of loops is not so pronounced as in the case
of oxidation, where the interstitials injected at the surface are "screened" by
the dislocation loops. The overall effects of dislocation loops shown in

University of Florida Home Page
© 2004 - 2010 University of Florida George A. Smathers Libraries.
All rights reserved.

Acceptable Use, Copyright, and Disclaimer Statement
Last updated October 10, 2010 - - mvs