Modeling of mechanical stress in silicon isolation technology and its influence on device characteristics


Material Information

Modeling of mechanical stress in silicon isolation technology and its influence on device characteristics
Physical Description:
x, 140 leaves : ill. ; 29 cm.
Rueda, Hernan A., 1970-
Publication Date:


Subjects / Keywords:
Integrated circuits -- Design and construction -- Simulation methods   ( lcsh )
Silicon-on-insulator technology   ( lcsh )
Electrical and Computer Engineering thesis, Ph.D   ( lcsh )
Dissertations, Academic -- Electrical and Computer Engineering -- UF   ( lcsh )
bibliography   ( marcgt )
non-fiction   ( marcgt )


Thesis (Ph.D.)--University of Florida, 1999.
Includes bibliographical references (leaves 133-139).
Statement of Responsibility:
by Hernan A. Rueda.
General Note:
General Note:

Record Information

Source Institution:
University of Florida
Rights Management:
All applicable rights reserved by the source institution and holding location.
Resource Identifier:
aleph - 030367677
oclc - 42843993
System ID:

This item is only available as the following downloads:

Full Text







Copyright 1999


Hernan Rueda

To my family


I would like to express my sincere gratitude to Dr. Mark Law,

chairman of my advisory committee, for his patience and guidance. He has

introduced me to the research area of process modeling, within which this

dissertation falls, and has provided valuable advice and direction throughout

my master's and doctorate studies.

Thanks also go to Drs. Gijs Bosman, Toshikazu Nishida, Kenneth K. O

and Kevin Jones for their interest and participation in serving on my

committee and their suggestions and comments.

I would also like to thank the University of Florida Graduate School

for its support during my master's program and the Semiconductor Research

Corporation for its support of my doctorate studies.

I have been very lucky to work with many industry mentors who have

provided invaluable assistance for my graduate studies. I acknowledge Drs.

Jim Slinkman and Dureseti Chidambarrao of IBM for their suggestions and

direction of the STI SKPM experiment. I also thank Dr. Len Borucki of

Motorola for assistance in design of the diode bending experiment. Thanks

also go to Drs. Paul Packan and Steve Cea of Intel for many valuable

discussions and suggestions on modern device concerns and strain modeling.

I wish to thank my friends who have made my time, over ten years at

the University of Florida, a very enjoyable experience. I've been very lucky to

meet so many good friends such as all the old TCAD research assistants, the

SWAMP group, my many roommates through all the years, and my 'old

school' friends back in the 'hood in Miami.

Last but not least, I express my love to my parents, Hernan Sr. and

Gloria, and my brother, Camilo, for their never-ending support, love and

encouragement throughout my whole life.



ACKNOW LEDGM ENTS................................................................................... iv

A B STR A C T ......................................................................................................viii


1 INTRODU CTION ......................................................................................... 1

1.1 M motivation ............................................................................................ 1
1.2 Stress-Induced Effects in Silicon Fabrication................................... 3
1.2.1 Oxidation Influences ................................................................ 4
1.2.2 Diffusion Influences ............................................................... 6
1.3 Stress-Induced Effects in Silicon Device Operation ......................... 7
1.3.1 Carrier Mobility Influences.................. .......... ............ 7
1.3.2 Energy Band Influences.............................................. ........... 11
1.4 G oals ........................................................... ................................ 12
1.5 Organization................................................................................. ....... 14


2.1 Continuum Mechanics ............................................................. 16
2.1.1 The Stress Tensor......................................................... 17
2.1.2 The Strain Tensor............................................................ 20
2.1.3 Stress-Strain Relationships ........................................ ........... 23
2.2 Strain Sources ....................................................... ................... 30
2.2.1 Film Stress................................................. ................... 31
2.2.2 Dopant Induced Stress ............................................ .......... .. 35
2.2.3 Oxidation Volume Expansion ....................................... .......... 39
2.3 Strain Computation Methods ............................................... .......... .. 40
2.3.1 Boundary Loading Method....................................... .......... .. 45
2.3.2 Fully-Integrated Method......................................... .......... ... 47
2.4 Sum m ary ...................................................................................... 48

3 APPLICATION EXAMPLES AND COMPARISONS............................... 50

3.1 FEM Com parisons......................................................................... 50
3.1.1 LO COS ................................................................................. 55
3.1.2 Post-STI Process Re-Oxidation............................................... 59
3.2 Raman Spectroscopy Measurements and Comparisons................. 65
3.2.1 Raman Simulation Method.................................... ........... 66
3.2.2 Nitride Film Edge-Induced Stress.......................................... 67
3.3 Boron-Doped Cantilever Bending Comparisons ............................... 70
3.3.1 Silicon Bulk Micro-Machining ............................................. 70
3.3.2 Cantilever Bending Simulations .......................................... 73
3.4 3D Boundary Loading Method Example (LOCOS)......................... 80
3.5 Sum m ary ....................................................................................... 82

4 KELVIN PROBE FORCE STI EXPERIMENT......................................... 85

4.1 Scanning Kelvin Probe Force Microscopy....................................... 86
4.2 Work Function Influence ............................................. ........... 88
4.3 STI Experiment.......................................................... ........... .. 93
4.4 SKPM Measurements ............................................................. 94
4.5 STI Strain Simulations............................................... .......... ..... 99
4.6 Results and Discussion .................................................................. 101
4.7 Sum m ary ....................................................................................... 105

5 STRESS INFLUENCES IN DEVICE OPERATION............................. 107

5.1 Uniaxial Stress Influence Experiment.......................................... 107
5.1.1 Stress-Inducing Apparatus Design..................................... 109
5.1.2 Stress-Wafer Deflection Relationship................................. 112
5.2 Experimental Procedure ................................................................ 114
5.3 Experimental Results .................................................................... 118
5.4 Sum m ary .......................................................................................... 124

6 SUMMARY AND FUTURE WORK ...................................................... 126

6.1 Sum m ary ....................................................................................... 126
6.2 Future W ork .................................................................................. 129
6.2.1 SKPM Strain Measurement Calibration Studies .......... 129
6.2.2 Effect of Stress on Dislocations............................................. 130
6.2.3 Silicidation-Induced Stress ................................................. 131
6.2.4 Three-Dimensional Modeling of the STI Process.............. 131

REFEREN CES........................................ ................................................ 133

BIOGRAPHICAL SKETCH........................................................................... 140

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



Hernan Rueda

May 1999

Chairman: Dr. Mark E. Law
Major Department: Electrical and Computer Engineering

One of the challenges the semiconductor industry faces as it attempts

and continues the scaling of silicon integrated circuits is understanding and

control of mechanical strain resulting from silicon fabrication technology.

High magnitudes of strain can be induced under standard fabrication

conditions and may produce deleterious effects in device behavior, such as

increased current leakage. Current leakage has been identified as a critical

device characteristic for future sub-micron dynamic random access memory

(DRAM) and complementary metal oxide semiconductor (CMOS)

technologies, as it is a limiting factor for increasing switching speeds and

decreasing power consumption. The following are various known sources of

stress in silicon technology: thermal expansion mismatch, intrinsic stress,

and oxidation volume dilation. This work results from an examination, by

modeling, experiment, and simulation, of the contribution of stress due to

these sources using the Florida Object-Oriented Process Simulator


The contributions of each source can be simulated using different

models that represent or approximate the physics involved. After the models

are described and presented, example applications are provided to

distinguish the advantages and limitations for each model.

Coupling experiment along with process simulation then validates the

results and allows for a better understanding of the problem. One such

problem examined in this work is the strain induced by the shallow trench

isolation (STI) process. STI has become an essential isolation scheme for

present and future sub-micron processes. It consists of several sequential

steps that exert stress in the silicon active area by each of the previously

described sources. Scanning Kelvin probe force microscopy (SKPM) is then

applied as a new technique to characterize the strain exerted from STI

processes through measurements of strain-induced work function variations

in silicon. Qualitative agreement is demonstrated between the SKPM

measurements and the work function influence due to finite element based

STI induced mechanical strain computations.

Finally, a wafer bending experiment is performed that quantifies the

influence of tensile and compressive uniaxial stress on forward current of pn-

junction devices. This effect is then modeled primarily through the strain

influences in the silicon bandgap.


1.1 Motivation

During the quest for increasing device density in integrated circuits,

many problems are encountered and need to be solved. As some problems are

alleviated, new issues emerge. One of the problems gaining importance in

silicon fabrication is process-induced mechanical stress. Many of the

processes used in silicon IC fabrication individually and cooperatively

contribute to the development of stress in the silicon active areas.

Of prime interest is the mechanical stress generated in the isolation

process flow. Isolation technology is continually challenged as design rules

are scaled. It is well known that high stress magnitudes in certain silicon

dioxide structures cause a decrease in oxidation rate [1, 2].

In LOCOS (LOCal Oxidation of Silicon) isolation processes, silicon

nitride is deposited over a thin pad oxide and patterned to mask the

oxidation of silicon. In order to cope with the necessary electrical isolation at

the ever-decreasing linewidths, increasing the stress during oxidation is

exploited [3]. This has been achieved by increasing the nitride thickness and

decreasing the pad oxide thickness. These techniques provide the sharp

transition from isolation field oxide to active area that is necessary in

decreasing pitch lengths for 0.35 micron technologies.

The main problem associated with this trend is that the silicon

substrate yields under the increased magnitudes of strain produced [4, 5].

Dislocations that are generated will degrade device performance [5, 6].

Dislocations serve as gettering sites that attract metal atoms introduced

during subsequent processes. Junctions are continually becoming shallower,

and therefore these nucleated dislocations have an even greater probability of

lying across the device junctions lined with metal atoms. Unacceptably

increased magnitudes of leakage current then result. Off-state currents are

critical design characteristics in logic and memory circuits, limiting the

switching speeds, causing increasing power consumption and limiting


Shallow Trench Isolation (STI) is steadily becoming the predominant

isolation technology as minimum feature sizes decrease below the minimum

attainable pitch lengths of LOCOS-based technologies. The general STI

process flow includes a nitride-pattered reactive ion etch, sacrificial sidewall

oxidation, oxide deposition and finally a chemical mechanical polish.

However, stress-induced dislocation generation is not exclusive to LOCOS

based isolation technologies and is also a factor in designing STI processes

since the steps in the STI process may cooperatively strain the enclosed

silicon active areas [7].

Aside from isolation processes, other fabrication processes also induce

strain in the silicon crystal. Processes such as thin film deposition and

dopant introducing processes induce strain that can generate dislocations. As

feature sizes decrease, all these different strain-generating sources become

closer together and their influences overlap each other. Under these extreme

circumstances, it is important to understand strain fields generated by

multiple sources neighboring the silicon active areas.

Strain fields may also play a role in affecting process and device

behavior other than generating dislocations. For example, in silicon micro-

machining, high residual stresses of boron induce bending sensor structures

composed of diaphragms and cantilevers [8]. Strained regions have also been

shown to affect diffusion of dopants [9]. Also, it is well known that crystal

strain affects the energy band structure [10] and carrier mobilities [11]. Both

are major parameters influencing the electrical properties of a device. These

are just a few of the many concerns related to the strain state of the crystal.

1.2 Stress-Induced Effects in Silicon Fabrication

Stress concerns in process design first became significant in LOCOS

isolation technologies. As thermal oxidation temperatures are reduced,

dislocation densities in silicon increased due to this isolation technology [4].

These effects were attributed to the increased stress magnitudes induced at

lower oxidation temperatures due to the increased viscosities of the silicon

dioxide. Since then, the correlation between increased dislocation densities

and high stress magnitudes have caused isolation technologies to be a major

stress related concern.

Since then the diffusion process has also been reported to be influenced

