POINTDEFECTBASED TWODIMENSIONAL MODELING OF
DISLOCATION LOOPS AND STRESS EFFECTS ON
DOPANT DIFFUSION IN SILICON
By
HEEMYONG PARK
A DISSERTATION PRESENTED TO THE GRADUATE SCHOOL
OF THE UNIVERSITY OF FLORIDA IN PARTIAL FULFILLMENT
OF THE REQUIREMENTS FOR THE DEGREE OF
DOCTOR OF PHILOSOPHY
UNIVERSITY OF FLORIDA
1993
Copyright 1993
by
Heemyong Park
ACKNOWLEDGMENTS
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, JungMee, for her love,
support, and everlasting inspiration.
I acknowledge the financial support from SEMATECH.
TABLE OF CONTENTS
ACKNOWLEDGMENTS..................................................................................... iii
A BSTR A C T ........................................................................................................... vii
CHAPTERS
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
II MODELING OF LOW DOSE SILICON IMPLANT DAMAGE
EFFECTS AND OXIDATION ENHANCED DIFFUSION....................... 19
2.1 DamageEnhanced and OxidationEnhanced Diffusion................... 19
2.2 A Model for Dopant Diffusion Based on Pairing Theory...................... 20
2.2.1 Diffusion Equations Accounting for the DopantDefect Pairs..... 21
2.2.2 Significance of Binding Energies in DamageEnhanced
D iffusion............................................................................................ 24
2.3 Defect Equations Incorporating Dual Reaction with Traps.............. 28
2.4 Simulation of DamageEnhanced 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
III MODELING OF THE EVOLUTION OF DISLOCATION LOOPS AND
THEIR EFFECTS ON OXIDATION ENHANCED DIFFUSION OF
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
IV MODELING OF THE STRESS EFFECTS ON DOPANT DIFFUSION
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 TwoDimensional Extension and the Effects on Threshold Voltage
of ShortChannel 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
POINTDEFECTBASED TWODIMENSIONAL MODELING OF
DISLOCATION LOOPS AND STRESS EFFECTS ON
DOPANT DIFFUSION IN SILICON
By
HEEMYONG PARK
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 dopantdefect pairing is found to be crucial in modeling
the implantation damage effects, and the effective binding energies for boron
defect and phosphorusdefect pairs are experimentally determined. The
extracted parameters provide an important reference for further modeling of
diffusion under high dose implantation conditions involving extended
defects.
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 looptoloop 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 pointdefectbased atomistic model for the stress effects on dopant
diffusion is developed by accounting for variation in formation enthalpy of
dopantdefect pairs due to the hydrostatic pressure. Binding energies and
diffusivities of dopantdefect 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 oxidepadded nitride film of various widths is caused by
the stress at the film edge. Twodimensional simulation of diffusion in the
pressure field leads to better prediction of threshold voltage shift in short
channel MOS transistors.
viii
CHAPTER I
INTRODUCTION
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
interest.
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 physicsbased 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, physicsbased modeling is strongly demanded. Doping
profile and junction depths in multidimensions 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 nonequilibrium
environments initiated by the other processing steps, typically by ion
implantation.
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
Highdose ion implantation is crucial in obtaining highlydoped regions
in silicon such as the source and drain of MOSFETs and DRAM cells. The
highdose 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 crosssection 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 selfenergy in
the crystalline structure. By intuition, we can visualize the stress distribution
I I I
II I
I Loop Radius R
I I
3D view
Figure 1.1. A crosssection 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 nonequilibrium
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 scaleddown 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 spacecharge 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 physicsbased 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 partiallyfilled band of edgestates analogous to twodimensional
surface state bands. The broken bonds in the extra halfplane 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 ntype Ge with
plastic deformation [4, 5]. He assumed that all conduction electrons are
expelled from the vicinity of a negatively charged dislocation in ntype
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 halffilled dislocation band in the neutral state,
which accounts for the donorlike behavior of dislocations in ptype 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
processing.
Recently, Ross et al. [9] measured the reverse leakage current of GeSi/Si
pn 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 generationrecombination 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 selfinterstitials [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 planview 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 postamorphization 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 cm2) 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 amorphouscrystalline 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 cm2,
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
1017
U Critical dose for Category I defects
"i amorphization
E 1o16
0o .. Category H defects
10
0
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 massdependent 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 cm2 for B, P, and Si; 5x1013
cm2 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 nonequilibrium diffusion mechanism of dopant atoms and
point defects in the absence of their complicated interaction with extended
defects.
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 emitterbase 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 emittercollector leakage. When they extend through the
narrow base, they become emittercollector pipes which allow significant
emittertocollector 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 stressinduced 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 4Mbit and 16Mbit
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 LOCOSbased 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 stressinduced defects. At the corner
locations of the trenches, larger stresses can be generated. Stress buildup
becomes worse as the oxidation temperature is decreased.
Figure 1.3 shows a simplified crosssection of a trench capacitor DRAM
cell with the substrateplatetrench structure [17]. Without proper control
over the processes, the leakage current may be a very significant problem in
Bit line
/ Word line
Dislocation loop
psubstrate
Th
"p
po+ poly
in oxide
SiO2
Ii0
SiO2
pepi

