UFDC Home  Search all Groups  UF Institutional Repository  UF Institutional Repository  UF Theses & Dissertations  Vendor Digitized Files   Help 
Material Information
Subjects
Notes
Record Information

Full Text 
MODELING OF MECHANICAL STRESS IN SILICON ISOLATION TECHNOLOGY AND ITS INFLUENCE ON DEVICE CHARACTERISTICS By HERNAN A. RUEDA A DISSERTATION PRESENTED TO THE GRADUATE SCHOOL OF THE UNIVERSITY OF FLORIDA IN PARTIAL FULFILLMENT OF THE REQUIREMENTS FOR THE DEGREE OF DOCTOR OF PHILOSOPHY UNIVERSITY OF FLORIDA 1999 Copyright 1999 By Hernan Rueda To my family ACKNOWLEDGMENTS 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 neverending support, love and encouragement throughout my whole life. TABLE OF CONTENTS page ACKNOW LEDGM ENTS................................................................................... iv A B STR A C T ......................................................................................................viii CHAPTERS 1 INTRODU CTION ......................................................................................... 1 1.1 M motivation ............................................................................................ 1 1.2 StressInduced Effects in Silicon Fabrication................................... 3 1.2.1 Oxidation Influences ................................................................ 4 1.2.2 Diffusion Influences ............................................................... 6 1.3 StressInduced 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 PROCESSINDUCED MECHANICAL STRAIN MODELS..................... 16 2.1 Continuum Mechanics ............................................................. 16 2.1.1 The Stress Tensor......................................................... 17 2.1.2 The Strain Tensor............................................................ 20 2.1.3 StressStrain 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 FullyIntegrated 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 PostSTI Process ReOxidation............................................... 59 3.2 Raman Spectroscopy Measurements and Comparisons................. 65 3.2.1 Raman Simulation Method.................................... ........... 66 3.2.2 Nitride Film EdgeInduced Stress.......................................... 67 3.3 BoronDoped Cantilever Bending Comparisons ............................... 70 3.3.1 Silicon Bulk MicroMachining ............................................. 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 StressInducing Apparatus Design..................................... 109 5.1.2 StressWafer 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 SilicidationInduced Stress ................................................. 131 6.2.4 ThreeDimensional 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 MODELING OF MECHANICAL STRESS IN SILICON ISOLATION TECHNOLOGY AND ITS INFLUENCE ON DEVICE CHARACTERISTICS By 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 submicron 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 ObjectOriented Process Simulator (FLOOPS). 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 submicron 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 straininduced 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. CHAPTER 1 INTRODUCTION 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 processinduced 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 everdecreasing 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. Offstate currents are critical design characteristics in logic and memory circuits, limiting the switching speeds, causing increasing power consumption and limiting reliability. Shallow Trench Isolation (STI) is steadily becoming the predominant isolation technology as minimum feature sizes decrease below the minimum attainable pitch lengths of LOCOSbased technologies. The general STI process flow includes a nitridepattered reactive ion etch, sacrificial sidewall oxidation, oxide deposition and finally a chemical mechanical polish. However, stressinduced 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 straingenerating 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 StressInduced 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 filminduced 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 submicron 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 SiSiO2 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 oxidation. 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 submicron 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 StressInduced 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, 1922]. 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 LOCOSinduced 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 (11) 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 relationship: E = pi. (12) 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 (13) P4 0 AP4 P5 0 AP5 P6 0 Ap6 The piezoresistive coefficients relate the stressinduced 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 ( * (14) 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 11 [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 increases. Table 11: Piezoresistive coefficients for silicon [17, 18]. Material p 711 7X12 744 (4cm) (1012 cm2/dyne) pSi 7.8 +6.6 1.1 +138.1 nSi 11.7 102.2 +53.4 13.6 The accepted explanation for this phenomenon is the manyvalley 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 pn 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: V V IF = I exp() + IRO exp( ) (15) VT 2VT where the saturation current is I,= qA po coth + Dn coth (16) L L" L. and the recombination current is qAniW IRo qAn (17) 270 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 2 no = (18) PpO The intrinsic carrier concentration is exponentially dependent to the stress dependent band gap Eg E ni =KT 3/2 exp( g). (19) 2kT 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 pn 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 processinduced 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 bandgap 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 submicron strain measurement. MicroRaman 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 microRaman'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 pnjunction 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 described. 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 edgeinduced stress solutions are validated with published microRaman measurements. Finally, threedimensional 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 pnjunction device operation is then described in chapter five. This wafer bending experiment addresses uniaxial stress influences on the forward current. This allows for 15 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. CHAPTER 2 PROCESSINDUCED MECHANICAL STRAIN MODELS 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 stressstrain 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 21) 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 F, F2 +z +Y V +x S Figure 21 Continuous body with external forces applied. 18 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 22). 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 22) of infinitesimal Txx +z Y" t 0 I / +X Figure 22 An infinitesimal cubic element located within a continuous body with stress tensor components shown. 19 dimensions located at point P in Figure 21. 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 xdirection is the following: T, = a, + f+ (21) Let the surface AS be the plane of the cube normal to the xdirection. The stress vector Tx is defined as the ratio of force acting on that surface area AS: AF dFP T = lim (22) 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 (23) In the definition above, ii are the normal stress components acting on the faces perpendicular to idirection and zj are the shear stress components oriented in the jdirection on the face with normal in the idirection. At mechanical equilibrium, it can be shown that the stress tensor is symmetrical [29], Tij = Tj . (24) 20 A column vector of six independent components can then describe the state of stress at a point: UT =[UXX ayy U z, Tu Ty ]. (25) 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 23 shows a onedimensional example of the deformation in a (a) (b) Unstretched Stretched Ax t Ax + Au v Figure 23 Onedimensional deformation of a spring: (a) original length (Ax), (b) deformed length (Ax+Au). 21 spring. The Lagrangian strain and Eulerian strain then, respectively, become the following over the length of the spring: increase in length Au E (26) S original length Ax increase in length Au Ei (27) 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 (28) Axo 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 = (29) e", C =i +) (29) dw 1(du dw cE E E i+ U U x 2 dx 22 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 EZ xy 'xzl Eyy Cyz SZ zz (210) A threedimensional example of a body undergoing deformation due to an externally applied force is illustrated in Figure 24. In the example a compressive force is applied to the infinitesimal cube in the xdirection. A negative displacement compressivee strain) results in the xdirection and positive displacements (tensile strains) result in the y and z directions. The dv /Zf 7\ /* +z + +x Before deformation Compressive Force in x direction S du + du/dx i du/dx < 0 dw/dz > 0 dw + dw/dz dv/dv > 0 < dv + dv/dy  After deformation Figure 24 Strain reference example displaying a compressive stress in the xdirection that generates normal strains in the x, y, and zdirections. 23 relationship of how the strain in each dimension results from the applied force will be discussed next. 2.1.3 StressStrain 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 pressure. 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 conditions. For a Hookean elastic solid, the stress tensor is linearly proportional to the strain tensor over a specific range of deformation: (211) T = CqklE k 24 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 = (212) 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 C11 Tc12 At2 =24.5x106 /Ko Cl2 Tc A4 55.5 x 10 /K. C44 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 1009000K [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 = (213) (1+ v)( 2v) Ev cL2 = (214) (1 + v)(1 2v) E c ( = v (215) (1+ v) where E represents the Young's modulus and v represents the Poisson's ratio. As was demonstrated in Figure 24, 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 26 applied in one direction and the other directions are free to deform as in Figure 24. 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 27 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 rebroken as the body deforms under an applied force. When the applied force is removed the bonds remain in their current configurations. 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 (216) 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 (217) 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. 28 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. (218) This relationship can then be expanded to formulate the wellknown Maxwellian viscoelastic relationship between stress and strain: e=+ (219) E p where E is the elastic modulus and g. is the viscosity of the material. Figure 25 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 timedependent. 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. 29 Therefore, for the Maxwellian viscoelastic relationship, the shortterm deformation is elastic and the longterm deformation is viscous. Applied Stress t2 t Hookecn Elcstic Solid 't1 't2 Newtonian Viscous Fluid ET Mcv wellicn Viscoelastic Fluid Figure 25 Comparison in loading response among Hookean elastic, Newtonian viscous and Maxwellian viscoelastic relationships. e, = co(t,t,)/ 30 In silicon IC fabrication, oxide (SiO2) is considered to behave nonlinearly viscoelastic at midrange (80011000C) 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: filminduced, dopantinduced, 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 nonplanar 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 difference: E, = a,, AT (220) 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 32 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.54.5x106 /K over the 300900K 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. (221) 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 33 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 nonplanar 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 = .' (222) f 6(1v,)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 34 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): t a, = 4Q f . (223) ts 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. Filminduced 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 26. 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. 35 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 27 exemplifies this in three and twodimensional 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 36 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 micromachining applications [8, 39, 40], where borondoped cantilevers have exhibited bending due to strain induced by the boron. 2D planar film These nodes' stresses are (7 oriented parallel to the  deposition interface ^7 ^Y / D 7,  +v +x 2D nonplanar film nx'yC, n, *,+ n,, * x ,\ ^ y These nodes' stresses are translated around covers by the no* n,* T, component of the unit surface n,*oax * normal Y Figure 26 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 27 Two and threedimensional 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 28). Boroninduced 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 38 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, (224) where asi is the silicon lattice constant (5.4295A at 25C) and Nsi is the density of Si atoms (5.02e22 cm3). The average concentration of Boron (CB) is computed for each element from its node quantities. Dopantinduced 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. 0.0009 0.0008 0.0007 S0.0006 E 0.0005 S0.0004 g 0.0003 9 0.0002 0.0001 0.0001 0 5E+19 1E+20 1.5E+20 Boron Concentration (cm3) 2E+20 Figure 28 Hydrostatic strain as a linear function of boron concentration [41] 39 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 cm3) 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 PSi bond lengths are of the same magnitude as the SiSi bond lengths. 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 40 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] do.. pa, =d + b (225) 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 following: q' = BTad(vol) J Nbd(vol) (226) 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 41 statement is valid quite generally for any stressstrain 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). (227) V' 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 (228) would be modeled by substituting the constitutive relation into the equivalent nodal forces equation: q = [B'DEd(vol) BDoDCd(vol) + BTd(vol). (229) V' V V' The strain e is related to the unknown displacement Aa through the B matrix e = BAae (230) and may also be substituted into the elastic equivalent nodal force equation: q' = LBTDBd(vol) Aae JBTDeod(vol) + B'Taod(vol). (231) v I V V, Under mechanical equilibrium, the elastic equivalent nodal force equation reduces to the following discretized form: [BTDBA]Aae = BTADE0A BTaoA (232) 42 where A is the area (volume in 3D) of the element. The lefthand 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 (233) G r has the following solution for the stress oa a= 7 1exp( t + Coexp(F (234) where r is the relaxation time constant and is the ratio of the viscosity r7 to the elastic modulus G =. (235) G The Maxwellian viscoelastic constitutive relation can then be substituted into the equivalent nodal force equation: q' = [fB T 1exp A Dd(vol) + f BT exp (vol). (236) The strain rate e is related to the unknown change in velocity Av through the B matrix 43 e = BAv' (237) and may also be substituted into the viscoelastic equivalent nodal force equation: q' =[ B T 1expt DBd(vol) Av' + fBT o exp (vol). (238) 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 1exp(1 DBA Ave = BTa exp( At (239) where A is the area (volume in 3D) of the element. The lefthand side represents the stiffness matrix of the element. The righthand 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~ (240) The viscoelastic formulation then becomes SBT tDBA Ave = B'o(1 At (241) f T I I T ) and reduces to the elastic formulation by allowing AtAve=Aae. 44 To model twodimensional problems, the plane strain formulation is used. This formulation can be used for problems where the strain component in the zdirection is zero or negligible [28]: EZ = E = E=c =0. (242) This can be approximated solving problems with infinitely long dimensions in the zdirection. Therefore, the strain in the zdirection 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. (243) Then plane strain presumption implies that stresses in the zdirection will still occur even if there is no zcomponent 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 (244) _C 0. where v is the Poisson's ratio. 45 The FEM formulation uses three noded triangular elements called 'faces' for twodimensional applications and four noded tetrahedral elements called 'volumes' for threedimensional 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 29. 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 1Calculate the stress in the upper surface films to compute the stress dependent oxidation growth Surface Film Forces Oxide Growth__  I Forces )k Step 2Calculate the stress induced in the substrate due to the forces generated during oxidation and strain generated from dopants. _ Strain  dopants, defects Figure 29 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* (245) 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 47 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 FullyIntegrated 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 210). 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 210 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 49 from each source were described. Finally the boundary loading (BL) and fullyintegrated (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 oxidationinduced strains in threedimensions, 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 filminduced and boroninduced strain models are demonstrated and they are compared and validated with experimental studies. CHAPTER 3 APPLICATION EXAMPLES AND COMPARISONS 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 oxidationinduced stresses for different processes in silicon fabrication. Next, solutions of nitride film induced stresses are compared with microRaman spectroscopy investigations. Afterwards, boroninduced cantilever bending simulations are performed and compared with previous experiments. Finally, a three dimensional example is demonstrated for simulating oxidationinduced stress using the boundary loading method. 3.1 FEM Comparisons To compare the two methods for computing oxidationinduced stress, two different oxidation processes are analyzed. The first is a standard LOCOS process. The second is a reoxidation process over an oxide filled shallow trench. These examples detail when the boundary loading method solutions fail to match the solutions computed by the fullyintegrated method. Each of the following simulations is twodimensional and use the plane strain approximation with reflecting boundary conditions. This means that the dimensions in the zdirection (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 twodimensional field. The hydrostatic pressure definition used is simply the negative average of the normal components or negative onethird of the trace of the stress tensor: CXX 0 0 P = 0 Y Y (31) 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 31. For the STI simulations where dry oxidations are also performed, the oxide material properties were calibrated to Kao's dry oxidation cylinder Table 31 Material coefficients at 1000C [44, 46] Elastic Parameters Viscosity Parameters material Elastic Poisson's Viscosity Vo Modulus Ratio (Dyn/cm2) (Dyn/cm2) oxide 6.6ell 0.17 2.0e14 450e24 nitride 3.89e12 0.3 7.0e14 100e24 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 31. The dry oxide parameters obtained from the calibration simulations are listed in Table 32 along with the wet oxide parameters obtained by Cea [44]. Dry Oxidation Calibrations Figure 31 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 32 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 Parameters 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 32 Calibrated dry oxide viscosity compared with the wet oxide viscosity calibrated by Cea [44]. 1017 106 .... S105 1014 10'13 0.72 Table 32 3.1.1 LOCOS The LOCOS process is simulated to compare between twodimensional stress solutions computed by each method. Figure 33 exemplifies simulations of a 10000C eightyminute wet oxidation process at atmospheric pressure. A 180nm thick nitride film with a twomicron 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 33 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 34. 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 34, 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 33 represents a matrix of simulations with varying nitride linewidths. The first set with linewidth of 2gm is demonstrated in Figure 33. 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 33 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 33 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 108 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 33 Comparisons for a LOCOSinduced stress computed by the boundary loading method (left) and the fullyintegrated method (right). FI Method Y Axis 59 Oxide Shape Comparisons x in microns 0.50 0.00 0.50 0.60 0.00 y in microns 0.60 Figure 34 Overlaying oxide shape solutions comparisons between the BL method and FI method for an example LOCOS simulation with 2gm nitride linewidth. 3.1.2 PostSTI Process ReOxidation Strain Solution Comparisons. An oxidation process over a previously processed oxide filled shallow trench isolation (STI) is simulated to compare between twodimensional stress solutions computed by each method. Figure 35 exemplifies a simulation of a 10000C threeminute wet oxidation at atmospheric pressure. The trench with 0.5pm depth and 0.5gm width is filled 60 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 35 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 reoxidation 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 0108 +108 y +108 o108 +108 1.0 1.0 S109 A Compression 10 10 7 1.OxlO7 1.5 109 1.5 Tension 1.0xlO 1o8 Compression 1.oxo 107 1.0xl0 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 35 Comparisons for post STI processinduced stress computed by the boundary loading method (left) and the fullyintegrated 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]. Reoxidation 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 reoxidation process [44]. Stressdependent models for oxide viscosity, reaction rate, and diffusivity are utilized for simulating the re oxidation anneals. Twodimensional plane strain simulations are performed for a reoxidation process to study the evolution of the radius of curvature of the upper STI corners. The reoxidation occurs after the planarization step of the STI process. Therefore the oxide fills the trench up to the surface of the nitride film (Figure 36). The STI dimensions have a depth of 300nm and a width of 400nm. A series of reoxidation simulations is performed varying the nitride thickness from 70nm to 180nm over a 30nm pad oxide. The re x in microns FLOOPS Si3N4 Nitnde Liftng__ ~,' ._ 0.00 i Bird' s Bak Le gt NI Raius of Curmtwe 0.10 Si02 Silicon 020 0.30 0.20 0.00 y in microns 0.40 Figure 36 10000C 300 minute dry reoxidation 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 37). 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 tradeoff against the thinner nitride process. 70o (a) 22 ......... (b)  6020 55 18 50 4 16 45 114 40 35 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 37 (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 38 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.15 0.10 0.00 0.10 950C wet (FI) 950C Wet (BL) 0.20 1100C Dry (BL) 0.30 0.00 in microns Figure 38 Overlay of the reoxidation 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 region. 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 edgeinduced stress with measured strain data. Currently, the best method for measuring local processinduced mechanical strains in the substrate at submicron dimensions is microRaman 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]. MicroRaman 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 Rcm1. 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 processinduced strain was introduced by Jones [51] and De Wolf [52]. This method is implemented here to validate the stress solution from the nitride film edgeinduced stress model with microRaman 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. (32) 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 (Pq)e12 2rE,3 (P q)eL2 E22+ q33+ q' El 2rE23=0 (33) 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 reading. 3.2.2 Nitride Film EdgeInduced Stress A simulation of the stress induced by a nitridepoly 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 of3e9 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 39 by stress contours of hydrostatic pressure. This represents the stress state of a PolyBuffered 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 310. 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 2.00 3.00 8 4.00 4.00 5.00 0.0 y in microns 10.0 Figure 39 Simulation of the stress induced by the intrinsic stress of a nitride poly stack for a linewidth of 5gm. Units are in Dyn/cm2. f0.3 0.1 0.3 0.3 Position (cm) Figure 310 Frequency shift simulation for the three Raman modes due to nitridepoly stack edge induced strain compared with Raman measurement [52]. 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 310. 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 filmrelated effects on the incident signal reflection, or convolution averaging errors due to light intensity decay and spot width averaging. 3.3 BoronDoped Cantilever Bending Comparisons Next, an application using the boroninduced strain model is presented in order to validate its solutions. The straininduced bending of boron doped cantilevers are simulated and compared with previous experimental investigations. 3.3.1 Silicon Bulk MicroMachining Silicon bulk micromachining is important for fabricating siliconbased sensors and transducers. Silicon sensors are often composed of thin membranes, bridges, cantilevers, and beams. These structures can be fabricated by various bulk micromachining 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 cm3. It is believed that the strong BSi 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 cm3. However due to these levels of boron concentration, high levels of residual stress are generated. Since micromachined 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 311 illustrates the bulk micromachining process performed by Yang [39]. # Si02 Oxidation Si Si02 Thermally Diffused Boron Backside Etch Staininduced Bending Figure 311 Bulk micromachining process flow for fabrication of the cantilever. 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 cm3). 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 311 [54]. First boron is introduced by thermal diffusion. The resulting profile has a peak concentration of 8x1019 cm3 (Figure 312). A backside wet chemical etch is then performed. A boron concentration of 6x109cm3 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 borondoped cantilever structures. The 2d cantilever has a length of 50pm. The grid spacing in the xdirection is 0.05gm and in the ydirection is 0.5gm. Since the xdirection 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 Boro 8x10 6x10 4x10 2x10 nCoi.c(o) FLOOPS____ 20? ____ )19 *,1 19 _________________ 19? _____________ ____ 200 Dep1 (micrni) 1 00 Figure 312 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 zdirection (width), 0.5am in the ydirection (length), and O.1pm in the xdirection (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 313a. An example of the deflection simulation for the 0.62pm thickness beam is shown in Figure 314. The structure is regridded after the nodal displacements are solved for. Figure 313 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 314). 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 315. 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 313b. The same general trend also results. (a) 20 Simulation 50 um Length Deflections 80E06 T 60E06 ..: .. . ,b .. " 4 rF 0 ., 3 2 1t .... .......... .............. ............. ... ...... .... 0 000i 00o0 0003 0004 YPosition along the Cantilever (cm) 00a0 (b) 40E07 3OE07 20E07 I OE07 E e  000.00 2 0E07 0 30E07 40E07 50E07 6.OE07 3D Simulation 10 um Length Deflections s^ 1~ " a' *092  1I7 T , 1 22 137 0 0.0002 00004 00006 0.0000 0001 YPosition along Cantilever (cm) Figure 313 (a) Two and (b) Threedimensional 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 solution. x in microns FLOOPS 1.40  12'' ... '" 'i2 i'K '" .. :;'7 Figure 314 Twodimensional 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 FLOOPS .. ._ ... . .. _ 0. 000393 1 0. 000000000 Figure 315 Threedimensional simulation of the bending of 1.37um thickness cantilever. The maximum deflection is at the bottom corner of the cantilever beam. cantilever beam. 77 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 316. 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  100 1 200 S300 N 0 0 I ~400  Various Thicknesses Figure 316 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 selfloading (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]: WL2 max = (34) 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 following: W = Apg (35) where A is the crosssectional area of the cantilever, p is the density of silicon (2.33 gm/cm3 for 5el9cm3 borondoped 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 relationship: 79 6M a A. (36) At The following gravityinduced stress relationship then results by substituting (34) and (35) into (36): 3pgL2 a 3 (37) t 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 5el9cm3. 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.170.88 = R 1.17 0.88 = 0.248. (38) 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 threedimensional strain simulations. For two dimensional oxidation simulations, the BL method averaged 25 times faster depending on the number of elements in the silicon. This time savings is even much greater for threedimensional simulations since the number of elements for a typical 3D simulation is an order of magnitude greater than in 2D. Figure 317 illustrates a threedimensional 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 317 Threedimensional simulation of LOCOSinduced stress in the substrate. film are 1.0pm x 1.5pm x 0.15pm. The isostress 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 isocontour lines. Figure 318 and Figure 319 are top views of the same simulation. The stress tensor component in the xdirection (directed normal to the surface) is Figure 318 YZ plane view of the tensile oxx component of stress induced during LOCOS oxidation. plotted as contour lines. Figure 318 demonstrates that the maximum tensile stress is at the corner of the nitride film. Figure 319 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 oxidationinduced stress computed by the BL method and FI Figure 319 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 reoxidation 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 edgeinduced stress of a nitride poly stack was demonstrated. The results of these simulations were then used to predict the Raman signal measured over the nitridepoly 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 boroninduced strain model was used to simulate the bending of borondoped cantilevers due to their boroninduced 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, threedimensional 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. CHAPTER 4 KELVIN PROBE FORCE STI EXPERIMENT 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 expansioninduced stress of the sidewall oxidation (Figure 41). 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 41 Major sources of stress in STI fabrication include the following: (a) Nitride film edgeinduced, (b) Insulator thermal expansion mismatch induced, and (c) Sidewall oxidation volume expansioninduced stress. Arrows indicate force vectors. generation process design, methods are needed to measure the stress quantitatively. Up until now, microRaman 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 stressinduced 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 [5762]. This has been achieved by detecting the electrostatic potential difference (EPD) AFM Cantilever Probe I I +z \V v V Rastered laterally FE(Z) Figure 42 Crosssectional sketch of SKPM measurement. between the SKP tip and sample (Figure 42) through their electrostatic and van der Waals force interaction. The method involves rastering a cantilever probe over a crosssectional 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 cantileversample 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 zdirection 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 42. Laboraory STI Array Sample Laboratory Cross Section Frame of Reference Z [iTO] Y [110] x [00T] Figure 43 Sketch of Kelvin probe scanning laterally (ydirection) and in depth (xdirection). 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 flatband voltage bias is then applied to the cantilever probe to null the electrostatic force. This DC bias is proportional to EPD of the cantileversample system. The DC flatband voltage is then recorded at each raster step and provides a twodimensional map of the work function variation of the sample (Figure 43). 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 44). It is similar to that of a metalinsulatorsemiconductor (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 surface. In SKPM mode, a DC flatband bias is applied to null the electrostatic force (Figure 45). This applied bias represents the EPD at each rastered TiSilicide Tip vacuum Silicon Figure 44 Energy band diagram of the cantileversample system at zero applied DC bias. point along the crosssectional surface of the sample. The measured VEPD can then be expressed as the sum of the flatband correcting voltage VFB and surface potential terms [61]: VEPD = Vf +V, + sE s,]. Coins (41) The flatband voltage is the difference between the work functions of the tip and sample: (42) The work function of the sample Osi is the energy difference between the vacuum level and the Fermi level: TiSilicide Tip vacuum Silicon qXsi I EC AX  Figure 45 Energy band diagram after VEPD is applied and the electrostatic force is nulled. V = TiSi Os,. 