by film-induced stress [12]. This affect has been explained by stress

influences on point defect concentrations as well as extended defect size and

concentrations such as dislocation loops. As junction depths continue to

decrease in the evolution of sub-micron technologies, the stress influences

generated by surface films have a greater effect.

1.2.1 Oxidation Influences

It is well known that the growth rate of thermal oxidation of silicon

has a stress dependence [13]. This stress dependence occurs in nonplanar

regions. In thermal oxidation, a volume of silicon reacts into a volume of

silicon dioxide that is 2.2 times in volume. For planar oxidation this does not

induce a large stress. This is due to the newly grown oxide lifting the old

oxide perpendicular to the surface. Since the surface of the oxide is not

constrained the oxide is free to move in this direction so negligible stress is

induced normal to the Si-SiO2 interface.

In nonplanar regions such as convex and concave corners of silicon, a

strain is exerted in the silicon dioxide. For the convex corners the strain is

laterally tensile since the old volume of oxide is stretched around a longer

periphery. In concave corners the strain is laterally compressive since the old

volume of oxide is compressed into a smaller periphery. Convex corners occur

at the top corners of a trench oxidation and concave corners result at the

bottom corners of a trench oxidation and at the Bird's Beak in a LOCOS


It was reported through cylinder structure oxidation experiments that

the stress induced in these structures reduced the amount of oxide grown as

compared to planar oxidation [1, 2, 13]. This behavior has been observed in

highly stressed LOCOS processes and in corners of trench oxidation. It has

been documented that the stress induced in the oxide alters the oxidation

reaction rate, oxidant diffusion of reactant, and the oxide viscosity [2]. All

these influences affect the oxidation growth rate.

The stresses generated in the silicon dioxide are relaxed by exerting

force on the neighboring films and the silicon substrate. This leads to

strained regions in the silicon due to the oxide growth. Silicon is an elastic

material for a wide range of strain. However, as the induced stress exceeds

the critical yield stress, dislocations in the silicon crystal are generated to

relax the strain [14]. The generated dislocations then present problems in

device behavior. Therefore, the strain induced in the silicon by oxidation

becomes a concern.

Trenches are popular isolation technologies that also exhibit

dislocation generation problems [7, 15]. The stress induced by trench

structures is due to other sources also. After the reactive ion etch, an

oxidation is performed to provide a low interface state density and a low

oxide trap density. This oxidation is the first source of stress in the silicon

substrate in the trench fabrication process. Next the trench is filled with a

deposited film. This film influences a stress due to its built in intrinsic stress

and thermal expansion mismatch. The three sources all act simultaneously

influencing strain in the surrounding silicon substrate. Again, the primary

concerns are dislocations generated to alleviate the strain in the substrate

induced by the trench process.

1.2.2 Diffusion Influences

It has been observed that both phosphorus and boron diffusion

behavior under silicon nitride films is different than in inert conditions [9]. It

has been established that boron and phosphorus diffusion is governed by

interactions with point defects, namely silicon interstitials. Ahn attributed

the measured diffusion reduction to a vacancy supersaturation and self-

interstitial undersaturation that may be due to the compressive stress under

the silicon nitride films [12]. Osada later confirmed this by performing

similar experiments [16].

These experiments were performed with junction depths of

approximately 0.8 microns. Scaling for sub-micron devices has lead to

decreasing junction depths. The magnitudes of strain in silicon due to film

stress are at a maximum at the film interface. Therefore, it is believed that

an even greater effect is encountered for the shallower junctions that are

getting more prevalent in modern technologies.

1.3 Stress-Induced Effects in Silicon Device Operation

Interest in mechanical properties of silicon was sparked with the

discovery of its piezoresistive effect [17, 18]. This led to the increased interest

of the use silicon as a pressure sensor. Sensor research also extended

investigation of the mechanical influences on other electrical parameters

such as the energy bands of semiconductor crystals [10, 19-22]. As a result,

there is satisfactory understanding of how applied mechanical forces affect

electrical device behavior. These principles are used mainly in the design of

semiconductor mechanical sensors. In the present age mechanical stress is

becoming a limitation in silicon device fabrication, these principles also need

to be incorporated into microelectronic device behavior.

1.3.1 Carrier Mobility Influences

The influence of mechanical stress on carrier mobility is starting to

gain importance in modern microelectronics. It has been shown that in

CMOS transistors the transconductance will vary dependent on the

magnitude of applied stress and the behavior is dependent on crystal

direction of the current flow as well as the type of carrier [23]. This effect is

was also noticed and believed to be due to LOCOS-induced stress in SOI

CMOS devices [24]. The variation in transconductance is attributed to the

piezoresistance effect on the carrier mobilities.

The piezoresistance effect of semiconductors describes how the

resistivity is influenced by mechanical stress. The electric field vector is

proportional to the current vector by a symmetrical resistivity tensor of rank

two with nine components:

El pA P4 P6 1
E2 P4 2 5 2 *
2 = p, p2 p (1-1)
E3 P6 P5 P3_ Ji35

If the system axis is aligned along the <100> crystal directions, the

normal resistivity components p1, p2, and p3 relate the e field vector

component to the current vector i component in the same direction. The cross-

resistivity components p4, P5, and p6 relate the E field vector component to the

current vector i component in a perpendicular direction. If in the unstressed

state, the normal components have the same magnitude p and the cross

components are equal to zero, this reduces to the following isotropic


E = pi. (1-2)

When the crystal is under mechanical stress, the resistivity

components change as the following:

P1 p App
P2 P Ap2

P3 P AP3
=3 + A3 (1-3)
P4 0 AP4
P5 0 AP5
P6 0 Ap6

The piezoresistive coefficients relate the stress-induced changes in the

resistivity components to the stress tensor influencing the change. This

matrix relating the six resistivity components to six stress components

consists of 36 piezoresistive coefficients nij. Due to the cubic crystal

symmetries in silicon, the piezoresistance coefficient matrix reduces to three

independent components, Rn1, 7C12, and x44:

Ap, ~n t12 X2 0 0 0 a
Ap2 2 711 r 12 0 0 0 02
1 Ap3 7_ 2 12 1 0 0 0 03 (
*- (1-4)
p Ap4 0 0 0 4 0 0 a4
Ap5 0 0 0 0 7C,4 0 5
_Ap6 0 0 0 0 0 7r44 0C6

The stress components are also referenced with the system axis oriented in

the <100> directions. Smith initiated the investigations of these

piezoresistive coefficients and found the following values for silicon at room

temperature displayed in Table 1-1 [17, 18].

The piezoresistance coefficients are also dependent on dopant

concentration as well as temperature. Later it was found that they would

decrease as the temperature increases and/or the dopant concentration


Table 1-1: Piezoresistive coefficients for silicon [17, 18].

Material p 711 7X12 744

(4-cm) (10-12 cm2/dyne)

p-Si 7.8 +6.6 -1.1 +138.1

n-Si 11.7 -102.2 +53.4 -13.6

The accepted explanation for this phenomenon is the many-valley

model [17, 18]. Anisotropic conditions exist when the mobility in one crystal

direction is different than the mobility in the other crystal lattice directions.

This results when the semiconductor is in a stressed state. The stress tensor

distorts the conduction energy bands of the unstressed semiconductor in

different magnitudes depending on direction. The energy levels and

curvatures of the band energies corresponding to the perpendicular directions

are influenced differently by the applied strain. The effective masses of the

carriers are proportional to the energy bands' curvatures in reciprocal k

space. Since the carrier mobilities are functions of the carrier effective

masses, the strain influence on the energy band level curvatures results in

directionally dependent influences on the carrier mobilities and therefore the

resistivities of the semiconductor. The energy band shifts are also influenced

on dopant concentration and temperature. Therefore the energy band's

sensitivity to stress will also be dependent on these influences.

1.3.2 Energy Band Influences

The mechanical stress state's influence on the energy bands also

affects the electrical behavior of p-n junction devices such as diodes and

bipolar transistors. In these devices the operation is governed by the flow of

minority carriers. Using a diode as an example, the forward bias current is

described by the following relation:

IF = I exp() + IRO exp( ) (1-5)

where the saturation current is

I,= qA po coth + Dn coth (1-6)
L L" L.

and the recombination current is

IRo -qAn (1-7)

The saturation current term is linearly related to the minority carrier

densities npo and pno. The minority carrier densities are directly dependent to

the square of the intrinsic carrier density

no =- (1-8)

The intrinsic carrier concentration is exponentially dependent to the stress

dependent band gap Eg

ni =KT 3/2 exp( -g). (1-9)

The stress induced shifts in the conduction and valence energy bands will

alter the band gap and therefore ultimately influence the saturation current.

Wortman initiated the quantification of the effects of uniaxial and

hydrostatic compressive applied external stresses to forward and reverse

biased diodes [10, 22].

Mechanical stress also influences the generation/recombination

current component of p-n junction devices. Rindner attributed the effect of

uniaxial compressive stresses to increased dislocation densities that

decreased the carrier lifetimes and therefore increased the

generation/recombination current component [25]. This effect becomes the

greater influence under higher magnitudes of stress due to the dislocation

generation to relax the applied stress.

1.4 Goals

The goals of this work are primarily to develop a system where strain

can be computed from multiple sources simultaneously. Silicon IC fabrication

involves a sequential flow of many processes that introduce and alter the

strain in the crystal. These process-induced strain fields influence the

behavior of processes later in the fabrication flow as well as device operation

once the process flow is completed. An accurate strain solution is necessary

for further investigation of its effects.

Once the strain in the system is understood, strain dependent models

can be developed to help understand unexplained behavior that has been

observed that may be due to strain. Such areas may include point defect and

extended crystal defect interactions, diffusion kinetics, and band-gap and

mobility influences.

Strain simulations also could aid in the development and analysis of

isolation process technologies. Often in fabrication process development, the

stress levels generated and dislocation densities produced may decide the

isolation process required. A strain field simulator could reduce the amount

of experimental work necessary for solving these problems.

Another goal of this work is to validate the process induced strain

models with experimental measurements. Currently, this is a major hurdle

due to the few characterization techniques available for localized sub-micron

strain measurement. Micro-Raman spectroscopy has been the most suitable

method for investigating localized strains [26]. But as technologies continue

to scale towards the 0.1pm generation, this may even surpass micro-Raman's

spatial resolution limits. Therefore, in this work Scanning Kelvin probe force

microscopy (SKPM) has been investigated as a new method for analyzing

localized strains through detecting influences in the silicon work function.

One last goal is to quantify the influence of tensile and compressive

stress on pn-junction device current. This would then provide some insight

into the mechanical strain influence on leakage currents.

1.5 Organization

Chapter two provides descriptions of the various process induced strain

sources and discusses how they are modeled in this work. Afterward, the

finite element methods that were developed for strain computation are then


Results and comparisons between the methods for various processes

are included in chapter three. Example applications are provided to

distinguish the advantages and limitations for each model. Next, the film

edge-induced stress solutions are validated with published micro-Raman

measurements. Finally, three-dimensional applications are demonstrated.

Scanning Kelvin probe force microscopy (SKPM) is then explored as a

technique for characterizing STI induced strain in chapter four. An STI

experiment is performed and strain influence is measured by SKPM. These

measurements are then compared with simulations of band gap influences

using the models described in chapter two.

An experiment relating mechanical stress to pn-junction device

operation is then described in chapter five. This wafer bending experiment

addresses uniaxial stress influences on the forward current. This allows for


quantifying the influence on the reverse leakage current through observation

in the forward bias.

Finally chapter six concludes with a summary of the research work

accomplished and addresses topics for future work.


The solution of strain present in a particular device technology is

computed using a finite element method (FEM) formulated to solve for the

strain induced by various sources in silicon technology [27, 28]. In a

fabricated device, the strain field in the silicon is generated due to various

processes at different steps along the fabrication flow. The most critical

sources for inducing stress are deposited and grown films. Sources in the

silicon crystal such as dopants and extended defects are becoming more

important as technologies advance.

In this chapter, a brief review of continuum mechanics is first

provided. Next the individual strain sources are then described. The

algorithms used to model the strain generated from each source are also

discussed. Afterward the finite element methods used to integrate the

various strain sources are described.

2.1 Continuum Mechanics

It is the intention that this review refreshes the reader with the theory

and notation of continuum mechanics. References are provided for a more

complete description. The stress tensor is first introduced in this section.

Next the strain tensor is described. Finally this section concludes with

different stress-strain relationships descriptions and how they may relate in

silicon processing.

2.1.1 The Stress Tensor

Stress is the distribution of internal body forces of varying intensity

due to externally applied forces and/or heat [29]. The intensity is represented

as the force per unit area of surface on which the force acts. To illustrate this

concept, consider an arbitrary continuous and homogeneous body (Figure 2-1)

under the applied external forces, F1, F2, F3, and F4. If the body is sliced into

two smaller volumes V1 and V2, then V2 exerts force on V1 at their surface


+z +Y V

+x S

Figure 2-1 Continuous body with external forces applied.


interface S to remain in equilibrium. The resulting force may be of varying

intensities along the surface.

At any point P on the surface between V1 and V2, the forces can be

reduced to a force and a moment, which can be described by a stress vector

acting on that surface. Three stress vectors acting on three mutually

orthogonal planes intersecting at that point can then determine the stress

state at any point P (Figure 2-2). The stress tensor is composed of the three

stress vectors and, according to Cauchy's equations of motion, is sufficient to

define the stress state in any element in a body [30].

To illustrate the tensor nature of stress present at point P in the

continuous body, consider a cubic element (Figure 2-2) of infinitesimal


+z Y" t 0

I /


Figure 2-2 An infinitesimal cubic element located within a continuous body
with stress tensor components shown.


dimensions located at point P in Figure 2-1. For simplicity of notation, let the

cube be aligned perpendicular with the system axis. The stress vector Tx

acting on the plane normal to the x-direction is the following:

T, = a, + f-+- (2-1)

Let the surface AS be the plane of the cube normal to the x-direction. The

stress vector Tx is defined as the ratio of force acting on that surface area AS:

T = lim (2-2)
s -4 ASO dSX

The stress tensor on that volume is defined by nine stress components acting

on the three surfaces of the cubic element, which make up the three stress

vectors Ti:

9,= x T xyy Ty Z (2-3)

In the definition above, ii are the normal stress components acting on the

faces perpendicular to i-direction and zj are the shear stress components

oriented in the j-direction on the face with normal in the i-direction. At

mechanical equilibrium, it can be shown that the stress tensor is symmetrical


Tij = Tj .


A column vector of six independent components can then describe the state of

stress at a point:

UT =[UXX ayy U z, Tu Ty ]. (2-5)

2.1.2 The Strain Tensor

The application of stress to a body in equilibrium causes it to undergo

deformation and/or motion. A measure of deformation is strain. Two common

measures for strain are the Lagrangian and the Eulerian definitions. Both

are functions of the initial and final dimensions. When the displacement

between the final and initial measurement is referenced to the original

position dimensions then it is known as the Lagrangian definition. The

Eulerian definition describes the deformation displacement referenced with

respect to the deformed position.

Figure 2-3 shows a one-dimensional example of the deformation in a

(a) (b)
Unstretched Stretched

Ax t Ax + Au v

Figure 2-3 One-dimensional deformation of a spring: (a) original length
(Ax), (b) deformed length (Ax+Au).


spring. The Lagrangian strain and Eulerian strain then, respectively, become

the following over the length of the spring:

increase in length Au
E- (2-6)
S original length Ax

increase in length Au
Ei (2-7)
deformed length Ax +Au (

For the case of infinitesimal deformation, the Eulerian and Lagrangian

descriptions become equivalc. ,. An infinitesimal deformation approximation

requires that the maximum deformations involved be much smaller than the

smallest dimension of the deformed body. A second requirement is that the

deformation gradient is much less than one. When these assumptions hold,

the infinitesimal strain can be defined by the following relationship for any

point in the spring [31]:

A. Au Qu
E=lim -(2-8)
Ax-o Ax ax

By expanding this definition in three dimensions, the strain is related to the

displacements by the following strain components [30]:

Sdu (du v
E =- E = E +--

E dv v = E = (2-9)
e", C =i +-) (2-9)

dw 1(du dw
cE E E i-+-
U U x 2 dx

where u, v, w are the displacements in the x, y, and z directions, respectively.

These components make up the strain tensor (s i) that is analogous to the

stress tensor:

Ek = yx

xy 'xzl
Eyy Cyz
SZ zz


A three-dimensional example of a body undergoing deformation due to

an externally applied force is illustrated in Figure 2-4. In the example a

compressive force is applied to the infinitesimal cube in the x-direction. A

negative displacement compressivee strain) results in the x-direction and

positive displacements (tensile strains) result in the y- and z- directions. The


/Zf 7\


+z +


Compressive Force in
x direction

S----- du +
i du/dx < 0
dw/dz > 0
dw +
------dw/dz dv/dv > 0
<- dv + dv/dy -


Figure 2-4 Strain reference example displaying a compressive stress in the
x-direction that generates normal strains in the x-, y-, and z-directions.


relationship of how the strain in each dimension results from the applied

force will be discussed next.

2.1.3 Stress-Strain Relationships

Bodies of different materials but of same dimensions may deform

differently under the same stress application. The relationship between the

stress tensor and the deformation is known as a constitutive relation. It may

vary for a given material depending on conditions such as temperature and


Elasticity. All structural materials possess, to some extent, the

property of elasticity. Elastic bodies possess memory during deformation. For

example, when a force is applied on an elastic solid, it will deform until it

reaches its elastic yield limit or until the load is released. Microscopically the

bonds between atoms that compose the solid 'stretch' during elastic

deformation. When the force is removed the body will return to its original

shape if it is an ideal elastic body and it had not reached its yield stress,

similarly to an ideal spring. When the load is removed the bonds return to

their unstressed lengths corresponding to the original environmental


For a Hookean elastic solid, the stress tensor is linearly proportional to

the strain tensor over a specific range of deformation:


T = CqklE k


by the tensor of elastic constants cykl [30]. In order to relate each of the nine

elements of the second rank strain tensor to each of the nine elements of the

second rank stress tensor, ijki consists of a fourth rank tensor of 81 elements.

However due to the symmetries involved for the stress and strain tensors

under equilibrium, ciki is reduced to a tensor of 36 elements.

Crystal silicon has diamond cubic crystal geometry resulting from its

strong directional covalent bonds. For such crystals, cijki has the following

form due to their cubic symmetry [32]:

Ct, c12 c12 0 0 0
C12 CtL C12 0 0 0
c,2 cL2 C1, 0 0 0
c = (2-12)
C*j =0 0 0 c4 0 0 ( )
0 0 0 0 c, 0
0 0 0 0 0 c4

Thus, for silicon the tensor of elastic stiffness constants reduces to the three

independent components: c1, c12, and c44. Due to the crystal's lattice

temperature dependence, the elastic constants are also thermally dependent.

At room temperature (25C) the elastic constants have been measured

as the following for silicon [33, 34]:

Ct, = 1.657 x 102dyn/cm2
C12 = 0.639 x 1012 dyn/cm2
c. = 0.7956 x1012 dyn/cmr2.
The elastic constants' thermal dependence has been documented as the

following linear relationship for silicon [35]:

Tc1 -L =-75X106 /0K
Tc12 -At2 =-24.5x106 /Ko
Tc A4 -55.5 x 10 /K.
From the above linear dependence, it can be seen that the elastic constants

do not change significantly for large temperature changes. Also the degree of

anisotropy does not change significantly for the range of 100-9000K [36].

These studies support that silicon acts as an anisotropic elastic material over

a wide temperature range frequently encountered in silicon IC fabrication.

Although silicon is an anisotropic crystal, it is sometimes desirable to

approximate it with isotropic elastic properties for simplification. When the

components of elastic constants for a material are equal for any rotation of

the reference axis, the material is said to be isotropic. This means that the

elastic properties of the material are the same in all directions. The tensor of

elastic constants for an isotropic material reduces to the following:

E(1- v)
cL = (2-13)
(1+ v)( 2v)
cL2 = (2-14)
(1 + v)(1 2v)
c ( = v (2-15)
(1+ v)
where E represents the Young's modulus and v represents the Poisson's ratio.

As was demonstrated in Figure 2-4, contractions in one dimension may be

accompanied by dilations in other dimensions. The Young's modulus is a

measure relating stress and strain for an elastic material when stress is


applied in one direction and the other directions are free to deform as in

Figure 2-4. The Poisson's ratio is a material property describing the ratio of a

strain perpendicular to the applied stress to the strain oriented in the

direction of the applied stress.

Due to its anisotropy, E and v vary depending in the direction of the

applied stress and plane acted upon for silicon. At room temperature E may

vary from 1.3e12 dyn/cm2 (for <100> directions) to 1.875e12 dyn/cm2 (for

<111> directions). Poisson's ratio also may vary from 0.06 to 0.34 for the

same conditions. The following are measured values for E and v in silicon at

room temperature [32, 33, 35]:

E[ooI =1.31 x 102 dyn/cm2
E11o1 = 1.69 x 10L dyn/cm2
EI L1 = 1.875 x 102 dyn/cm2
voo = 0.279.
Viscosity. Although silicon deforms elastically over a wide temperature

and load range, other materials used in silicon fabrication behave differently.

Some deform elastically at temperatures near room temperature and flow at

higher temperatures with fluid behavior. Silicon dioxide (Si02) and silicon

nitride (Si3N4) are examples of this these type of materials. Fluids resist

deformation with viscous behavior.

As was previously mentioned, Hookean elastic solids will return to

their original shape after an applied stress is removed. Materials with

viscous constitutive properties may not. Viscous bodies relax or minimize the


strain associated with an applied stress. As the strain is relaxed the stress

also reduces. Due to the reduction in strain, when the applied stress is

removed, the body will retain its currently deformed shape.

Microscopically, in bodies with viscous mechanical properties, the

bonds between atoms that compose the solid break during deformation. New

bonds are then formed and re-broken as the body deforms under an applied

force. When the applied force is removed the bonds remain in their current


A common constitutive relationship for a viscous body is the

Newtonian fluid. In a Newtonian relationship, the shear stress on the surface

is linearly proportional to the rate of deformation [30]:

Crj oc pkijki (2-16)

where /jki is the tensor of viscosity coefficients. The components of ijkl for

fluids are not as well known as cijki for elastic solids. However, most fluids

appear to behave isotropically. An isotropic approximation with the

restriction of incompressibility (constant density) allows for the following

constitutive relationship known as a Stokes fluid [30]:

aj = -pS, +pij (2-17)

where p now is a scalar that represents the viscosity of the fluid. The normal

stress components are dependent on the static pressure p of the fluid where

y& represents the Kronecker delta function.


Viscoelasticitv. It has been recognized that some materials deform

with a combination of elastic and viscous properties. There are various

models that have been formulated to describe the mechanical behavior of

viscoelastic bodies. In the Maxwell model of viscoelasticity, the total strain is

simply the sum of the strain due to elastic deformation and the strain due to

viscous deformation:

E = EE +C. (2-18)

This relationship can then be expanded to formulate the well-known

Maxwellian viscoelastic relationship between stress and strain:

e=--+-- (2-19)
E p

where E is the elastic modulus and g. is the viscosity of the material. Figure

2-5 illustrates the differences in deformation responses among a Hookean

elastic, Newtonian viscous and a Maxwellian viscoelastic material given the

same applied pulsed stress. Notice that the viscous and viscoelastic response

are time-dependent. Instantaneously after the stress is applied, the

viscoelastic body deforms elastically. Later when the applied stress is

constant, the viscoelastic body begins to display a linear viscous deformation.

As the applied stress is removed, the viscoelastic body then elasticallyy'

deforms towards its original strain state. However due to the viscous

relaxation component, it does not return to its original strain state.


Therefore, for the Maxwellian viscoelastic relationship, the short-term

deformation is elastic and the long-term deformation is viscous.

Applied Stress

t2 t

Hookecn Elcstic Solid

't1 't2

Newtonian Viscous Fluid

ET Mcv wellicn Viscoelastic Fluid

Figure 2-5 Comparison in loading response among Hookean elastic,
Newtonian viscous and Maxwellian viscoelastic relationships.

e, = co(t,-t,)/-


In silicon IC fabrication, oxide (SiO2) is considered to behave

nonlinearly viscoelastic at midrange (800-11000C) anneal temperatures [2,

14]. The nonlinearity is due to the fact that the degree of its viscosity is also

stress dependent. Nitride has been recognized to behave as a viscous body in

this same temperature range [37].

Not very many materials behave exactly as a Hookean elastic solid, a

Newtonian viscous fluid, or a Maxwellian viscoelastic fluid. However in

limited ranges of stress, strain, and temperature, these constitutive relations

can approximate their deformation behavior quite well. As an example silicon

behaves elastically for a wide temperature and stress range. However under

high stresses, silicon will yield and nucleate dislocations and defects in the

crystal to relax the stress present. Therefore it important to learn under what

conditions this will result.

2.2 Strain Sources

The individual strain sources included for strain computation are

discussed next. Three different strain sources are modeled and integrated:

film-induced, dopant-induced, and oxidation volume expansion induced

strain. Each strain source is discussed and explanations of how they induce

strain in the bulk are included. Following each is a discussion on how the

strain induced is modeled using the finite element method.

2.2.1 Film Stress

Physics. As silicon technology progresses, layers upon layers of

different materials are grown and deposited on top of and adjacent to each

other. When materials that have different structural, mechanical and

thermal properties are attached to each other, strains in each of the

materials can result. Adjacent material films relax or expand differently

based on their material properties. This causes one film to stretch or contract

the other film in a manner that will cause a local strain in each film. Large

localized stresses can result due to discontinuous films in regions such as at

the film edges and in non-planar sections.

According to Hu [14], stresses result in thin films due to two different

mechanisms. The first is referred to as an 'extrinsic' stress and is primarily

due to thermal expansion mismatch of neighboring materials. The process

used to deposit the film is done at an elevated temperature ranging from 150-

12000C. After the process is over and the thin film is deposited, subsequent

thermal cycles will cause the film to expand and contract. If the film was

attached to or restrained by a rigid material that does not thermally deform,

the amount of strain induced is linearly proportional to the temperature


E, = a,, AT (2-20)

where acth is the linear thermal expansion coefficient of the film material. The

body experiences an increase in volume in each normal direction. No shear


components result from the thermal difference. For many materials the

thermal expansion coefficient is not necessarily constant or linear over a wide

temperature range. For silicon, ath varies 2.5-4.5x10-6 /K over the 300-900K

temperature range. Oxygen content and dopant concentration are factors

leading to ath variations.

If the material that the film is attached to also expands due to a

thermal increase then the local strain at the interface that is produced is due

to the difference in thermal expansion coefficients:

,t = (a,i -a,,).AT. (2-21)

Thermal mismatch stress is often incorrectly shortened and referred to

as thermal stress. However, thermal stress is due to thermal gradients

within a material. This often occurs in the wafers during before and after

temperature cycles. As the wafers cool down, the maximum stress is due to

surface tension. However to maintain force equilibrium the interior of the

wafer must be in compression. As the wafer heats up, the surface proceeds to

expand due to thermal expansion. However, the cooler interior of the wafer

restricts this expansion causing a compressive stress at the surface and a

tensile stress in the interior. Thermal stress primarily occurs in the substrate

as its thickness allows for a greater thermal gradient between the surface

and the interior. The thermal stress in the substrate becomes a problem

during temperature ramp rates encountered in Rapid Thermal Annealing

(RTA). At higher temperature ramp rates, the thermal stress built up from


the high temperature gradient may exceed the yield stress and cause the

wafer to shatter. Normally the thicknesses of grown and deposited films in IC

fabrication are so thin that a negligible thermal gradient exists across them.

The other source of stress encountered in thin films is the 'intrinsic'

stress. Several researchers attribute this stress as due to growth mechanism

of the material during the process. For grown oxides, the intrinsic stress

results from the planar volume expansion resulting from the oxidation

reaction. This stress should not be confused with the stress induced at

isolation edges, which is non-planar and is a multidimensional problem

discussed later. Other material films such as polysilicon, silicon nitride, and

silicides exhibit intrinsic stress also.

After a film is deposited or grown, the wafer will warp according to the

total stress in the film. A highly tensile film will bend the wafer's edges

towards the film and the reverse for a compressive film. A popular technique

for measuring film stress is to measure the amount of wafer curvature

optically. The total film stress is then proportional to the radius of curvature

by the following relation [38]:

E t 2
= -.' (2-22)
f 6(1-v,)R t,

where Es and Vs are the elastic properties of the substrate and t. and tf are

the thicknesses of the substrate and film respectively. The film stress

measured is the total stress due to thermal mismatch and its intrinsic


components. The intrinsic stress can then be derived from this measurement

and the known thermal mismatch stress from the previous relationship.

Even though high stresses can result from the sum of thermal

mismatch and intrinsic stress in a film, if the film is uniformly planar then

the stress resulting in the substrate due to the film will be orders of

magnitude smaller. This is due to the large difference in thicknesses between

the film (tf) and the substrate (ts):

a, = -4Q f -. (2-23)

Higher local stresses result in the substrate due to discontinuities in the film.

Such discontinuities include etched film edges from masking and

nonplanarities as in trench fill depositions.

Film-induced Strain Model. The stress due to deposited films is

modeled as an initial condition before deformation. The planar

measurements of intrinsic stress are used as the initial condition for the

finite element solution. The stresses are input at the nodes as the film is

deposited. They are directed in a biaxial tangential orientation along the

growing interface as is shown in Figure 2-6. To handle the stress components

in nonplanar interfaces, the stresses are translated from the planar system

axis to the axis perpendicular to the normal of the growing film. After the

film is deposited or grown, the stresses at the nodes are then averaged to

their neighboring triangular (2D) or tetrahedral (3D) elements.


The thermal mismatch stress is modeled as a hydrostatic stress. Each

normal component of strain exerted is equal in magnitude. No shear

components result from thermal mismatch. For each element in the film, the

three normal components of strain are added by superposition to the previous

state of stress due to other sources. This presents a problem in plane strain

FEM formulations since the strain in one direction is set to zero. This

problem will be addressed in section 2.3.

2.2.2 Dopant Induced Stress

As dopants are introduced to the silicon substrate, the mechanical

state of the substrate will change. The dopants may substitute for the silicon

positions in the lattice. Silicon atoms are displaced forming extended defects

that are lodged in the crystal lattice. Different dopants have varying atomic

sizes and therefore have different mechanical behavior in the crystal.

Precipitates and other atoms present such as oxygen and carbon also alter

the mechanical properties of the crystal.

Boron. Boron is well known as a substitutional dopant. Its atomic size

is smaller than that of silicon. When it locates into a substitutional site,

lattice contraction results due to its smaller size. This presents an atomically

localized strain in the crystal due to each boron atom. Figure 2-7 exemplifies

this in three- and two-dimensional illustrations. For high concentrations of

boron in silicon, these atomic strains add up significantly and result in a

larger localized region of strain in the boron doped silicon lattice. The non-


boron doped silicon region will resist the diffused boron layer from

contracting and thus result in a tensile strain field. This effect has been

demonstrated in silicon micro-machining applications [8, 39, 40], where

boron-doped cantilevers have exhibited bending due to strain induced by the


2D planar film

These nodes' stresses are (7
oriented parallel to the -
deposition interface
-^-7- ^--Y / D --7, -

2D nonplanar
film nx'yC,

n, *,+ n,, *---- x ,\ ^ y

These nodes'
stresses are
around covers
by the no* n,* T,-
component of
the unit surface n,*oax *
normal Y

Figure 2-6 Intrinsic film stresses are oriented parallel to the interface on
which the film is grown or deposited.

2D approximation
0 0 000 00 00 0
0 0 0 0 0 0 0000

o o o o O o oOo o
00000 00000
0 0 000 0 0000

Figure 2-7 Two- and three-dimensional illustrations of lattice deformation
due to boron substitution.

Densitometric studies have been done with boron doped silicon

crystals. Horn measured the silicon lattice constant variation as a function of

boron concentration [41]. From his measurements, the induced strain (Aa/a)

was extracted and given as function of boron concentration (Figure 2-8).

Boron-induced Strain Model. The empirical relationship introduced by

Horn is used as the contribution of strain in silicon due to substitutional

boron dopant. According to his measurements, 0.0141 A of lattice contraction


results per percentage of B in Si at room temperature. Using this figure, the

following relationship is derived:

Aa 0.0141 C8
E, = E= = Eu -100
asi as, Ns,


where asi is the silicon lattice constant (5.4295A at 25C) and Nsi is the

density of Si atoms (5.02e22 cm-3). The average concentration of Boron (CB) is

computed for each element from its node quantities. Dopant-induced strain is

also hydrostatic (as in thermal mismatch strain) and again presents a

problem in plane strain FEM formulations since the strain in one direction is

set to zero. This problem will be addressed in section 2.3.

E 0.0005
g 0.0003
9 0.0002

0 5E+19 1E+20 1.5E+20
Boron Concentration (cm-3)


Figure 2-8 Hydrostatic strain as a linear function of boron concentration


Other dopants. The strain contributions of other common dopants such

as phosphorus and arsenic are not as great per atom as that due to boron.

Because arsenic has a much larger atomic mass than silicon, one would tend

to believe that arsenic would induce a large compressive strain in the silicon

lattice. However, it has been reported that heavily doped arsenic (5x1021 cm-3)

only induces a lattice compression (Aa/asi) of approximately 0.0019 [42]. The

atomic size of phosphorus is more closely matched to that of silicon and

therefore the P-Si bond lengths are of the same magnitude as the Si-Si bond


2.2.3 Oxidation Volume Expansion

The oxidation process also induces strain in the substrate due to the

net volume expansion of the oxidation reaction. It is well known that silicon

oxidizes to a volume of oxide that is 2.2 times larger. For planar oxidation,

this presents no problem since the newly acquired volume pushes the old

oxide upward towards the unconstrained surface perpendicular to the

interface. However, in nonplanar regions such as in trench corners and in

constrained regions such as in LOCOS edges, this presents a problem. For

these regions, the boundaries are constrained and therefore the newly

acquired oxide volume compresses against the earlier grown oxide. Since the

oxide has no place to move, large compressive strains build up in these

regions. The strains are somewhat relaxed by applying pressure to the silicon


interface. The forces applied to the silicon are high enough to surpass the

yield stress of silicon.

2.3 Strain Computation Methods

Two different finite element based methods are developed to solve and

compute the previously described models of the various strain sources. These

algorithms are developed in the process simulator FLOOPS [43] and are

integrated with and derived from methods developed to model stress-

dependent oxidation and silicidation [44].

Newton's second law of motion governing deformation is stated as the

following [30]

pa, =d + b (2-25)

where p is the density of the body, ai is the acceleration, oiy is the local stress

tensor, and bi is the body force.

The equivalent nodal force for each element (qe) may be represented as the


q' = BTad(vol) J Nbd(vol) (2-26)
V' V,

where b represents local body forces (e.g. gravitational or electromagnetic

forces), the B matrix relates the strain rate to the displacement rate (velocity)

of element and N represents the shape functions of the element. This


statement is valid quite generally for any stress-strain relationships. The

assumption of negligible body forces and negligible acceleration for each

element allow the equivalent nodal force equation to reduce to the following:

q' = fB'd(vol). (2-27)

For mechanical equilibrium where the body is not under rigid body motion,

net force is equal to zero.

An Hookean elastic element with the following constitutive relation

a = D(e- e,)+ o (2-28)

would be modeled by substituting the constitutive relation into the

equivalent nodal forces equation:

q = [B'DEd(vol) BDoDCd(vol) + BTd(vol). (2-29)
V' V V'

The strain e is related to the unknown displacement Aa through the B matrix

e = BAae (2-30)

and may also be substituted into the elastic equivalent nodal force equation:

q' = LBTDBd(vol) Aae JBTDeod(vol) + B'Taod(vol). (2-31)
v I V V,

Under mechanical equilibrium, the elastic equivalent nodal force equation

reduces to the following discretized form:




where A is the area (volume in 3D) of the element. The left-hand side

represents the stiffness matrix of the element. The right hand side represents

the initial stress and strain state of the element. From this relation, the

displacement Aa is solved for globally and the current stress and strain state

can be derived from the displacement Aa.

A viscoelastic body is handled in the same manner. The Maxwell

viscoelastic constitutive relation

-+-= E (2-33)
G r

has the following solution for the stress oa

a= 7 1-exp( -t + Coexp(F (2-34)

where r is the relaxation time constant and is the ratio of the viscosity r7 to

the elastic modulus G

=--. (2-35)

The Maxwellian viscoelastic constitutive relation can then be substituted into

the equivalent nodal force equation:

q' = [fB T 1-exp A Dd(vol) + f BT exp (vol). (2-36)

The strain rate e is related to the unknown change in velocity Av through the

B matrix

e = BAv' (2-37)

and may also be substituted into the viscoelastic equivalent nodal force


q' =[ B T 1-expt DBd(vol) Av' + fBT o exp (vol). (2-38)
Lv,. I v ).lJ _T v] I T

Under mechanical equilibrium, the viscoelastic equivalent nodal force

equation reduces to the following discretized form:

B {l 1-exp(1 DBA Ave = -BTa exp( -At (2-39)

where A is the area (volume in 3D) of the element. The left-hand side

represents the stiffness matrix of the element. The right-hand side represents

the initial stress and strain state of the element. From this relation, the

unknown change in velocity Av is solved for globally and the current stress

and strain rate can then be derived from the new velocity change Au.

The viscoelastic formulation reduces to the elastic formulation for

large relaxation time constant z. If z >>At, then

exp =l- A~ (2-40)

The viscoelastic formulation then becomes

SBT tDBA Ave = -B'o(1- At (2-41)
f T I I T )

and reduces to the elastic formulation by allowing AtAve=Aae.


To model two-dimensional problems, the plane strain formulation is

used. This formulation can be used for problems where the strain component

in the z-direction is zero or negligible [28]:

EZ = E = E=c =0. (2-42)

This can be approximated solving problems with infinitely long dimensions in

the z-direction. Therefore, the strain in the z-direction will approach zero. As

was mentioned before in the earlier sections, a problem arises using the

plane strain approximation while computing thermal mismatch and dopant

induced stress. These sources include a hydrostatic strain field Eo described as

the following:

D, = EC = E. (2-43)

Then plane strain presumption implies that stresses in the z-direction will

still occur even if there is no z-component of strain. These stresses occur due

to dopant and thermal expansion and are affected by the elastic constants. To

account for this using a plane strain approximation, the following expression

is used instead for an elastic relationship [28]:

o = E] =(l+v) e (2-44)
_C 0.

where v is the Poisson's ratio.


The FEM formulation uses three noded triangular elements called

'faces' for two-dimensional applications and four noded tetrahedral elements

called 'volumes' for three-dimensional applications. Linear shape functions

are used for interpolation of the strain solution within each element.

Reflecting boundary conditions are used to handle the strain solution

at the boundary of the simulation field. The normal component of the velocity

field or displacement rate is set to zero across this interface. Physically this

corresponds to a 'mirror reflected' symmetrical structure across the boundary.

2.3.1 Boundary Loading Method

The boundary loading method (BL) uses the strain solution calculated

during the oxidation to drive the elastic solution in the silicon. This technique

has been used in the SUPREM IV process simulator to calculate the elastic

strain in silicon due to oxidation [2]. Chidambarrao also used this technique

to investigate strains resulting from isolation trenches [45]. It is fully

decoupled and is performed in two sequential steps. This is demonstrated in

Figure 2-9. In the first step, the nonlinear viscoelastic oxide flow is computed

with the surface films modeled as viscoelastic materials and with silicon

acting as a rigid body [44]. This assumption allows for a more efficient

technique for solving the oxide growth, since only the surface films are

iterated over in the nonlinear stress dependent oxidation solution.

The stress tensor in the silicon dioxide elements along the silicon -

silicon dioxide interface are averaged to interface nodes. The averaged node

Step 1-Calculate the stress in the upper surface films to compute the stress-
dependent oxidation growth
Surface Film Forces
Oxide Growth__ -- I
Forces )k
Step 2-Calculate the stress induced in the substrate due to the forces
generated during oxidation and strain
generated from dopants. _

Strain -- dopants, defects

Figure 2-9 The boundary loading method is composed of two separate
steps for computing the strain in silicon.

stresses are then converted into boundary forces along the interface. The

boundary forces are simply the product of the stress tensor with the unit

surface normal of each node:

F, a,, aT zXr n.
F: = 1 (Y a n* (2-45)
z o' yz azJ _z

These forces are then input as boundary loading forces that drive the

calculation of the elastic substrate strain solution. Silicon is modeled as an


isotropic elastic material in the second step. Strain may be included that is

exerted from dopants and defects within the silicon for this solution.

2.3.2 Fully-Integrated Method

The other technique developed to investigate the strain in silicon is

similar to the finite element method implemented by Senez [46]. This fully-

integrated method couples the strain solution in the silicon along with the

surface films to compute the oxidation growth. This method allows for the

oxide stresses to relax by exerting forces on the silicon substrate. In the BL

method, oxide stresses can't relax due to the rigid body boundary condition

imposed on the silicon. Therefore, lower magnitudes of stresses are expected

in using this method. This method is extended to include strains from other

sources simultaneously (Figure 2-10).

This technique becomes very computationally intensive because the

silicon elements are also assembled into the nonlinear oxidation equation.

This usually results in an order of magnitude more elements included to

adequately define the substrate and reduce the effect of the reflecting

boundaries on the solution. Therefore, the FI method's Newton iteration

involves a much larger matrix than the BL method and results in much

slower performance overall.

As in the previously described method, the surface films are modeled

as viscoelastic bodies and the silicon substrate is modeled as an isotropic

elastic material. The boundary condition at the oxidant reaction interface is

6_- -_Surface Film Forces

Oxide Growth, -~
Forces \_

Strain -- dopants, defects

Figure 2-10 The fully integrated method couples the solutions of strain in
the upper films along with the substrate simultaneously.

similar to a polysilicon boundary condition: the silicon flows and is consumed.
Additionally, strain from dopants and intrinsic film stress may be exerted in
the silicon substrate and also is included as a stress source in this method.

2.4 Summary

One of the goals of this work are to develop a system where strain can
be computed from multiple sources independently and simultaneously. In
this chapter, a brief review of continuum mechanics was first provided to
refresh the reader with the necessary terminology. Next the individual strain
sources and the algorithms used to model the mechanical strain generated


from each source were described. Finally the boundary loading (BL) and

fully-integrated (FI) finite element were introduced and discussed for

modeling the strain due to oxidation volume dilation. The optimized BL

method is enforces a rigid body assumption for silicon. This allows it to be

useful for analyzing oxidation-induced strains in three-dimensions, as will be

discussed in the chapter three. However, this assumption may produce

influences in the oxide growth and stress solutions. Therefore, in chapter

three, the two methods will be analyzed for several applications to find if and

when there are influences. Also, in chapter three, applications using the

film-induced and boron-induced strain models are demonstrated and they

are compared and validated with experimental studies.


Due to the fully integrated method's computational intensiveness, the

boundary loading method is more commonly used. This brings into question

how well the solution from each method agrees. Comparisons are performed

between each method for computing the oxidation-induced stresses for

different processes in silicon fabrication. Next, solutions of nitride film-

induced stresses are compared with micro-Raman spectroscopy

investigations. Afterwards, boron-induced cantilever bending simulations are

performed and compared with previous experiments. Finally, a three-

dimensional example is demonstrated for simulating oxidation-induced stress

using the boundary loading method.

3.1 FEM Comparisons

To compare the two methods for computing oxidation-induced stress,

two different oxidation processes are analyzed. The first is a standard

LOCOS process. The second is a re-oxidation process over an oxide filled

shallow trench. These examples detail when the boundary loading method

solutions fail to match the solutions computed by the fully-integrated


Each of the following simulations is two-dimensional and use the plane

strain approximation with reflecting boundary conditions. This means that

the dimensions in the z-direction (which extends outward perpendicular to

the plane of the page) are infinitely long. The reflecting boundary conditions

indicate that the structure is periodic along the right and left boundaries of

the field.

Hydrostatic pressure contours are plotted for each process to compare

solutions computed by each method. This allows for a scalar interpretation of

the stress tensor's variation in the two-dimensional field. The hydrostatic

pressure definition used is simply the negative average of the normal

components or negative one-third of the trace of the stress tensor:

CXX 0 0
P = 0 Y Y (3-1)
0 0 CTZ_

For each wet ambient simulation, material properties calibrated at

1000C by Cea and Senez [44, 46] are utilized. These parameters are listed in

Table 3-1.

For the STI simulations where dry oxidations are also performed, the

oxide material properties were calibrated to Kao's dry oxidation cylinder

Table 3-1 Material coefficients at 1000C [44, 46]

Elastic Parameters Viscosity
material Elastic Poisson's Viscosity Vo
Modulus Ratio (Dyn/cm2)
oxide 6.6ell 0.17 2.0e14 450e-24
nitride 3.89e12 0.3 7.0e14 100e-24
silicon 1.7e12 0.3 -___

experiments [13], similarly to how the wet oxidation parameters were

calibrated [44].

In Kao's experiment, cylinders and holes of varying radii were oxidized

to grow approximately 3500A at various temperatures in a dry ambient. In

both the cylinders and holes the oxide growth was retarded when compared

to plane wafer oxide growth. It was concluded that this decrease in oxide

growth was due to the stress dependence on the diffusivity of oxidant

reactant, reaction rate, and viscosity of oxide. As the radii of the cylinders

and holes decrease, the amount of oxide growth also decreases. A standard

calibration procedure adopted by Rafferty [2] and Cea [44] is used to obtain

the oxide parameters that fit the curves for each temperature and radius

geometry. The results of the cylinder and hole simulations are plotted in

Figure 3-1. The dry oxide parameters obtained from the calibration

simulations are listed in Table 3-2 along with the wet oxide parameters

obtained by Cea [44].

Dry Oxidation Calibrations

Figure 3-1

0 1 2 3 4 5
1/R (Si) (um)
Results of dry oxidation calibrations for holes and cylinders
calibrated to Kao's measurements [131

The calibrated low stress viscosity of dry oxide is plotted over

temperature in Figure 3-2 along with the wet oxide viscosities calibrated by

Cea [44]. It can be seen that the dry oxide low stress viscosities are about an

order of magnitude less than the corresponding viscosity at each

temperature. This is expected since dry oxidations produce a denser oxide.

Also, the viscosities follow along the same slope in the Arrhenous plot with

similar activation energies.

Calibrated dry and wet [44] oxide parameters.

Oxide Viscosity Diffusivity and
Parameters Reaction Rate
T(C) Tlo Vo VD VK

(dyn/cm2) (A3) (A3) (A3)
1100 W 3.3e13 1100 75 10
1000 W 2.0e14 450
900 W 1.5e15 300
1100 D 3.5e14 250
1000 D 2.2e15 150
900 D 2.5e16 6.25 "

Oxide Viscosity versus Reciprocal Temperature

0.74 0.76 0.78 0.8
1000o r (1/K)

0.82 0.84 0.86

Figure 3-2 Calibrated dry oxide viscosity compared with the wet oxide
viscosity calibrated by Cea [44].


106 ....




Table 3-2

3.1.1 LOCOS

The LOCOS process is simulated to compare between two-dimensional

stress solutions computed by each method. Figure 3-3 exemplifies

simulations of a 10000C eighty-minute wet oxidation process at atmospheric

pressure. A 180nm thick nitride film with a two-micron linewidth is

patterned over a 20nm thick pad oxide. The grid spacing used for the

simulations is 0.05am at nitride edge, 0.15 m at the left and right side

boundaries, and 0.25pm at the bottom boundary.

Figure 3-3 displays the stress solution for each method. By analyzing

the average pressure contours displayed, it is observed that there is

qualitative agreement between the two methods for this particular LOCOS

process and geometry. Corresponding regions are in tension (negative

pressure) and in compression (positive pressure). However the relative areas

occupied by the corresponding pressure contours indicate that their

magnitudes disagree. This is expected since the fully integrated method

allows for some stress relaxation through a 'flexible' silicon base. The

boundary loading method does not allow oxide stress relaxation through the

silicon because the silicon is modeled as a rigid body. The presence of both

compressive and tensile regions in the silicon is due to its elastic behavior.

When a compressive force is applied to an elastic body (such as at the bird's

beak in silicon), a neighboring region in the material tends to deform in

tension due to the Poisson's ratio effect.

The solutions in the oxide agree more closely in magnitude than in the

silicon. It is for this reason that the bird's beak heights (BBH) are almost

identical as can be seen in Figure 3-4. This indicates that the rigid body

boundary condition used for computing the oxide growth, as modeled in [44],

does not affect the oxide stress solution enough to alter the oxide growth

significantly at longer nitride linewidths. However at shorter linewidths, it

becomes evident that the BB curvature underneath the nitride edge is

influenced by the strain computing method. In Figure 3-4, an overlay of

simulations using both strain computing methods is shown. The FI method

solution shows a little more oxidation under the nitride edge. This indicates

that the BL method computed a higher stress in this region that retarded the

oxide growth. A matrix of LOCOS simulations is performed to analyze the

differences of the stress and oxide growth solutions for shorter nitride

linewidth simulations using each FEM method.

Table 3-3 represents a matrix of simulations with varying nitride

linewidths. The first set with linewidth of 2gm is demonstrated in Figure 3-3.

Depending on the amount of elements in the field, the boundary loading

method is about two to three times faster than the fully integrated method.

Table 3-3 demonstrates that the oxide thickness and BBH is consistent

between the two methods and the variation in oxide shape is small until the

nitride linewidth decreases to 0.5pm. Below this linewidth, the oxide punch-

through effect begins to occur and the BBH varies more significantly

depending on the strain computing method used. The maximum tensile and

compressive stress magnitudes in the silicon are also noted for each

simulation and tend agree.

Table 3-3 Matrix of LOCOS simulations depicting how the two FEM
solutions correlate with decreasing nitride linewidth.

linewidth BBH Max Oxide Max Max CPU
(urn) (um) Growth (um) Tension Compression time (s)

2.0 FI 0.276 0.499 -1.5e9 4.9e9 1077.4

BL 0.268 0.499 -3e9 5e9 633.9

1.0 FI 0.233 0.453 -3.7e9 7.5e9 1638.5

BL 0.235 0.46 -3.4e9 5.1e9 552.5

0.7 FI 0.215 0.415 -6.5e9 6.9e9 895.6

BL 0.215 0.42 -8e9 5e9 343.0

0.5 FI 0.179 0.367 -5e8 8e9 821.6

BL 0.190 0.368 -5e8 8e9 300.2

Y Axis BL Method

0.00 +,09 0.00
+108 1 9 108 +109

1.o00 -108 1.00 +108 1

.................... 10%
2.00 1 8j Compression 2.00

-_1.0x10o Tension ___ 10io8 t Compression

_i O Tnin8 8 Tension
---l -1.0x10

3.00 __ 3.00

-1.00 0.00 X Axis 1.00 -1.00 0.00 X Axis 1.00

Figure 3-3 Comparisons for a LOCOS-induced stress computed by the boundary loading method (left) and the
fully-integrated method (right).

FI Method

Y Axis

Oxide Shape Comparisons

x in microns




-0.60 0.00 y in microns 0.60

Figure 3-4 Overlaying oxide shape solutions comparisons between the BL
method and FI method for an example LOCOS simulation with 2gm nitride

3.1.2 Post-STI Process Re-Oxidation

Strain Solution Comparisons. An oxidation process over a previously

processed oxide filled shallow trench isolation (STI) is simulated to compare

between two-dimensional stress solutions computed by each method. Figure

3-5 exemplifies a simulation of a 10000C three-minute wet oxidation at

atmospheric pressure. The trench with 0.5pm depth and 0.5gm width is filled

with deposited oxide. The nitride mask in this simulation has been previously

stripped as in a CMP process. The grid spacing used for the simulation is

0.02pm at the trench edges, 0.20am at the left and right side boundaries, and

0.10pm at the bottom boundary. The initial field stress solution is set to zero

to compare the stress induced by the oxidation alone.

Figure 3-5 is an example of a 'wedge driving effect' that is an oxidation

performed over a fabricated trench that is not masked. A vertical bird's beak

forms at the upper corners of the trench due to reactant diffusion into the

trench. This produces a very high compressive stress in the silicon substrate

along the side of the trench as oxide is grown on the trench sidewall.

Therefore this is analogous to a wedge being driven in at the trench sidewall.

By observing the average pressure contours displayed, it is observed

that there is qualitative agreement between the two methods for this process.

The tensile and compressive regions again correspond between the two

method's solutions. The relative areas occupied by the corresponding pressure

contours indicate that their intensity variations disagree. The BL method

computes that the entire field is in compression with the exception of the

silicon at the bottom trench corner. Therefore, the two methods agree

qualitatively in the strain solution computed for this process. However, for

STI re-oxidation processes the strain solution is only part of the concern.

Y Axis BL Method

Y Axis FI Method

0.00 '00 1 -. 0.00

+109 +109 +109 +109

0.50 0.50
0. -10 109

y +108 -o108 +108

1.0 1.0

Compression 10
10 7
1.5 109 1.5 Tension -1.0xlO
1-o8 Compression -1.oxo
S -1 .x i o Tension
_1 Ox10-'

2.0 21 2.0

0.00 0.50 XAxis 1.00 0.00 0.50 X Axis 1.00

Figure 3-5 Comparisons for post STI process-induced stress computed by the boundary loading method (left) and
the fully-integrated method (right).

Corner Rounding Comparisons. One of the most critical design aspects

of the STI is the top corner shape. It is desirable for the top corners to be

rounded for several reasons. The most important being that sharp corners

lead to high electric fields and result in undesirable electrical characteristics

[47, 48]. Another reason for avoiding sharp corners is they lead to higher

mechanical stress magnitudes that generate dislocations also deleterious to

device operation [7].

Re-oxidation after fabrication of the STI has been introduced as a

method to round the upper trench corners [48]. The increased radius of

curvature at the upper corners alleviates these problems. Selecting the most

appropriate oxidation process for this technique can be simplified through

oxidation process simulation [49].

Viscoelastic flow models are used for the grown and deposited films

during the re-oxidation process [44]. Stress-dependent models for oxide

viscosity, reaction rate, and diffusivity are utilized for simulating the re-

oxidation anneals. Two-dimensional plane strain simulations are performed

for a re-oxidation process to study the evolution of the radius of curvature of

the upper STI corners. The re-oxidation occurs after the planarization step of

the STI process. Therefore the oxide fills the trench up to the surface of the

nitride film (Figure 3-6). The STI dimensions have a depth of 300nm and a

width of 400nm. A series of re-oxidation simulations is performed varying the

nitride thickness from 70nm to 180nm over a 30nm pad oxide. The re-

x in microns FLOOPS

Nitnde Lift-ng__ ~,' -._
0.00 i Bird' s Bak Le gt
NI Raius of Curmtwe




-0.20 0.00 y in microns 0.40

Figure 3-6 10000C 300 minute dry re-oxidation simulation.

oxidation processes simulated and analyzed are a 300 minute 10000C dry

anneal, 100 minute 11000C dry anneal, and a 30 minute 9500C wet anneal.

The oxidation processes result in about 150nm of oxide grown on a blank

wafer. Each process is examined for each of the nitride film thicknesses.

The radius of curvature relationship to nitride thickness is consistent

with the experimental results of Chang [48] (Figure 3-7). The amount of

nitride lifting also varies linearly with the nitride thickness. The bird's beak

lengths are fairly constant over this range so there is little trade-off against

the thinner nitride process.

70o (a) -22 -----.......--.---.------- (b) --

55 18
50 4 16
45 114
30 50 70 90 110 130 150 170 190
50 70 90 110 130 150 170 190
50 70 iD NiNde Thickne. (m) 0 1 1 Nitride Thickess (nm)

Figure 3-7 (a) The radius of curvature of the upper STI corner and (b)
nitride lifting at the bird's beak increases with decreasing nitride mask
thickness for the 1100C dry oxidation.

Figure 3-8 depicts how the various process simulations examined

resulted for a given nitride thickness using each oxidation strain computation

method. It can be seen that for the higher temperature 11000C dry oxidation

resulted in the greatest amount of corner rounding and stress relaxation.

x n microns FLOOPS 98.1




0.10 950C wet
950C Wet
0.20 1100C Dry


0.00 in microns
Figure 3-8 Overlay of the re-oxidation profiles for each anneal process
with 120nm nitride film thickness.

Also, due to less stress relaxation in the lower temperature 9500C wet

oxidation, less corner rounding is resulted by the BL method and an

overhang results by the FI method. The overhang also agrees with re-

oxidation experiments [48]. However its shape is different since the

simulated shape may be locally grid dependent. The higher stresses reduce

the oxidation rate at the upper silicon corner just below the nitride edge and

result in excessive compressive stress magnitudes in the neighboring silicon


From this application, it is found that for the low temperature wet re-

oxidation simulations, the FI method results in an overhang at the trench

corner due to the severe oxide growth retardation. However for the same

simulation, the BL method returns a rounded corner. Since in experiment

[48], this process does result in a slight overhang, it is believed that including

the silicon deformation while calculating the oxide growth plays an important

role in modeling this effect at the trench corner.

3.2 Raman Spectroscopy Measurements and Comparisons

In the previous section qualitative comparisons were performed for

different oxidation processes between the two FEM methods implemented. In

this section, efforts are made to compare the nitride film edge-induced stress

with measured strain data. Currently, the best method for measuring local

process-induced mechanical strains in the substrate at sub-micron

dimensions is micro-Raman spectroscopy [26].

Raman spectroscopy is most commonly utilized for crystallographic

and chemical composition analysis. Mechanical strain or stress affects the

frequencies of the Raman modes of the crystal. The Raman spectrum of the

silicon crystal is sensitive to strain in the lattice and therefore has been used

to measure the stress state induced by microelectronics fabrication processes

[50]. Micro-Raman stress measurements involve focusing laser light onto a

sample. The scattered light of the sample is then collected and analyzed. The

light scattered from a (100) unstressed silicon sample has a Raman peak

frequency of about 520 Rcm-1. Stressed samples will influence a shift in the

Raman peak. The amount and direction of the shift corresponds to the

magnitude and sign compressivee or tensile) of the stress in the laser spot

region of the sample. A negative shift corresponds to a tensile stress and a

positive shift indicates a compressive stress.

3.2.1 Raman Simulation Method

A method for determining the expected frequency shift in Raman

modes due to process-induced strain was introduced by Jones [51] and De

Wolf [52]. This method is implemented here to validate the stress solution

from the nitride film edge-induced stress model with micro-Raman

measurement data. The method involves first evaluating the shift in Raman

modes for a given strain tensor solution. For an elastic material, the strain

tensor is related to the stress tensor through the compliance tensor sijkl:

E- = Sik kk. (3-2)

The compliance tensor sijkl is simply the inverse of the stiffness tensor cijkl

introduced in chapter two. Once the strain solution is known, the three

Raman mode shifts (AoI, Ao)2, Ao3) are computed by solving the following

eigenvalue problem [52]:

PE +q' E22 +q33 (P-q)e12 2rE,3
(P- q)eL2 E22+ q33+ q' El- 2rE23=0 (3-3)
2rE,3 2rE23 p33 +q(eL +E22)-

where Ao = 22oo.

Next, the Raman mode shifts are convoluted to account for the

penetration of the laser light into the silicon crystal and the width of the

laser spot. An exponential decay of light intensity and gaussian intensity

distribution over the spotwidth is assumed. These corrections average the

Raman mode shifts over a volume of the crystal corresponding to each spot


3.2.2 Nitride Film Edge-Induced Stress

A simulation of the stress induced by a nitride-poly stack is used to

compare with a Raman measurement experiment. A 0.24gim thick nitride

film with intrinsic stress of 1.2e10 dyn/cm2 is deposited over a 0.05pm poly

film with intrinsic stress of-3e9 dyn/cm2. The stack rests over a 0.01 m pad

oxide and is then patterned for a linewidth of five microns. The stress

induced by the film edges is then solved for and illustrated in Figure 3-9 by

stress contours of hydrostatic pressure. This represents the stress state of a

Poly-Buffered LOCOS structure before oxidation. Finally using the method

described in the previous section the three Raman mode shifts are solved for

and are displayed in Figure 3-10.

This simulation is then compared to Raman measurement data

introduced by De Wolf [52] for the same process. Since the magnitudes of the

x in microns

.oo /

1.00 8 I x -1.0x108
-1.Ox10 i


3.00 8



0.0 y in microns 10.0

Figure 3-9 Simulation of the stress induced by the intrinsic stress of a
nitride poly stack for a linewidth of 5gm. Units are in Dyn/cm2.





Position (cm)
Figure 3-10 Frequency shift simulation for the three Raman modes due to
nitride-poly stack edge induced strain compared with Raman measurement

Raman mode shifts are small, only one shift can be recovered in experiment,

which is an averaged function of the three Raman mode shifts. This

simulated Raman shift is plotted along with the measured Raman shift [52]

in Figure 3-10. It is noticeable that frequency shift predicted correlates well

with the De Wolfs measurements. The resulting compressive strain induced

shift is consistently lower than the measured Raman shifts. This indicates

that the compressive strain magnitudes under the nitride film may be under-

predicted. Other reasons for the discrepancy may be film-related effects on

the incident signal reflection, or convolution averaging errors due to light

intensity decay and spot width averaging.

3.3 Boron-Doped Cantilever Bending Comparisons

Next, an application using the boron-induced strain model is presented

in order to validate its solutions. The strain-induced bending of boron doped

cantilevers are simulated and compared with previous experimental


3.3.1 Silicon Bulk Micro-Machining

Silicon bulk micro-machining is important for fabricating silicon-based

sensors and transducers. Silicon sensors are often composed of thin

membranes, bridges, cantilevers, and beams. These structures can be

fabricated by various bulk micro-machining methods. Anisotropic wet

chemical etching is often used to develop sensor structures due to its

simplicity and convenience as well as providing very accurate dimensional

control [18].

Boron etch stops are often used as a method for controlling etch depth

in silicon substrates. Thin silicon film structures can be fabricated by

thermally diffusing or implanting boron on one surface of the silicon wafer

and then by etching through a mask window on the other side of the wafer.

For wet chemical etchants such as KOH, the etch rate decreases significantly

as the etch front approaches boron concentrations greater than 7x1019 cm-3. It

is believed that the strong B-Si bond tends to bind the crystal more

stringently, therefore requiring more energy to release the silicon atom [53].

It is then possible to design thin silicon film structures with the

desired thickness by controlling the diffusion of the boron dopant profile so

that the etch stop will occur at depth where the boron concentration

approaches ~7x1019 cm-3. However due to these levels of boron concentration,

high levels of residual stress are generated. Since micro-machined thin

membranes are critical components of silicon sensors and transducers,

residual stresses in these structures may lead to mechanical failure of the

device and/or deteriorate its performance. Figure 3-11 illustrates the bulk

micro-machining process performed by Yang [39].

# Si02



Thermally Diffused Boron

Backside Etch


Figure 3-11 Bulk micro-machining process flow for fabrication of the

First boron is introduced by thermal diffusion and is masked by an

oxide layer on the backside of the wafer. Next the front side is masked and

the backside is chemically etched. The etch slows as the surface approaches

the high concentration of boron (~6x1019 cm-3). Next the front side mask is

stripped off and the cantilever is developed. In Yang's experiment, a series of

front side etches then performed to study the bending behavior of cantilevers

of different thicknesses doped with the same boron profile. The relative depth

of the peak of the boron profile then shifts for each cantilever etch.

The high concentrations of boron necessary to produce the etch stop

behavior results in residual tensile stress with magnitudes approaching and

exceeding levels of 1x109 dyn/cm2. To relieve these high levels of stress the

silicon crystal yields and may generate dislocations that may be deleterious

to device and sensor performance [8]. This is one of the main reasons for

studying residual stress and its origins.

The residual stress resulting is dependent on the gradient and

maximum magnitude of the boron dopant profile as well as the thickness of

the cantilever resulting after the backside etch. Since the boron dopant

profile is not uniform the stress distribution varies with depth causing the

cantilever to bend in order to relieve the resulting residual stress. This is

evident in previous studies analyzing positive and negative bending of boron-

doped cantilevers under varying diffusion conditions [39, 40].

3.3.2 Cantilever Bending Simulations

Simulations are performed using the finite element models previously

described for the cantilever process described in Figure 3-11 [54]. First boron

is introduced by thermal diffusion. The resulting profile has a peak

concentration of 8x1019 cm-3 (Figure 3-12). A backside wet chemical etch is

then performed. A boron concentration of 6x109cm-3 is chosen as the stop for

this etch. The resulting cantilever thickness is about 1.4gm.

Next a series of 2d plane strain and 3d elastic deformation simulations

is performed for the resulting boron-doped cantilever structures. The 2d

cantilever has a length of 50pm. The grid spacing in the x-direction is 0.05gm

and in the y-direction is 0.5gm. Since the x-direction spacing is limited by the

necessary resolution to represent the boron profile, it becomes a challenge to

preserve element quality as the cantilever length is increased while





nCoi.c(-o) FLOOPS____
20? ____


19 -_________________

19? _____________ ____

-200 Dep1 (micrni) -1 00

Figure 3-12 Resulting cantilever boron profile with a thickness of 1.4pm.

maintaining the number of elements constant. The element aspect ratio

problem is magnified further in three dimensions. The 3d cantilever has

dimensions of 1am width and 10m length. The 3d cantilever grid spacing is

0.3pm in the z-direction (width), 0.5am in the y-direction (length), and O.1pm

in the x-direction (thickness).

The following mechanical material properties were utilized for all

simulations performed: Young's modulus (Ep+si = 1.22x1012 dyn/cm2) [55],

Poisson's ratio (v = 0.3).

A series of varying front side etch simulations were performed to

analyze how the shift in boron profile and decrease in thickness of the

cantilever affected the deflection solution. These series of etch simulations

model the experiment performed by Yang [39].

The results of the 2d plane strain simulations are displayed in Figure

3-13a. An example of the deflection simulation for the 0.62pm thickness beam

is shown in Figure 3-14. The structure is regridded after the nodal

displacements are solved for. Figure 3-13 displays the nodal displacements

along the top surface of each cantilever structure simulated. The differences

between each structure thickness relates to the amount etched off the top

surface. Notice that all the cantilever beams except the thickest deflected in

the negative x- direction (upward in Figure 3-14). Generally as the

cantilevers were etched thinner, the amount of deflection increased.

A 3d simulation for the 1.37pm thickness cantilever is demonstrated in

Figure 3-15. It is more difficult to examine the deflection visually in the 3d

simulations due to the shorter length of the cantilevers. The deflection

solution plot for the 3d simulations is displayed in Figure 3-13b. The same

general trend also results.

(a) 20 Simulation
50 um Length Deflections
80E-06 T

..: ..

.- ,b- .. -"
4 -rF -0

., 3 2-

1-t .... .......... .............. ............. ... ...... ....
0 000i 00o0 0003 0004
Y-Position along the Cantilever (cm)



I OE-07
- 000.00

-2 0E-07
0 -30E-07

3D Simulation
10 um Length Deflections



--" -a'-
-- 1I7 T ,
--1 22

0 0.0002 00004 00006 0.0000 0001
Y-Position along Cantilever (cm)

Figure 3-13 (a) Two- and (b) Three-dimensional simulation results for
deflection curves for cantilevers of varying thickness.

The simulation lengths of the cantilever beams are much shorter than

those fabricated in various experiments [39, 40]. Typically cantilever beams

are fabricated with lengths up to 1mm in order to have an accurate

measurement of the deflection. The element quality problem limits the length

of the cantilevers simulated before a significant error is resulted in the elastic


x in microns FLOOPS
-1.40 -
12'' ... '-" 'i2 i'K '" .. :;'7

Figure 3-14 Two-dimensional plane strain simulation of the boron strain
induced bending of a cantilever with thickness 0.62pm.

To compare the simulations performed with experiments in the

literature, a parabolic curve fit was used to extrapolate the expected

deflections for longer cantilever beams. This is possible because beams

processed the same with varying lengths would all have the same bending

.. ._- ... . .. _
0. 000393


0. 000000000

Figure 3-15 Three-dimensional simulation of the bending of 1.37um
thickness cantilever. The maximum deflection is at the bottom corner of the
cantilever beam.
cantilever beam.


moment [40]. Therefore, the cantilever deflection simulation results are then

extrapolated to amounts corresponding to 1250km to compare with

experiments performed by Yang [39]. These results are presented in the

histogram shown in Figure 3-16.

Several points can be deduced from the results obtained. First is the 2d

and 3d simulations resulted in roughly the same amount of deflection. This

confirms that the plane strain approximation does not affect the result of the

simulation and that cantilever width is not a factor in the stress solution.

Second, both the simulations agree with the published experiment in the

direction of the deflection for each beam thickness. However the measured

Cantilever Deflection Comparisons
200 ........................ ............................... ..... .. ..................................................................... .............
M Measured [6]
e 100-- Extapolated 2D Simulation
3 OExtapolated 3D Simulation
c -


1 -200


N 0 0 I
~400 ------------------

Various Thicknesses

Figure 3-16 Results of both sets of simulations compared with Yang
experiment [391.

deflection results are consistently about two to three times the magnitude of

the simulations. Also the measured quantities have a relative maximum

negative deflection for the 0.92pm thick beams, while the simulations showed

relatively constant deflections for cantilevers of less thickness. It is believed

that these differences are due to errors in the boron profile and a variation in

the lattice strain coefficient.

Cantilever self-loading (due to its weight) has also been suggested as a

possible source of extra strain. However the maximum bending moment due

to a constant distributed load on a cantilever supported on one end is

described as the following relationship [56]:

max = (3-4)
max 2

where L is the length of the cantilever and W is the uniformly distributed

load per unit length. For a body force such as gravity, W is represented as the


W = Apg (3-5)

where A is the cross-sectional area of the cantilever, p is the density of silicon

(2.33 gm/cm3 for 5el9cm-3 boron-doped silicon [41]) and g is gravitational

acceleration constant. For a cantilever with rectangular cross section and

thickness t, the longitudinal stress a is proportional to M by the following



a -A. (3-6)

The following gravity-induced stress relationship then results by substituting

(3-4) and (3-5) into (3-6):

a -3 (3-7)

It can then be shown that for a cantilever of 1000m length and lgm

thickness, the maximum stress produced is about 7e5 dyn/cm2, which is at

least three orders of magnitude less than the stress magnitudes induced by a

concentration of boron at 5el9cm-3.

The most probable cause for the difference between the simulations

and experiment is differences in the boron profile and lattice strain

parameter. Since the strain is computed directly from the local concentration

of boron and the lattice strain parameter, a shift in these results in a linear

shift in the strain and therefore the stress. Therefore scaling the lattice

parameter by a factor of three results in much closer agreement with the

deflection magnitudes in Yang's experiment. This also results if the boron

concentration is scaled by the same factor. It is therefore believed that the

differences in both factors result in the observed disagreement in deflection

magnitudes. This is highly probable since the accuracy in the concentration

of boron reduces as the depth increases and the lattice parameter may be

different for the thermally diffused case.

Another method to estimate the lattice contraction parameter, is to

subtract the covalent radius of boron (0.88A) from that of silicon (1.17A) and

then normalizing the difference to the silicon atomic radius. This results in

an estimate of the normal radius contraction when a boron atom substitutes

for a silicon atom:

R R 1.17-0.88
= R 1.17 -0.88 = 0.248. (3-8)
Rs, 1.17

This lattice contraction parameter results very close to that extracted from

Horn's measurements [41] (0.014A/5.43A)*(100).

3.4 3D Boundary Loading Method Example (LOCOS)

The computational advantage of the boundary loading technique

allows for less CPU intensive three-dimensional strain simulations. For two

dimensional oxidation simulations, the BL method averaged 2-5 times faster

depending on the number of elements in the silicon. This time savings is even

much greater for three-dimensional simulations since the number of

elements for a typical 3D simulation is an order of magnitude greater than in


Figure 3-17 illustrates a three-dimensional simulation of the stress

induced in a LOCOS process computed by the boundary loading method. The

process underwent a 10000C wet oxidation for 15 minutes at atmospheric

pressure. The pad oxide thickness is 10nm and the dimensions of the nitride

Figure 3-17 Three-dimensional simulation of LOCOS-induced stress in the

film are 1.0pm x 1.5pm x 0.15pm. The iso-stress contour lines demonstrate

how the stress levels are dependent on the geometry of the nitride film. In

the 1.5pim length dimension the stress level is greater than in the 1.0pm

dimension. This is evident by comparing the area occupied by corresponding

iso-contour lines.

Figure 3-18 and Figure 3-19 are top views of the same simulation. The

stress tensor component in the x-direction (directed normal to the surface) is

Figure 3-18 YZ plane view of the tensile oxx component of stress induced
during LOCOS oxidation.

plotted as contour lines. Figure 3-18 demonstrates that the maximum tensile

stress is at the corner of the nitride film. Figure 3-19 demonstrates that the

maximum compressive stress is located at the edges of the nitride film and is

greatest along the longer 1.5pm side of the nitride film.

3.5 Summary

The purpose of this chapter was to demonstrate through application

examples the capabilities of the strain models described in chapter two.

Another intention was to validate the results of various simulations with

experiments that have been reported. In the first section of this chapter,

solutions of oxidation-induced stress computed by the BL method and FI

Figure 3-19 YZ plane view of the compressive on component of stress
induced during LOCOS oxidation.

method were analyzed The LOCOS simulation example applications showed

that at higher nitride linewidths, the two method's solutions agreed.

However, at short linewidths, the two method's solutions began to disagree.

This lead to believe that the stress relaxation in the silicon becomes a critical

factor in solving for the strain as dimensions are scaled. Also, by analyzing

the re-oxidation simulations it was found that the simulations performed

with the FI method computed an overhang in the low temperature wet

process at the trench corner. The BL method simulation did not. Since in

experiment it was found that an overhang is resulted, it is believed that the

stain solution computed using the FI method was more accurate in retarding

the growth for that process.

Next, an application solving for the film edge-induced stress of a

nitride poly stack was demonstrated. The results of these simulations were

then used to predict the Raman signal measured over the nitride-poly stack

structure using a method outlined by De Wolf [52]. The expected frequency

shift in Raman modes then agreed both qualitatively and quantitatively with

Raman measurement experiments for the same structure.

Next, the boron-induced strain model was used to simulate the

bending of boron-doped cantilevers due to their boron-induced residual stress.

The simulations were then compared with studies of this effect and produced

qualitatively favorable results. The simulations agreed with experiment in

comparative bending between the different thickness cantilevers. The

difference in deflection magnitudes can be attributed to the lattice strain

parameter used in the model. By simply scaling this parameter, simulation

deflection magnitudes would coincide better with the measurements. One

possible explanation may be that the lattice strain parameter used may not

apply since it was measured in bulk silicon and may have a cantilever

thickness dependence.

Finally, three-dimensional simulations of oxidation induced stress

were demonstrated for a LOCOS process using the BL method. The BL

method's efficiency allows for three dimensional strain simulations to be

performed rapidly.


In the previous chapter, several applications were simulated using the

strain models described in chapter two. Those applications analyzed

individual sources of mechanical stress and how the modeled strain fields

compared with various experimental studies. In this chapter, the strain field

due to the shallow trench isolation (STI) process sequence is investigated.

STI has become an essential isolation scheme as CMOS technologies

are scaled down below the 0.25pm generation. However the basic STI process

sequence involves several sources of mechanically straining the enclosed

silicon region. Such sources include the thermal mismatch strain between the

isolation dielectric and the silicon substrate, intrinsic stress of the nitride

mask, and volume expansion-induced stress of the sidewall oxidation (Figure

4-1). It is therefore important to study each source and their combined strain

influences. As pitch lengths continue to scale, the distance between transistor

active areas and the STI sidewall also narrows. Therefore they continue to

approach the more highly stressed regions at the sidewall interface. This in

turn influences the device characteristics such as leakage through higher

dislocation densities and band gap deformation. Therefore for next

Figure 4-1 Major sources of stress in STI fabrication include the following:
(a) Nitride film edge-induced, (b) Insulator thermal expansion mismatch-
induced, and (c) Sidewall oxidation volume expansion-induced stress. Arrows
indicate force vectors.

generation process design, methods are needed to measure the stress

quantitatively. Up until now, micro-Raman spectroscopy has been the most

suitable method for investigating the stress [26]. But as technologies

continue to scale towards the 0.lpm generation, this may even surpass micro-

Raman's spatial resolution limits. In this chapter, scanning Kelvin probe

force microscopy (SKPM) is investigated as a new method for measuring

stress-induced influences in the silicon work function AO and for inferring the

stress magnitudes through coupling with mechanical strain simulation.

4.1 Scanning Kelvin Probe Force Microscopy

Scanning Kelvin probe force microscopy (SKPM) has previously been

researched as method for profiling 2D dopant concentration [57-62]. This has

been achieved by detecting the electrostatic potential difference (EPD)

AFM Cantilever Probe
+z \V

v V Rastered laterally


Figure 4-2 Cross-sectional sketch of SKPM measurement.

between the SKP tip and sample (Figure 4-2) through their electrostatic and

van der Waals force interaction.

The method involves rastering a cantilever probe over a cross-sectional

surface of the sample. At each step, both a surface topography measurement

and a surface work function measurement is performed. For measuring

surface topography the technique simply works as an atomic force microscope

(AFM). Under the AFM mode of operation, the mean distance between

cantilever probe and sample is kept constant over each step position. The

cantilever probe then is oscillated mechanically at a particular frequency

[58]. As the cantilever-sample mean distance changes, the van der Waals

force also changes causing a change in the oscillation amplitude. The

cantilever position is then corrected in the z-direction to return the cantilever

to its original oscillation amplitude. The corrected cantilever position is then

recorded as the cantilever is rastered across the sample and therefore maps

the surface topography as demonstrated in Figure 4-2.

Laboraory STI Array Sample
Laboratory Cross Section
Frame of Reference

Z [iTO]

Y [110]

x [00T]

Figure 4-3 Sketch of Kelvin probe scanning laterally (y-direction) and in
depth (x-direction). Variation in work function (EPD) is measured by
monitoring the electrostatic force contribution to the cantilever deflection.

Surface work function measurement is performed under the SKPM

mode of operation. First, an AC bias is also applied to the cantilever along

with the mechanical oscillation. An electrostatic force then results between

the cantilever tip and sample. A DC flat-band voltage bias is then applied to

the cantilever probe to null the electrostatic force. This DC bias is

proportional to EPD of the cantilever-sample system. The DC flat-band

voltage is then recorded at each raster step and provides a two-dimensional

map of the work function variation of the sample (Figure 4-3).

4.2 Work Function Influence

The EPD scanned by the SKPM method is a measure of the difference

in work functions of the cantilever tip and sample. This can best be described

with the energy band diagram of the system (Figure 4-4). It is similar to that

of a metal-insulator-semiconductor (MIS) device band diagram, with the

airspace between the cantilever and sample representing the insulator. In

this application, the cantilever tip is constituted of silicon with a layer of

titanium silicide grown on the surface for improved signal detection. When

the tip and sample are electrically connected, their Fermi levels become

aligned. The existing electrostatic force between the tip and the sample

causes the silicon bands to bend at the surface. Additional band bending is

also due to the presence of surface states and charges on the semiconductor


In SKPM mode, a DC flat-band bias is applied to null the electrostatic

force (Figure 4-5). This applied bias represents the EPD at each rastered

Ti-Silicide Tip vacuum Silicon

Figure 4-4 Energy band diagram of the cantilever-sample system at zero
applied DC bias.

point along the cross-sectional surface of the sample. The measured VEPD can

then be expressed as the sum of the flat-band correcting voltage VFB and

surface potential terms [61]:

VEPD = Vf +V, +- sE s,].


The flat-band voltage is the difference between the work functions of the tip

and sample:


The work function of the sample Osi is the energy difference between the

vacuum level and the Fermi level:

Ti-Silicide Tip

vacuum Silicon


A-X ---

Figure 4-5 Energy band diagram after VEPD is applied and the electrostatic
force is nulled.

V = Ti-Si Os,.