Figure 1.3. A simplified crosssection of a trench capacitor DRAM cell; ref.
[17].
Bssssssssssssssa
this structure due to superimposed effects of the implantationinduced
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 highdose ion implantation as well as the stress
induced edge dislocations during the oxidation processes. Figure 1.4 shows a
simplified crosssection of a typical nchannel lightlydoped drain (LDD) MOS
structure, which features the selfaligned phosphorus doped n regions
between the channel and the n+ source and drain regions doped by highdose
arsenic implantation. During the annealing cycle for activating the
source/drain region, a layer of dislocation loops is produced near the
source/draintosubstrate junctions. The dislocation loops, which can cross
the junctions, may trap the fastdiffusing impurities and provide generation
recombination centers at their sites. This can lead to large leakage currents
through the reversebiased source/drainbody 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
Sptype substrate
LDD region
Figure 1.4. A simplified, unsealed crosssection view of an nchannel 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 arsenicimplanted region occurs during the subsequent anneal. The
local redistribution of dopant in the substrate can make a significant
difference in the VT rolloff in the short channel MOSFETs [18]. The observed
boron pileup was exactly centered around the arsenic implantinduced
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 postamorphized 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 fastdiffusing 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. Oxidationenhanced diffusion (OED) is the most
widely measured and understood nonequilibrium 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 shortchannel effect [21, 22]. It was observed that the
threshold voltage initially increases as the channel length decreases until the
final VT falloff begins. The anomalous VT behavior was attributed to the
OED effects on redistribution of the channel dopant (boron in nchannel
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 shortchannel
MOSFETs, as mentioned in the previous section.
It has been reported [23, 24, 25] that low dose silicon implantation also
provides nonequilibrium diffusion conditions by creating point defects or
defect busters. If the silicon implant dose is below 2x1014 cm2, 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 shortchannel devices. The lightlydoped drain regions shown
in Figure 1.4 are made typically through phosphorus implantation with a
dose range of 1013 cm2, 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 lowdose 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, physicsbased models are developed for the
interactions among dopant atoms, point defects, and extended defects.
Materialrelated 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
thesis.
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 SUPREMIV
is reexamined 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 physicallymeaningful quantity that governs the transient
diffusion. A set of point defect parameters including effective binding
energies of dopantdefect 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 implantationinduced 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 ObjectOriented 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 pressuredependent dopant diffusion equation is
derived by accounting for the variation in enthalpy, binding energy, and
diffusivity of dopantdefect 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.
CHAPTER II
MODELING OF LOW DOSE SILICON IMPLANT DAMAGE EFFECTS AND
OXIDATION ENHANCED DIFFUSION
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 dopantdefect pairing is found to be crucial in
modeling the implantation damage effects, and the effective binding energies
for borondefect and phosphorusdefect pairs are estimated from the
simulations.
2.1 DamageEnhanced and OxidationEnhanced 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 implantationinduced
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 nonequilibrium 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 nonequilibrium, 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 dopantdefect 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 dopantdefect reaction such as Frank
Turnbull mechanism recombinationn of vacancy with dopantinterstitial pair)
and dissociation mechanism recombinationn of interstitial with dopant
vacancy pair). Those mechanisms have been used in explaining diffusion of
fastdiffusing 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 FrankTurnbull
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 DopantDefect Pairs
When it is assumed that local equilibrium is attained between the
dopant, the defects, and the dopantdefect pairs at each point in the diffusing
dopant profile, the following equations can be applied to modeling of
damageinduced transient diffusion. The dopant equations proposed by
Mathiot and Pfister [33] are modified to limit the maximum number of
dopantdefect pairs as described in Law et al. [45]:
JAX dAcKAxGxCcxc))Cl C Vlog CA+ C 21
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)1c, 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
noninert condition in terms of C~x as follows [46]:
Cx1 = Gxc x 22
Cx
where Cx = GxC'xc, and Cx = C Cxe 23, 24
C C
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 dopantdefect
pairs can be expressed by a simple mass action relationship with a pairing
coefficient KAXC:
CAX = KAX CX1 CA 25
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 26
C C
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. 25 is
expressed as [20, 32, 47]:
0 EbAx(\
AXc bAXc
KAXC = Cs exp WkTJ 27
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 dopantvacancy and dopantinterstitial pairs in
silicon. In fact, Eq. 27 is based on diffusion theory via a vacancy mechanism
without considering the charge states of the defects. However, it has been
extended selfconsistently 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 CX1
CAX= KAX CX' C1 KAXCX CA= Cx CIA 28
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 dopantdefect pairs AX in the
same way as Eq. 27:
OAAX (EbAX 29
KAX = S exp kT 29
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 lightlydoped 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
pair.
2.2.2 Significance of Binding Energies in DamageEnhanced Diffusion
Since ion implantation creates such a large number of excess defects, it
becomes critical to account for the dopantdefect 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. 21, 26, and
28 can be used to derive an equation for the instantaneous diffusivity
enhancement of a dopant A under intrinsic conditions:
DA CAI CAV
fAI + ( 1 AI )
A I C + +n t IKAIC +KAVC] 210
/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. 210
becomes
DA KAICI +(CI/C 211
DI KAICI+ 1
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
OED.
(b) When CI / Cl exceeds the threshold (KAI C ) 1, i. e., when
KAI CI > 1, Eq. 211 is approximately
DA
S= 1+ (KAIC)1 212
Eq. 212 characterizes the damageenhanced diffusion during the initial
transient period. Under intrinsic doping conditions, Eq. 212 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. 211 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 dopantdefect
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=5x1017 c3
 KAI=5x1016 cm3
10  _ KAI=5x1015 cm3
 with CAI neglected
0< 104
C/ I* 
102
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
C,/C,*
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 dopantdefect pairs.
considering the pairing model actually showed excessively enhanced
diffusion of phosphorus and boron with implant damage.
We can draw from Eq. 212 a very important fact about the equilibrium
concentration of interstitials and the binding energy of dopantinterstitial
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 timeaveraged 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. 212.
From Eq. 212, 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. 29. 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  213
where Cs is the number of available lattice sites in the crystal, HI is the
formation enthalpy of silicon selfinterstitial, 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. 213 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. 29 and Eq. 213, we obtain the relationship
of the factor KAI C}I with temperature and energies as follows:
KAI CI = OAI 1 exp ) exp (. .EbAI k 214
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 dopantinterstitial 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 nonequilibrium 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 dopantvacancy binding energy EbAv should be smaller than the
vacancy formation enthalpy Hy. These relationships provide a crucial
restriction on physicsbased 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 215
The above reaction accounts only for the reverse reaction as the dissociation
of the filled trap IT (i.e., interstitialtrap 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
trapmediated recombination of point defects:
V + IT T 216
In Eq. 216, 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
terms:
aCET
S= RTI + RTV 217
at
where RTI and RTV are the rate of empty trap concentration change due to the
interstitial reaction in Eq. 215 and the vacancy reaction in Eq. 216,
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)] 218
RTI =tCT ET
RTV = KtrapV (CT CET) CV CT ET C* CET] 219
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. 218 is similar to that given in Law [28],
but the trap equation (Eq. 217) 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 220
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. 220,
the trap reaction rates can be formulated as follows:
47C aExI _AEET I 221
KtrapI = CT1 D exp AEET 221
4xc asrv AEFTV
Ktrap =4 a DV exp AkT) 222
Kap TCs
where AEETI and AEFrV are energy barriers of recombination of empty trap
interstitial and filled trapvacancy, respectively. The corresponding capture
radii aETI and aFTV can be assumed to be equal to aiv for direct IV
recombination.
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 223
3t
(Cv + CA) = V(Jv + JAV) KR (CI Cv CCv) RTV 224
3t
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. 218. In
Eq. 224, Jv, JAV, and RTV are similarly defined. Assumption of detailed
balance shown in Eqs. 218 and 219 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(CIC) = g 225
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 reinjected into the crystal as interstitials. The
surface recombination velocity, KI, can also depend on the surface growth
rate. The surface recombination is expressed as:
226
KI = KImax 0ox +1 Kmin
Vox
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 oxidationenhanced 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
lowtemperature CVD (LTCVD). As a tentative means of testing the trap
mediated diffusion model in this section, two of their results are simulated
with SUPREMIV 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 asgrown profile is a rough approximation to the SIMS data shown in
[51]. The total trap concentration CT used in the simulation is 1x1017 cm3, 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 cm3. Both Figure 2.2
1019
1018
1017
1016
1015
1014
1013
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. SUPREMIV simulation of the boron diffusion at 9000C with total
trap concentration equal to 1x1017 cm2, emulating the data [51] for the case of
MBEgrown epi silicon.
 Boron at 0 min.  Interstitial at 0 min.
1 019 Boron at 20 min.  Interstitial at 20 min.
1018
S10L A 
0 15_ L L J 
a 1017
0
o 1014
1013 _
1012 I
11
0.0 0.2 0.4 0.6 0.8 1.0 1.2
Depth (glm)
Figure 2.3. SUPREMIV simulation of the boron diffusion at 9000C with total
trap concentration equal to 1x1019 cm2, emulating the data [51] for the case of
FGCVDgrown 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
depthdependent boron diffusion for the two cases. In the MBEgrown
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 FGCVDgrown epilayers. Although the assumed value
of CT (1x1019 cm3) 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]. Nonequilibrium
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 DamageEnhanced Diffusion and OED
The above model of dopantdefect pairing was implemented into
SUPREMIV, 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 damageenhanced 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 Preexponential Activation
DI
a
C*l
Dv
Cv
KImin
KImax
Kvmin
Kvmax
0
aXI
(Xv
KR
CT (Floatzone)
CT (Czochralski)
CET / CT
Ktrapl
KtrapV
DIB
fBI
EbBI
EbBV
DAI (w/ neutral I)
fPI
Ebpi
Ebpy
600.0 cm2/sec
5.0 x 1022 cm3
13.0 cm2/sec
2.31 x 1021 cm3
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
0.84
1.79 x 107
8.16 x 104 cm3/sec
1.0 x 1017 cm3
1.0 x 1018 cm3
2.39 x 105
8.16 x 104 cm3/sec
1.77 x 105 cm3/sec
1.7 cm2/sec
0.8
1.52 eV
0.9 eV
3.85 cm2/sec
1.45
1.49 eV
0.6 eV
2.44
2.36
2.92
1.08
1.67
2.73
2.50
5.39
1.82
0.05
 1.91
3.19
1.57 eV
2.94 eV
3.92 eV
3.56 eV
3.66 eV
0.05 eV
__________________________ ________________________________ I __________________________
All the damageenhanced 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 SUPREMIV, which numerically solves the dopant
and defect diffusion equations (Eqs. 21, 223, and 224) in the subsequent
diffusion steps. The program used for the initial defect calculation in this
work is a Monte Carlo model [52] implemented in SUPREMIII [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 asimplanted 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 selfinterstitials 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
energy.
 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
0
 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
1012
0 0.2 0.4 0.6 0.8
Depth (micron)
Figure 2.4. Asimplanted initial distribution of interstitials after the 28Si
implantation. Monte Carlo simulation from SUPREMIII 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 damageenhanced diffusion predicted by any point defectbased
model is contingent on the simulation of initial asimplanted 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 lowdose 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 boroninterstitial 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
(keV)
0 50 100 150 200 250
100 I, I, 1 .
75 
E A
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 cm2 as a function of implant energy and anneal temperature. The
data [24] and the SUPREMIV 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 asimplanted 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
dopantinterstitial pairs. DI and C* are the most important parameters to be
decided along with the binding energy. In this boron experiment, floatzone
substrate was used [23, 24], and the total trap concentration in that material is
assumed to be so small (= 1.0x1017 cm3) 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]. SUPREMIV 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. 212, 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 timeaveraged 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. 214
demonstrated. The EbBI for boroninterstitial 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 damageinduced 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. 225 and 226 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. 212. 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
42
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
50
2 5 . .. I . I I .
1011 1012 1013 1014 1015
Dose (cm2)
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 SUPREMIV simulations are compared. The anneal times are 30
minutes at 8000C, 10 minutes at 9000C, and 4 minutes at 10000C.
(Eqs. 223 and 224) under the assumption that the bulk recombination is the
only dominant process during the short time [48].
The simulation for 1x1012 cm2 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
43
Boron, 8000C, 29Si implant at 200 keV
1 7 5 , . .. 1 . . 1 'I ,
0 5x1013 cm2 Data  5x1013 cm2 Simu.
150 1x1014 cm2 Data I x1014cm2 Simu.
S ~ 2x1014cm2 Data  2x1014cm2 Simu.
125
C100  
75
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
SUPREMIV simulations are compared.
damage amount, because the inert diffusion tends to dominate the junction
motion.
The time dependence of the damageenhanced 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
:n8^^
S75100
50
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 cm2 with different energy conditions.
The data [24] and the SUPREMIV 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 SUPREMIV 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 IV recombination [20, 57,
58]. Moreover, dopantmediated 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 boronvacancy 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 nonGaussian
with values of EbBv larger than 1.0 eV. According to the physicsbased
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
dimensions.
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 A9000C Simu.
E
S. [] 1000C Data 1 1000C Simu. /
o 60 /
1100C Data 1100C Simu.
40
20 .
B
0.1 1 10 100
Time (min)
Figure 2.9. Phosphorus profile motion /Defft due to 28Si implants of dose
1x1014 cm2 at 60 keV as a function of anneal time and temperature. The data
[25] and the SUPREMIV 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 damageenhanced 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 floatzone material and Park and Law [25]
used Czochralski. The different number of traps contributes to different time
behavior.
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 dopantdefect 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 phosphorusinterstitial 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 phosphorusvacancy
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)
1019
1018
1017
1016
1015
1014
1013
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 nonGaussian dopant diffusion shown in
SUPREMIV 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
phosphorusvacancy 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 cm3 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 cm3 reduces the junction depth at 15 and 30 seconds, while at 60
minutes the results are the same as when CT is 1.0x1018 cm3.
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 diffusionlimited with an energy
barrier AEFTV 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
AEFTv. 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., trapmediated 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 interstitialtrap 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 EFTV = 1.9eV
  EFTV = 1.0 eV
4
104  Without traps
103
0102
101
0.1
103 102 101 100 101 102
Time (min)
Figure 2.11. Interstitial supersaturation at 0.1 im depth position for different
AEFTV values in the simulation.
interstitial supersaturation C / C* is simulated at one depth position.
Compared to the case without traps, the trapmediated 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 AEFrv 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, temperaturedependent 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 phosphorusinterstitial pair DP*I was kept as the default from Fair [65], and
vacancy contribution Dp* was added to match the shorttime 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 boroninterstitial binding energy
within its error range (about 0.02 eV). Therefore, all the extracted
parameters are selfconsistent in modeling both phosphorus and boron
diffusion.
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
interstitialvacancy equilibrium, which affects the time dependence of the
nearsurface 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
Oxidation
I I . .i
11000C 0
0* 10000C
Simulation
Simulation
* Data
*0
90
* ^^0
S I TI 1 I
103
Time in Minutes
00C
Sr 1 04
104
Figure 2.12. Lateral decay length of interstitial diffusion as a function of time
and temperature from the measured data [49] and the SUPREMIV
simulations.
2D
20
15
10
5
102
i1 r
I I I I .
I I
Interstitial Supersaturation @1100C
5 1 Packan
20A Initial oxide M Ahn
4.5 
Griffin
4 o Ahn 201pm membrane
A Ahn 55pm membrane
3.5  Simulation
S400A initial oxide
2.5
2
1.5 backside enhancement/ 0
O A
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
[54,67,68].
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
model.
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 interstitialvacancy 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 dopantdefect pairing.
Vacancy Undersaturation @1100C Oxidation
0 .5 I I.' ' . .. ' 1 . 11
 Simulation + Mizuo
SA Packan A Antoniadis
1
1.5 A
>
2 A
0* m
2.5
3
3.5
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 SUPREMIV 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 SUPREMIV 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 thermallydetermined quantity
KAI CI is found to quantify the enhanced diffusion of intrinsicallydoped
boron and phosphorus due to implant damage. The intrinsic doping
56
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.
CHAPTER III
MODELING OF THE EVOLUTION OF DISLOCATION LOOPS AND
THEIR EFFECTS ON OXIDATION ENHANCED DIFFUSION OF BORON
IN SILICON
High dose ion implantation is an essential technique for obtaining
heavilydoped 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
amorphouscrystalline 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 endofrange 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 crosssection TEM and planview
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)
Z*
ro (1+e)
Ab
h
x'
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 physicsbased 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
1.0x1010
"E 5.0x109
2. o
0
)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. 31 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() +K31
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 crosssection viewed in a twodimensional 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 crosssection 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
crosssection of a dislocation loop
L 
L L
t L
, /
R i
( t
[A]
,  
(0,0): x [D]
  L
I I
I I
I I
I I
I I
r
t
[
',[B]
I L
0
IIr
i I
I I
I I
I I
I I
,I..
6
SPloop(i) dx dy
i=l
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 equalsized 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 crosssection TEM pictures, where most of the category II end
ofrange 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
 
64
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=LR,y=0) [D]
0 0* /
a_
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 equalsized 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. 31 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
halfside 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 crosssection 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
magnitude.
Analytic functions for the average pressure are extracted in this work to
replace the timeconsuming numerical integration of the pressure
components. The average pressure from onesized 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
2.0x10
  R = 200 A; L = 1000A
E
> 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 bestfit 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
min
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
1.2x1010
9.0x109
6.0x109
3.0x109
0
3.0x109
" n 4 n

max ; R=100A
S
max; R=200A
  <>max R=400A
S
max; R=800A
: 
min ; R=10OOA
S
min ; R200A
 
min ; R400A
S
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.
1

. I t t
3.1.2 Asymmetric Triangular Density Distribution of Loop Radius
A preliminary model for the equalsized loops was reported previously
by Park and Law [82]. However, since TEM data show nonuniform 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 ac2; 9000C Oxidation
1x1010 I II I
o 0.5 hour
E
0 o 1 hour
8 8x109 a hours
a
SA 4 hours
S6x109  triangular distribution
O
.2 4x109 
SOO A
0
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 32
( Rmax Rmin ) ( Rp Rmin )
d(R) = 2 Dall ( Rmax R)33
fd(R (Rmax Rmin ) ( RmaxRp wen Rp < < max 33
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:
Rmax
Dall = fd(R) dR 34
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:
AR
rR+
D(R) = 2 fd(R') dR' = fd(R) AR 35
fR
2
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
Dp
o AR
oD2*
DiD
D1
Rmin Rp Rmax Radius
,4 time = tl
Q.
0
0
0
CO
"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. 32, 33, and
35:
2 Dall AR 2 Dali
Dp = fd(Rp) AR = Rmax Rmin m 36
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 fRax
Rave = C RiDi / Di = R/d(R) dR fd(R) dR
i i Rmtn JRmin
1
= 3 (Rmin + Rp + Rmax ) 37
Therefore, the quantities such as Rave and Dall, which are measurable from
the TEM pictures by pattern recognition without discretization, can be given
by closedform 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:
Rmax
P=
=
dR
i AR J
Rmax
= 2 DRal f
dR 38
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 onesized 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
4.0x109
3.0x109
$ 2.0x109
t
>1.
2 1.0x109
2 0
a_
1.0x1
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
averagesized loops
, where Dp = 0.5/Lp2 and Dall =
0.5/Lall2.
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 39
(PAVv\
Cv*(p) = Cv*(P=O) exp (I ) 310
where P is the total effective pressure given by Eq. 38. The above expressions
for the pressuredependent 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 pressuredependent equilibrium interstitial
concentration as shown in Eq. 39:
S= V [DI CI(p) V )] KR ( CI CV CI*(p) CV*(P))
KIL ( CI CIb) at loop layer boundaries 311
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) 312
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 selfforce of a dislocation loop during
emission and absorption processes at its edge [77]:
AEf= gb n 8Rave 2V 1 313
47x(lv)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. 37. Vacancy continuity equation can also be
derived in an expression similar to Eq. 311, with a different reaction rate KVL
and a boundary condition as:
CVb = gbc1 CV*(P) exp( AEf / kT) 314
which in conjunction with Eq. 312 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 Lp2 :
Rmax
N = na a R7Di = 7 na R2fd(R) dR
i JRai
= m L2 [(Rp + Rmax)(Rp + Rmin) + Rax + Rmin] 315
where na is the atomic density of Si atoms on (111) plane ( = 1.5x1015 cm2 ),
m is the ratio of Dall and Dp as defined in Eq. 36. Timederivative of Eq. 315
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 316
where a is an effective crosssection 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. 311, 312, and 314. It is
noticeable in Eq. 316 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. 315 as
closedform expressions in terms of those distribution parameters. It is
assumed in Eq. 316 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
zero.
Thus, Eq. 316 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 looptoloop 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 dopantdefect
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 oxidationinduced evolution of the endofrange dislocation
loops introduced by Si ion implantation at 50 keV into (100) Si wafers. The
two Si implant dose conditions ( 2x1015 and 5x1015 cm2 ) 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 fourhour oxidation at 9000C as obtained from the simulation
for the 5x1015 cm2 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 xray rocking curve analysis. The xray
5.0x109  Pat 2 hours
Si implant dose = 5x1015 cm2
 P at 4 hours
4.0x109
E 3.0x1 09
 2.0x109 
1.0x10 
1.0x109
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 selfforce. Growth of the loops can be
visualized from the change in the interstitial distribution near the layer
? 5x1013
E
0
C
o
0
I2x1013
12
_ 5x1012
Qvin12
 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 cm2
'  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 logscaled 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.
82
900C Oxidation
5 0 .0 , , ,1 , ,,,, ..L
 . .without dislocation loops
0
 with dislocation loops
C ( Si implant dose = 5x1015 cm2)
E
1.
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 fourhour oxidation, with and without the dislocation
loops.
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 loopinterstitial 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 8Simulation; 2xlo15cm2
0 1.4x1015 Data; 5x1015 cm2
S O Simulation; 5x1015 cm
1.2x1015
5 1.0x1015
.oxilo14
S8.0x1014
.0
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 cm2.
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 planview 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 noncircular 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 twodimensional
simplification can lead to an effective model for the interaction of extended
defects and point defects, which is an essentially threedimensional
phenomenon.
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 noncircular 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 cm2 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 equalsized loops led to a larger discrepancy at the four hour condition
[82].
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
5x1010
E
0 4x1010
S3x1010
Z 2x1010
C x
0D
900C
, , I , I ,
Oxidation
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
I I I I I
1
(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 cm2
E Simulation 2x1015 cm2
Data; 5x1015 cm2
 O Simulation 5x1015 cm2
Data; 2x1015 cm
Simulation; 2x1015 cm2
Data; 5x1015 cm2
 Simulation; 5x1015 cm"2
700
600
500 
400 
300
200
100
: : : : , ,
'
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 nontriangular 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 pn junction region and from causing detrimental leakage
currents. More systematic data extraction from TEM pictures will lead to
9000C Oxidation
,3000 
P :
S2500 
 2000
C
: 
2 1500
0
S1000
500
,)
n
Lp data
 Lp simulation
Rmax data
 O Rmax simulation
0 Rp simulation
 A Rmin simulation
Si implant dose = 2x1015 cm2
A A  A     A
^' 00~
_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. Asgrown
boron profiles were obtained by SIMS. The endofrange dislocation loops
were introduced into some of the samples by Si implantation at 50 keV to a
dose of 2x1015 cm2. 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.
Parameters
KIL (loopinterstitial reaction constant) @9000C
KVL (loopvacancy reaction constant) @9000C
a (capture and emission crosssection)
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
Value
4.3x103 sec1
1.0x104 sec1
7.5x107 cm
1.11 A
2.47 A
0.5
2.3
0.0
0.9
I
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
1019
1018
1017
1016
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
A
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.
90
 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
E
9
I 7
(1)
"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
E
m  Simulation: With loops; oxidation
o 100
g 80
60
CL
0 40
2 0 . i.i .
0 1 2 3 4 5
Time (hours)
Figure 3.22. Representation 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 represents 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 bestfit 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