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

Full Text 
PRAGMATIC AND RELIABLE DEVICE/CIRCUIT SIMULATION FOR DESIGN IN ADVANCED SILICONBASED TECHNOLOGIES By DUCKHYUN CHANG 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 1997 ACKNOWLEDGEMENTS I would like to express my deep gratitude and appreciation to my advisor, Professor Jerry G. Fossum, for his devoted guidance, patient encouragement, and unceasing support in my academic and personal affairs throughout my graduate study. His invaluable guidance made this work possible. My sincere gratitude goes to Professor Mark E. Law for his helpful discussions concerning FLOODS. I also would like to thank the members of my supervisory committee, Professors Sheng. S. Li, Kenneth. K. O, and ChenChi Hsu, for their willing service and guidance on my committee, and Mary Turner for all of her help throughout my graduate career. I would like to acknowledge and thank IBM, NCCOSC, and the Semiconductor Research Corporation for their technical and financial support throughout my research. Especially, I feel obliged to express my deep thanks to Mr. Mario Pelella and Scott K. Reynolds at IBM for providing technological information invaluable to the integrity of this work. I would like to recognize in this achievement Dr. GhyBoong Hong and Mr. MingYeh Chuang, who helped me in many ways with the sincere and profound discussions on the related topics. Gratitude is also extended to many of my friends, who helped me through technical discussions or by making my life one of the most cheerful. I cannot mention all, but I would like to mention Drs. Dongwook Suh, Hyunsoo Kim, Ping Chin Yeh and MinKyu Lee, and Messrs. Yochul Ho, Kihong Kim, Soonchul Park, Namkyu Park, Minbo Shim, Keunwoo Kim, Glenn Workman, Doug Weiser, MengHsueh Chiang, and Jonathan Scott Brodsky. ii I also thank my parents, fatherinlaw, Dr. Jong Duk Lee, motherinlaw, brothers and sisters for their loving encouragement through the long period. Simple words seem to diminish their unconditional and unending love and support. Finally, I would like to share this work with my lovely wife and friend, Sohee, and my dearest son Alex, without whom the long years of my graduate study might have been so simple and barren. iii TABLE OF CONTENTS ACKNOW LEDGM ENTS .................................................... ...................... ii ABSTRACT ............................................................. vi CHAPTERS 1 INTRODUCTION ..................................... ................................... 1 2 SIMPLIFIED ENERGY BALANCE MODEL FOR PRAGMATIC MULTIDIMENSIONAL DEVICE SIMULATION ...................... 8 2.1 Introduction .............................................................................. 8 2.2 HotCarrier Transport Theory ...................................... 11 2.2.1 Carrier Transport Equations ...................................... 11 2.2.2 NonLocal carrier Heating ...................................... 15 2.2.3 Hierarchy of Transport Models ..................................... 16 2.3 Simplified Energy Balance Model ..................................... 19 2.3.1 Model Development and Its Implementation in FLOODS ........................................ ............ 19 2.3.2 HotCarrier Effects .......................................... ....... 26 2.4 M odel Verification ................................................................ 27 2.4.1 Bulk M OSFET ........................................ .......... 27 2.4.2 SOI M OSFET.................................... .............. 37 2.5 M odel Extension............................ ...... ................. 39 2.6 Sum m ary ..................................... ......................... 43 3 IMPLEMENTATION OF TEMPERATURE DEPENDENCE IN MMSPICE HBT/BJT MODEL ......................................... 45 3.1 Introduction............................................... 45 3.2 Si and SiGebase BJT Characteristics in Temperature....... 47 3.2.1 Collector Current ......................................................... 47 3.2.2 Base Current.............................. ..... ............... 51 3.2.3 Current Gain................................. ..... ............. 52 3.2.4 Output Conductance............................. ........... 53 3.2.5 Base Transit Time ....................................... ..... .. 54 3.2.6 Emitter Transit Time .......................................... 55 3.2.7 Cutoff Frequency ....................................... ...... .. 57 3.3 MMSPICE Implementation.......................... ........... 58 iv 3.3.1 Temperature Dependences of HBT/BJT Model Parameters..................... ............ 58 3.3.2 Model Implementation in MMSPICE3.0.................. 68 3.4 Verification and Discussion ................................................. 69 3.5 Sum m ary ..................................... ......................... 80 4 PHYSICAL SOI MOSFET CHARGE MODEL FOR LOWVOLTAGE MIXEDSIGNAL CIRCUIT SIMULATION ....... 82 4.1 Introduction.................................................................. 82 4.2 SOI MOSFET Charge Models ...................................... 83 4.2.1 NFD/SOI ......................... ............. 85 4.2.2 FD/SOI .................................................... 95 4.2.3 Charges and Transcapacitances ................................ 97 4.3 Assessment of SOI Technologies Based on Unified FloatingBody Model ..................................... 100 4.4 Summary ...................................... 110 5 SOI CMOS CIRCUIT SIMULATION AND DESIGN ISSUES..... 112 5.1 Introduction............................ 112 5.2 Model Calibration ..................................... 113 5.3 10bit SOI DAC Design ....................................................... 119 5.3.1 Capacitance Characteristics of the SOI MOSFETs ... 119 5.3.2 Current Cell Design............................... 123 5.3.3 DAC Simulation Results ..................................... 130 5.4 FD DRAM Cell Design ......................................................... 138 5.5 Summary ...................................... 140 6 SUMMARY AND SUGGESTIONS FOR FUTURE WORKS ....... 145 APPENDIX AN AUGMENTED DRIFTDIFFUSION MODEL FOR EFFICIENT NUMERICAL DEVICE SIMULATION ................... 150 REFERENCES ........................................ 165 BIOGRAPHICAL SKETCH ..................................... 172 V 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 PRAGMATIC AND RELIABLE DEVICE/CIRCUIT SIMULATION FOR DESIGN IN ADVANCED SILICONBASED TECHNOLOGIES By DUCKHYUN CHANG AUGUST 1997 Chairman: Dr. Jerry G. Fossum Major Department: Electrical and Computer Engineering This dissertation describes pragmatic applicationspecific mixedmode device/circuit simulation approaches for design in contemporary Sibased technologies. First, a methodology for improving the efficiency of conventional numerical device simulation is proposed, where a new simplified energybalance (SEB) model is developed and implemented in FLOODS as a pragmatic option. In the SEB model, the energyrelaxation length is estimated from a preprocess driftdiffusion simulation using the carriervelocity distribution predicted throughout the device domain, and is used without change in a subsequent simpler hydrodynamic (SHD) simulation. The new SEB model is verified by comparison of twodimensional SHD and fullHD DC simulations of submicron bulkSi and SOI MOSFETs. The most noteworthy feature of the new SEB/SHD model is its computational efficiency, which results from reduced Newton iteration counts due to the enhanced linearity. Next, a seminumerical approach is considered for mixedmode device/circuit simulation where physical compact models are developed based on device structures. The approach is first applied to advanced SiGebase HBTs through MMSPICE to investigate their temperature dependences. MMSPICEsimulated and measured DC and vi AC characteristics are shown to be in good agreement for SiGe HBTs, over a wide temperature range, thus proving the practical utility of the seminumercial approach. The approach is next applied to submicron fully depleted (FD) and nonfully depleted (NFD) SOI MOSFETs via SOISPICE to analyze floatingbody effects on circuit performance. New continuous charge models for scaled devices, having continuous derivatives in all regions of operation as well, are presented. The floatingbody effects are physically accounted for in the models and unified for DC, AC, and transient simulations. Several applications of SOISPICE are described to exemplify unique SOI circuit behavior due to floatingbody effects and to show how it necessitates physical models for reliable circuit simulation and design. FD/SOI DRAM design is first examined as an alternative to avoid possible failure mechanisms of NFD/SOI DRAM, such as transient leakage currents in the memory cell and data instability in the sense amplifier. In another circuit application, a simulationbased design of a 10bit currentsteering digitaltoanalog converter (DAC) is described to demonstrate a novel approach to ensure kinkfree operation of floatingbody NFD/SOI analog circuits. vii CHAPTER 1 INTRODUCTION Lowpower and highspeed integratedcircuit (IC) technology is the key to future portable and wireless communication systems [Ish97]. Scaling of conventional complementarymetaloxidesemiconductor (CMOS) and bipolar devices has been the common way to advance the technology. However, as the devices are scaled towards tenthmicron and smaller dimensions, a significant tradeoff between performance and reliability has to be made. For example, scaling of the MOS device to deepsubmicron regimes (e.g., 0.1 gm) imposes tremendous challenges to device designers mainly due to hotcarrier effects (HCEs) and shortchannel effects (SCEs). Hot carriers generated by the very large electric fields in scaled devices cause severe reliability problems and present major constraints on design, even under normal operating conditions [Ago94]. The hot carriers can be injected into the oxide and trapped, which leads to significant thresholdvoltage shift and transconductance degradation as well as oxide breakdown [Hwa89]. The substrate current measured in MOSFETs provides a proof of the occurrence of such HCEs. The hotcarrierrelated reliability problems can be improved with reduced power supply voltage, but then SCEs become a dominant factor limiting device performance in lowvoltage/lowpower MOSFET applications. A traditional approach for alleviating SCEs is to increase the channel doping concentration, which increases the threshold voltage and suppresses the subthreshold conduction. However, such high concentration degrades device 1 2 performance due to decreased mobility and increased junction capacitance [Yan92]. Use of thinfilm SiliconOnInsulator (SOI) technology has been widely investigated as a promising alternative to ease the tradeoff of scaled CMOS design. The SOI technology shows superior performance to bulkSi technology due to the reduced source/drain parasitic capacitance, simpler device isolation, absence of latchup, increased radiation hardness, feasibility of resistors and capacitors free of parasitic junction effects, and less noise coupling afforded by highresistivity substrates [Yeh95, Suh95]. In fact, various circuits utilizing deepsubmicron SOI technology have been recently fabricated, and its advantages have been demonstrated [Ued97]. SOI as a mainstream technology is coming in circuit applications like DSP and microprocessor used in personal communication systems, gigabit DRAM and SRAM, and microwave RF circuits. Also, the recent advances in bipolar technology, such as sustained device scaling, lowtemperature epitaxialbase process [Gan90], and Si/SiGe heterojunction [Kon96, Jos96] have significantly extended the utility of the Sibased bipolar junction transistor (BJT). The graded SiGebase heterojunction bipolar transistor (HBT) provides improvement in many key analog design parameters (e.g., P, ft and fma) compared to the conventional Si BJT, and is thus of interest for advanced microwave RF circuits. With these advanced Sibased technologies like SOI CMOS and SiGebase HBT, the performance of the scaled circuits can be significantly enhanced beyond the limitations imposed by the conventional technology. To maximize the performance improvement and reduce the development cost, an efficient and reliable IC technologyCAD (TCAD) system is required. The conventional IC TCAD system consists of purely numerical process (e.g., SUPREME) and device (e.g., MEDICI) simulators, and a circuit simulator with empirical 3 compact device models like BSIM3 (for MOSFETs) and GummelPoon (for BJTs). Numerically intensive device simulation, based on the process simulation, is done for device characterization, and then modelparameters are optimized for the circuit simulation through rigorous parameterextraction processes. To effect a design, this algorithm is iterated until optimal circuit/ device performance is predicted from a given process flow. Severe drawbacks of the conventional TCAD system are the extremely long turnaround time and loss of correlation between predicted circuit performance and device structures and process flow. For example, the thinfilm nature of the SOI MOSFET can underlie mechanisms that complicate circuit simulation and portend equivocation in design [Fos97]. Mainly because of floatingbody effects, which in fact obtain even in bodytied structures because of unavoidable high resistance, empirical compact models like those used for designing bulkSi CMOS circuits may not be useful for SOI. Not only are such models deficient, but they furthermore cannot be easily calibrated using common parameteroptimization techniques because of equivocal data acquisition implied by SOI device selfheating in DC measurements and the floatingbody effects in pulse measurements designed to avoid the selfheating [Fos97]. The efficiency and reliability of the IC TCAD system could be greatly enhanced by a pragmatic mixedmode device/circuit simulator based on physical, compact models which have processbased parameters that relate directly to the device structure. With this objective, MMSPICE (for bulkSi and SiGebase bipolar and CMOS) [Jeo89, Jin92, Hon94] and SOISPICE (for SOI CMOS) [Vee89, Suh95, Yeh95] were created. Here, we further extend the processbased device models for general circuit applications to demonstrate the utility of this seminumerical mixedmode simulator. Further, the resulting enhanced models are used in various device/circuit simulations in actual technologies to 4 demonstrate viability and to reveal unique simulation and design issues for scaled SOI CMOS and SiGebase HBT circuits. In Chapter 2, a purely numerical way to improve the efficiency of the conventional TCAD system is examined. To pragmatically account for nonlocal carrier heating and hotcarrier effects such as velocity overshoot and impact ionization in multidimensional numerical device simulation, a new simplified energybalance (SEB) model is developed and implemented in FLOODS [Lia94a] as a hierarchical option. In the SEB model, the energyrelaxation length is estimated from a preprocess driftdiffusion simulation using the carriervelocity distribution predicted throughout the device domain, and is used without change in a subsequent simpler hydrodynamic (SHD) simulation. The new SEB model is verified by comparison of twodimensional SHD and fullHD DC simulations of submicron bulkSi and SOI MOSFETs. The SHD simulations yield detailed distributions of carrier temperature, carrier velocity, and impactionization rate, which agree well with the fullHD simulation results obtained with FLOODS. The most noteworthy feature of the new SEB/ SHD model is its computational efficiency, which results from reduced Newton iteration counts due to the enhanced linearity. Relative to fullHD, SHD simulation times can be shorter by as much as an order of magnitude since larger voltage steps for DC sweeps and larger time steps for transient simulations can be used. The improved computational efficiency can enable pragmatic threedimensional SHD device simulation as well. In Chapter 3, a detailed analysis of temperature dependences in Si and SiGebase HBT/BJTs is presented and implemented in MMSPICE to demonstrate the utility of seminumerical coupled mixedmode device/circuit simulation for bipolar ICs, and to gain physical insight regarding the performance of advanced bipolar transistors in temperature. Conventional Si bipolar 5 transistors are generally not considered for operation at low temperature because of reduced current gain, increased base resistance, and degraded transit time due to lower mobility [Dum81]. The dominant reason for the current gain reduction is larger bandgap narrowing in the heavily doped emitter than in the base. However, for the scaled Si BJTs, the base doping level must be high to prevent punchthrough as the base width is thinned to minimize the base transit time. The high base doping level combined with the polysilicon emitter contact leads to the improvement in current gain at low temperature [Woo88]. Also, by incorporating Ge into the base, further improvement in current gain can be obtained since Ge introduces additional bandgap narrowing in the base. In fact, for the recently developed SiGe HBT having larger bandgap narrowing in the base than in the emitter, the current gain actually increases as temperature decreases. It is shown that key device characteristics such as current gain and cutoff frequency have different temperature dependences according to the bandgap narrowing and doping gradient in the base. Although most concerns are given to lowtemperature operation, the device operation above room temperature is also examined for hightemperature applications. For verification, MMSPICE simulations in temperature are performed and compared with measurements based on recently developed HBT/BJT technologies. In Chapter 4, the seminumerical approach to mixedmode simulation is applied to submicron fully depleted (FD) and nonfully depleted (NFD) SOI MOSFETs via SOISPICE to analyze the floatingbody effects on circuit performance. SOISPICE was previously developed mainly for digital applications and verified as a reliable simulation tool for lowvoltage and highspeed digital circuit design [Fos97]. Its physical, processbased models can be effectively used to characterize FD and NFD SOI MOSFETs. However, to apply the 6 models to analog and mixedsignal circuit design, model enhancement is necessary because, unlike digital circuits, where performance depends on relatively few circuit parameters such as threshold voltage and drive current, the performance of the analog circuit is quite sensitive to the smallsignal parameters as well as largesignal characteristics of the SOI MOSFET. New charge models for scaled SOI devices are presented, where all current and charges and their derivatives, which form the basis of the SOISPICE quasistatic formalism, are continuous in all regions of operation. The characterization of body charge in the FD model accounts for accumulation charge, which means that the model is hybridFD/NFD with regard to charge dynamics. In the NFD model, the body charge is defined by charge neutrality, which means that the charge dynamics properly reflects all the intrinsic transconductances. The unified floatingbody model is described based on the charge models, and is shown to predict well the frequencydependent floatingbody effect as well as DC and transient floatingbody effects. In Chapter 5, the upgraded SOI MOSFET models are calibrated and used in circuitsimulation applications in actual technologies to reveal simulation and design issues for scaled SOI CMOS. First SOI DRAM design is described. SOI has been getting attention for the nextgeneration DRAM technology (e.g., 1 Gb DRAM) due to its lower power consumption and high speed which are critical for the application to portable electronics. NFD/SOI is dominant technology for this application, but dynamic floatingbody effects such as transient leakage and threshold hysteresis can cause potential failure of the normal DRAM operation [Fos97]. Here the FD/SOI DRAM circuit is examined, checking if it can avoid the dynamic floatingbody effects. The dynamic floatingbody effect associated with the sense amplifier can be solved using the FD device because it effectively decouples the body voltage from the 7 channel current. However, like the NFD device, transient leakage current is observed due to the capacitance coupling. Next a 10bit currentsteering DAC design is presented. FD or bodytied NFD SOI devices have been used almost exclusively for analog/mixedsignal applications because they suppress the kink effect, which degrades the gain of amplifiers and the constancy of current sources [Red96]. However, the small process margin of FD devices with regard to thresholdvoltage control, and the area penalty of body ties, as well as the high inherent body resistance and capacitance introduced by the body contact, are problematic [Cha95]. The floatingbody NFD/SOI MOSFET has not been seriously considered as an option for analog circuits since it suffers from the kink effect. Recently, NFD/SOI circuitdesign approaches to avoid the kink effect, i.e., highfrequency operation [Red96, Sin97], cascode configuration [Red96, Red90], and linked bodies [Che96, Koh97] have been reported, but none has been demonstrated to be generally useful for analog and mixedsignal circuits. Here a 10bit currentsteering DAC based on floatingbody NFD, bodytied NFD, and FD SOI MOSFET technologies is presented, and a novel approach to ensure general kinkfree operation of floatingbody NFD/SOI analog circuits is proposed. SOISPICE, with its device models upgraded and refined for highfrequency analog applications as described in Chapter 4, is used to simulate the DAC and verify the design approach. In Chapter 6, the main accomplishments of this dissertation are summarized, and future research tasks are suggested. In the Appendix, with the aim of Chapter 2, an efficient onedimensional hotcarrier transport model is developed and implemented in FLOODS to account for nonlocal carrier heating and velocity overshoot for device simulation. The model is verified through simulation of several advanced device structures. CHAPTER 2 SIMPLIFIED ENERGYBALANCE MODEL FOR PRAGMATIC MULTIDIMENSIONAL DEVICE SIMULATION 2.1 Introduction In contemporary deepsubmicron devices, nonstationary and nonlocal electron transport phenomena make prominent impact on device performance and reliability. Since they are fundamentally described by the electron energy distribution, not by the electric field, the conventional driftdiffusion (DD) model is not suitable for device design in spite of its maturity and computational efficiency. The Monte Carlo (MC) analysis of microscopic carrier transport is considered as the most appropriate simulation method since it solves the Boltzmann transport equation (BTE) directly and thus provides the detailed electron energy distribution [Ste93]. But its statistical nature leads to extremely long computation time which limits its utility for engineering applications [Scr94]. The hydrodynamic (HD) model, in which carrier dynamics is described by averages of energy and velocity, has been widely investigated as an alternative to the MC method due to its capability for describing nonlocal carrier transport in less computing time [Str62, Blo70, Tan93, Tho91]. However, the overall accuracy of the HD model is questionable due to the equivocation in the implicit physical modeling, e.g., in the characterization of momentum and energy relaxation times and their dependencies on carrier energy [Go188], and to the difficulties in the proper model incorporation for nonparabolic energy 8 9 band [Bor90], nonMaxwellian energy distribution [Che92], and higherorder moments of the BTE (e.g., heat flux) [Bor91]. Also, it suffers from relatively much longer simulation runtime compared with the DD model. Indeed users of a contemporary HD numerical device simulator can wait hours for a simulation to executed, and have little confidence in the results. Alternatively, the augmented driftdiffusion (ADD) model (Che91], in which nonlocal effects are accounted for by introducing additional analytic fielddependent terms in driftdiffusion current equations, has been proposed as a way of efficiently extending the utility of driftdiffusion simulation. Since it only solves the carrier continuity equations with Poisson's equation, its simulation run time is comparable to the DD simulation. In the ADD model, the electron temperature is calculated from the local electric field and field gradient. Even though this model has been partially verified as an efficient method to predict velocity overshoot [Che91], it loses its validity when the second, or higherorder derivatives of the electric field are significant, which is not uncommon in small devices. In the Appendix, we present an enhanced ADD model in onedimensional form and verify it for numerical simulation of scaled devices. To obtain an adequate description of the average electron energy more efficiently than via the HD model, simplification of the energybalance equation has been recently proposed for postprocessing use. The postprocessing technique, in which electron energy is calculated from the electric field distribution obtained from a DD simulation, has been primarily used to predict the energydependent impactionization rate [Cra90, Slo91, Sou93, Ago93]. Slotboom et al. [Slo91] simplified the energybalance equation by assuming a constant (saturated drift) velocity throughout the device, and then used it to model the impact ionization along a single current path in MOSFETs. 10 However, the model does not provide detailed information about the twodimensional distribution of the carrier energy within the devices. Recently, Agostinelli and Tasch [Ago93] demonstrated an energydependent twodimensional impactionization, or substrate current model for the simulation of submicron MOSFET's, in which the simplified onedimensional energybalance equation proposed by Cook [Coo82] is applied to many current contours in order to generate a twodimensional representation of carrier energy. However the average energy predicted by this postprocessing technique is not correct in highly scaled devices since it decouples the electric field from the energy. For example, in deepsubmicron MOSFETs, velocity overshoot influences the electricfield distribution for a given bias condition, and effectively defines a higher drain saturation voltage, which in turn defines a higher current and affects gate charge, or capacitance, as well. This effect further reduces the impactionization current because of the lower electric field near the drain. The significant electric fieldenergy coupling hence tends to invalidate the postprocessing method of characterizing the substrate current in scaled MOSFET's. Furthermore, issues related to the implicit simplified modeling of the carrier velocity [Slo91, Coo82] in the energybalance equation prohibit use of these methods in general device simulation. The simplifications are not appropriate for regions in which the net flow of the carriers is small. According to these models, the carriers are heated, even at equilibrium, due to large builtin electric fields, for example at the baseemitter junction of a scaled BJT. This is contrary to the definition of thermal (or quasi) equilibrium. In such cases, the carriers are not heated because the drift current is balanced out by the diffusion current. For a MOSFET, this error can be more serious for all bias conditions. In the vicinity of the inversion layer, the noted simplifications 11 predict high carrier energy due to the large transverse electric field and field gradient; but there is no carrier heating due to that field. Hence, the inadequate modeling of the carrier velocity in the energybalance equation prohibits reliable device simulation, and indeed could cause numerical errors in simulation as well. The purpose of this chapter is to develop a reliable simplified energybalance (SEB) model which accounts for the nonlocal effects with reasonable accuracy. The new SEB model is formulated for practical multidimensional device simulation and is implemented in the twodimensional tool FLOODS (FLorida ObjectedOrient Device Simulator) [Lia94a] as a pragmatic option. The physical models underlying the SEB formalism are also discussed. Nonlocal effects in the submicron bulksilicon and SOI MOSFETs, such as velocity overshoot and impact ionization, are analyzed based on simulation results, and are compared with those of the full HD model in FLOODS. The superior computational efficiency of the new model, relative to the HD model, which renders it pragmatic and attractive for threedimensional and transient device simulation, is addressed as well. 2.2 HotCarrier Transport Theory 2.2.1 Carrier Transport Equations The energy distribution of hot carriers plays a crucial role in predicting the longterm reliability of advanced devices. For example, in a deepsubmicron MOSFET, the substrate current, the gate current, and the oxide breakdown are in general determined by the high energy tail of the carrier energy distribution [Fie93]. The solution of the Boltzmann transport equation (BTE) gives the detailed information on the carrier distribution, and thus momentum or energy 12 distribution. Averages of these distributions provide carrier concentrations, current densities, and carrier temperatures throughout the device. Once the distribution function is found, a complete device description is possible. The BTE can be derived from the semiclassical dynamics, under the assumptions that carriers behave as classical particles and collisions are instantaneous [Lun90]: af +* (Vrf)+Fe (pf) = afc (21) where f is the carrier distribution function, u is the carrier velocity, Fe is the external force, and p is the momentum of the carriers. The terms on the left side are often referred to as the drift terms, and the term on the right side as the collision integral. The drift terms represent the change of the distribution function with time, the flow of carriers in the real space, and the flow of carriers in the momentum space; and the collision integral represents the variation rate of the distribution function due to the scattering, which can be calculated by integrating probabilities of transition. Since it is very difficult to solve the BTE directly and sometimes it is undesirable, simpler expressions such as balance equations are often used to analyze the performance of the devices [Lun90]. Virtually an infinite series of balance, or moment equations can be derived from the BTE. If we want to obtain the same information that the BTE gives, the infinite set of balance equations should be solved. However, fortunately, the first few moments, i.e., the zeroth, first, and secondorder ones that relate to the average carrier density, velocity, and energy, respectively, are usually sufficient to characterize semiconductor devices. The zerothorder moment equation is the familiar carrier continuity equation. In the steadystate condition, this equation for electrons is given by 13 Ve(n ,) + (Rn Gn) = 0 (22) where n is the electron concentration, and (RnGn) is the recombination/ generation rate; un now represents the average electron velocity. Note that (22) physically describes the conservation of the electron density. The firstorder moment equation is usually referred to as the momentumbalance equation. In the steady state, this equation for electrons becomes kB kTnVn (kBn) (23 4n  n E n + V (23) where Rn is the electron mobility, E is the electric field, and Tn is the average electron temperature. It states that the change of momentum is due to the forces exerted by the electric field and the electron pressure, nkBTn. Note that if we neglect the temperature gradient term in (23) and assume the fielddependent mobility model, the momentumbalance equation reduces to the well known driftdiffusion equation which has been widely used in device simulation. However, since scattering processes are fundamentally described by the carrier energy, not by the electric field, the driftdiffusion model is not adequate for investigating modern scaled devices. A better approach is to include the temperature gradient term in the momentumbalance equation and to use the energydependent mobility model, as is done in (23). For this purpose, we should solve the first three moments of the BTE rather than just the first two. The additional equation is the energybalance equation, which is derived after the BTE is multiplied by mi1 and integrated over momentum space. This equation for electrons, under the steadystate condition, is given by S* Vwn = (n E) V(nkgTnn + Q) n (24) n 'r." 14 where wn is the average electron energy, Q is the heat flux, and Tw is the energy relaxation time. In (24), wn and wo are defined as 3 1 2 Wn = kBTn + m n (25) wo = kBTL (26) where TL is the lattice temperature and m* is the effective electron mass. The lefthand side of (24) is the total rate at which energy changes within a differential volume element of electron gas as it moves along the direction of electron flow. The righthand side describes how individual forces, such as the energy supplied by the electric field, the work performed by the electron gas, the energy consumed by the conduction of heat, and the energy lost by the collisions contribute to change the total energy of the volume element. Note that no matter how many moment equations we write, they always contain one or more unknown higherorder moments. The zerothorder moment equation includes the firstorder moment, the average velocity. The firstorder moment equation includes the secondorder moment, the carrier average energy, or temperature. And the secondorder moment equation includes the heat flux, the thirdorder moment. Therefore, to analyze semiconductor devices by using balance equations, the system of equations must be truncated and the remaining equations simplified. For example, the heatflux term, which represents the contribution of the third moment to the energybalance equation, is commonly assumed as the phenomenological Fourier law, Q = KVTn (27) 15 where K is the thermal conductivity. However, this truncation and simplification often cause nonphysical numerical results such as spurious velocity spikes [Lun90]. And thus, it should be noted that even though this momentequation approach is far simpler for characterizing the devices than the directsolution approach of the BTE, the accuracy and reliability of this approach is strongly dependent on how the infinite series of balance equations are truncated and how the resulting higherorder unknowns are physically modeled. The DD and ADD models are other examples of truncation and simplification, in which the carrier energy is implicitly calculated from the first two moment equations. 2.2.2 NonLocal Carrier Heating Free carriers in a semiconductor are heated by the electric field driving a current so that their temperatures are substantially different from that in the thermal equilibrium state, the lattice temperature. The conventional fielddependent hotcarrier model has long been used to examine the transport of semiconductor devices. However, the rapid spatial variation of the electric field prevents the average carrier energy from following a local dependence on the electric field. In fact, the average carrier energy is not only a function of the local electric field but also depends on how the electric field varies in space. This nonlocal carrier heating can be easily explained quantitatively by using a simplified onedimensional energybalance equation proposed by Slotboom [Slo91]. The equation can be derived from (24) by assuming that the effect of the heat flux and diffusion are negligible and v, = U1sat: dwn W O ( 3 + = E (28) dx kw 5 16 where the energy relaxation length X which is defined as 5t, sat/3, is set to a constant. If we iteratively solve (28) for the carrier energy, we obtain 3 d (3 d d W (WwEE XWx( E wwE+(... 3 3 2 dE 3 3 d2E (3) = E 5 wx + 5Xw dx2 +0 (29) As can be seen in (29), the average carrier energy is described by an infinite series of the field derivatives as well as the local electric field. The local fielddependent transport model (which retains only the first righthandside term in (29)) widely used in current device simulation is hence strictly valid only when the local electric field is much larger than the sum of its derivatives. 2.2.3 Hierarchy of Transport Models The incorporation of the physics in carriertransport models involves a tradeoff between efficiency and accuracy. Different transport models are distinguished by how the underlying physics is modeled [Lin93]. Fig. 2.1 shows the hierarchy of the transport models. The top of the hierarchy is the Monte Carlo simulation, the most accurate numerical technique available for analyzing transport in devices. It solves directly the BTE with Poisson's equation. In the MC model, fullenergy band structure is physically modeled and scattering processes are simulated via random number generation. But, in spite of its physical accuracy, it has several limitations when applied to global device simulation due to its well known statistical uncertainty. For a semiconductor under low applied field, the expected drift velocity can be much smaller than the statistical noise. This is the reason why the MC method is mainly used in investigating the highfield carrier transport in devices. It is also difficult to 17 Monte Carlo Full Energy Band Structure] Scattering Boltzmann Transport Eq. Poisson Eq. Physical Accuracy NonMaxwellian NonParabolic Band Hydrodynamic SHeat Flux Carrier Continuity Eq. Degeneracy Momentum Balance Eq. Energy Balance Eq. Energy Relaxation Time Poisson Eq. EnergyDependent Mobility Computational Simplified Energy Balance Efficiency Augmented DriftDiffusion FieldDependent Mobility Carrier Continuity Eq. Current Eq. Poisson Eq. Fig. 2.1 The hierarchy of numerical simulation. 18 treat the carrier recombination/generation because the time scales involved are much longer than the relaxation times that characterize transport in MC simulation [Lun90]. Because of such limitations, the hydrodynamic model, in which the first three moments of the BTE are solved with Poisson's equation, is frequently used as an alternative to the MC model. By accounting for the average carrier energy explicitly, the hydrodynamic model can effectively describe nonlocal transport and hotcarrier effects. The underlying physics associated with band structure and scattering mechanism is absorbed into two parameters, energy dependent mobility, or momentum relaxation time, and energy relaxation time. However, the simplifying assumptions in the conventional hydrodynamic model, such as parabolic band structure and Maxwell energy distribution, should be modified in highly scaled devices since for highenergy electrons, the energyband structure becomes nonparabolic and the energy distribution is not Maxwellian. These modifications introduce additional terms into the hydrodynamic models since they basically depend on higherorder moments of the distribution function. Note that these terms considerably increase the nonlinearity of the HD model, and thus computation time. As a pragmatic option of the HD model, simplified energybalance models, in which the carrier temperature is calculated from a given electric field distribution, have been widely investigated in the last few years. However, most of the works in this area rely on onedimensional post processing in which the coupling between the carrier energy and the electric field is ignored. We believe this decoupling increases the peak electric field and thus fails to predict the peak carrier temperature, which determines the energydependent impactionization rate. Therefore, in the next sections, we develop a new coupled twodimensional simplified energybalance model and demonstrate its utility. 19 2.3 Simplified EnergyBalance Model 2.3.1 Model Development and Its Implementation in FLOODS We here propose a new simplified energybalance (SEB) model with reliable accuracy, still maintaining the computational efficiency like the postprocessing model proposed by Slotboom [Slo91]. We first assume that the drift energy is only a small part of the total kinetic energy; hence the average electron energy is approximated as wn = kBTn. Under most conditions in MOS devices, this assumption is generally valid, although near the source region, for example in a 0.2gm device, the drift energy can be nearly 15% of the total kinetic energy [Ste93]. We next assume that the heatflux term Q is negligible in determining the electron temperature. This term in (24) represents the third moment of the BTE and is thus an additional unknown [Bor91]. It is generally approximated by the Fourier law, but in large temperaturegradient regions, this approximation greatly overestimates the heat conduction and thus leads to nonphysical results [Lun90]. To resolve this problem, a more sophisticated model for the heat flux has been presented [Ste93]; however, numerical difficulties caused by it provide an additional computational burden. Because of this modeling complexity and the fact that Q is often relatively small, the heatflux term is frequently and justifiably neglected in the energybalance model [Go188, Slo91]. Now focusing on (24), we assume that V*(nn) 0, which, as (22) shows, is equivalent to neglecting carrier recombination/generation in regions where carriers are heated. Then (24) simplifies to (VTn TL) 2qn E. (210) n n 5 rw 5k n 20 Expressing the electric field on righthand side by the gradient of the conduction band edge EC, which renders (210) applicable for homostructures as well as heterostructures, we rewrite the steadystate SEB equation as V Tn + + "E~ ] (TnTL) = 0. (211) In (211), we have implicitly assumed that the effect of diffusion is negligible since the electrons are heated mainly by drifting. For numerical device simulation, (211) can be spatially discretized using the box method. For the twodimensional case on a triangular mesh structure as illustrated in Fig. 2.2, the discretization gives lij( E + Tnwn, ij+A(T, To)i = 0 (212) edge 5k B where 1, is the length of the box side bisecting the triangle edge between mesh points i and j, Ai is the box area around point i, and Xwn, ij is the electron energyrelaxation length defined on triangle edge (ij) of the grid, which is perpendicular to lij. The energyrelaxation length, first defined by Slotboom [Slo91] with constant uD = usat, is defined here as 5 1wn = 3Twn7n (213) wan =3 where un is spatially varying. In the first term of (212), Ec and Tn are expressed as (Ec, i + E, j)/2 and JT, T, respectively, on each triangle edge [Lia94]. In FLOODS, the velocity v, defined as Jn/qn, is commonly determined by solving the momentum balance equation. The incorporation of the electron velocity into the energybalance equation greatly increases the nonlinearity of the entire system, and thus the computational time during 21 j=3 j=2 j=4 j=I j=5 j=6 Fig. 2.2 A twodimensional structure with triangular meshes. 22 simulation. To reduce this computational burden, Slotboom assumed a constant velocity throughout the device and applied the model to predict the energydependent impactionization rate along a single current path in such devices as BJT and MOSFET. However the model does not provide detailed information about the twodimensional distribution of the electron energy within the devices. Recently, Agostinelli and Tasch [Ago94] demonstrated an energydependent twodimensional substrate current model for the simulation of submicron MOSFETs in which the simplified onedimensional energy balance equation proposed by Cook [Coo82] was applied to many current contours in order to generate a twodimensional representation of electron energy, and thus impactionization rate. However, as we noted in the introduction, these models are only available as postprocessing in device simulation, and thus exaggerate the electron energy, and hence overly predict the impactionization rates. In this section, we present an energybalance model with reliable accuracy, still maintaining the computational efficiency like the model proposed by Slotboom. Before we present the model, some numerical methods are mentioned to facilitate our model development. There are two basic numerical approaches to solve the discretized semiconductor transport equations [Ban86]. One is Gummel's method in which the equations are decoupled and solved sequentially, and the other is Newton's method in which the equations are solved simultaneously, thus accounting for all of the coupling between state variables. The Newton algorithm is very stable and the solution time is nearly independent of the bias condition, even if highlevel injection obtains. For the HD model, it is very difficult to directly implement the full Newton method because of the large Jacobian size and the nonlinearity introduced by the energybalance equation. In MEDICI [Tma92], a decoupled approach for the 23 solution method of the HD model has been chosen. The algorithm consists of a solution of the DD equations for potential W, electron concentration n, and hole concentration p followed by a solution of the energybalance equations for the carrier energies. The carrier energies are used to reevaluate I, n, and p. This procedure is repeated until the system converges. However due to the decoupled nature between the energybalance equation and DD equation, convergence of the HD simulation in MEDICI can be slow for highly scaled devices. Alternatively, FLOODS has chosen the full Newtonlike approach with the DD simulation as a preprocessor. From the DD simulation, the solver defines the initial guesses for variables of the HD simulation. This approach has been verified as a reliable method to characterize the performance of deepsubmicron devices [Lia94], albeit less computationally efficient than MEDICI. Our new SEB model is developed based on the numerical HD method used in FLOODS. The key assumption is that the coupling between the carrier energy and the electric field is more significant than the one between the carrier energy and carrier velocity, which is reasonable because the carrier velocity tends to saturate with the carrier energy. The algorithm of the DC model is flowcharted in Fig. 2.3. The coupling between the electric field and carrier energies, which is strong in scaled devices, is faithfully accounted for, in contrast to postprocessing analytical methods for estimating carrier temperatures, e.g., the one proposed by Slotboom [Slo91]. According to Slotboom's model, the electron energyrelaxation length is assumed to be spatially constant, empirically set at 65 nm. However, we found that even though this assumption has some merit for predicting the peak impactionization rate in simple onedimensional devices, it is not physical and it is not sufficiently accurate for use in general multidimensional device simulation. 24 Set Bias DD Simulation (W, n,p) Converged? No Yes Calculate wn, which define on each edge SHD Simulation; Poisson's Eq. (W) Carrier Continuity Eqs. (n, p) SEB Eqs. (Tn Tp) No Converged? Yes Set New Bias Fig. 2.3 The flow chart of the SEB/SHD model for a specific bias. 25 Alternatively in our SEB model, the energyrelaxation length in (213) is estimated from a preprocess driftdiffusion simulation using the carrier velocity distribution predicted throughout the device domain, as can be seen in Fig. 2.3, and is used without change in a subsequent simpler hydrodynamic (SHD) simulation where Poisson's equation, the carrier continuity equations, and the carrier SEBs, given in (211) for electrons, are solved simultaneously. The carrier energies are thus directly coupled to and defined by the electric field distribution; the predefined energyrelaxation length distribution in the SEB equation is representative, yet eliminates dependences on the carrier densities, thereby substantially enhancing the linearity of the system of SHD equations. The SEB/SHD model is reasonably accurate, and it can substantively increase computational efficiency because of the enhanced linearity. The SEB/SHD model has been implemented as a pragmatic option in FLOODS without substantial code modification due to the high portability of the objectoriented device simulator [Lia94a]. The SHD algorithm and the discretization scheme used are the same as those implemented in the HD model of FLOODS, but the original energybalance equation is replaced by the SEB equation, i.e., (212) with Xwn in (213) evaluated from the preceding DD simulation. The Jacobian matrix elements are calculated accordingly. Note that the Jacobian for SHD simulation is much simpler than that of the conventional HD simulation due to the loosely coupled nature of the carrier energies and densities. All physical models such as energy bands, SRH and Auger recombinations, impactionization rates, and energyrelaxation times have been used without change. The energy dependence of the carrier mobility is incorporated through Klaassen's model [Kla92]. 26 2.3.2 HotCarrier Effects In deepsubmicron MOSFETs, hot carriers generated by the very large electric field cause severe reliability problems and present major constraints on device design, even under normal operating conditions [Ago94]. The hot carriers can be injected into the oxide and trapped, which leads to serious threshold voltage shift and transconductance degradation [Hwa89]. The substrate current measured in a MOSFET provides a proof of the presence of such hot carriers. The SEB/SHD model developed in the previous section can be used efficiently to predict these hotcarrier effects since the effects are fundamentally described by the carrier energy, not by the electric field. In this section prominent nonlocal effects caused by hotcarriers such as velocity overshoot and impactionization are briefly discussed. When the electric field increases rapidly, carrier velocity can overshoot the value corresponding to the local electric field [Jin92]. This enhanced transport occurs because the carrier energy lags the field and remains relatively small. Such a nonlocal effect has long been recognized as significant in scaled MOSFET [Hun77, Sha90]. The SEB/SHD model utility will be assessed through the overshoot analysis in the next section. Many expressions for the impactionization rate exist in the literature [Oga65, Bar91, Sou93]. In FLOODS, the SchollQuade model [Sch87] has been used to describe the process. The model is derived from the BTE under the assumption that the spherical part of the electron distribution function can be represented by a Maxwellian form. The resulting net generation rate per unit volume is given by n (214) Gn = jexp erfc (214) To 24) 0o 7 w ; 27 where z'= kBTT,/wth and to is a characteristic time which depends on the effective mass of the conduction electron and dielectric constant of the semiconductor. In (214), the threshold energy Wth and the time constant to are treated as adjustable parameters according to the measurement. In FLOODS, this carrier generation due to the impact ionization is implemented as a part of the carrier transport equation (22). Alternatively, for computational efficiency, it can be calculated in post processing, which is adequate for the common case of weak impact ionization, in contrast to the case of avalanche breakdown, which necessitates the former exact calculation. The substrate current of a MOSFET is calculated by the integration of the generation rate over the simulation mesh. We stress here that our SEB/SHD model enables simulation of impactionization as well as other hotcarrier effects directly in the coupled numerical iteration process. The SEB/SHD model is not merely decoupled postprocessing. 2.4 Model Verification 4.3.1 Bulk MOSFET The SEB/SHD model option has been used in DC simulation of the 0.5 gm nchannel MOSFET used by Liang [Lia94]. Fig. 2.4 shows the simulated device structure. The gate oxide thickness is 15 nm, the S/D junction depth is 0.3 gm, and the channel doping is 1.4 x 1017 cm3. An LDD with 0.1 gim depth and 1.0 x 1018 cm3 doping density is part of the device structure. Fig. 2.5 illustrates the electric field distributions along the channel simulated with the DD and fullHD models of FLOODS at VDS=5.0V and VGS=5.0V. The electric field and field gradient near the drain predicted by the 28 G S D n+ poly n+ n 0.5m n+ ptype 1.4 x 1017 cm3 Fig. 2.4 The simulated MOSFET structure. The gate oxide thickness is 15 nm, the S/D junction depth is 0.3 gm. 29 350 O DD * HD 6 250 S 150 50 50 0 0.2 0.4 0.6 0.8 1.0 1.2 Distance (pm) Fig. 2.5 Electric field distribution along the channel. The electric field and field gradient simulated by the DD model are larger near the drain. 30 DD model are larger, and the field in the channel is smaller. This is because, as we mentioned before, the velocity overshoot influences the field distribution for the given bias, and defines a higher drain saturation voltage, which translates to the lower electric field near the drain as can be seen the figure. This is a proof that the electron energy evaluated directly from a DD simulation is not correct. Fig. 2.6 shows the electron temperature distribution along the channel simulated with the SEB/SHD model option in FLOODS and compares it with the full HD simulation result. To facilitate the direct comparison between the two models, equal values were used for model parameters like energyrelaxation time (tw = 0.35ps). Note that the strong nonequilibrium behaviors are predicted near the drain end of the channel by both models. The hotcarrier temperature distribution simulated with the SEB/SHD model shows good agreement with the one from the HD model. Fig. 2.7 shows the corresponding electron velocity obtained form the SHD and full HD simulations; the velocity, evaluated as on = Jn/qn, was averaged over the channel depth. The SHD and HD velocity predictions are in close agreement. The velocity distribution obtained from the DD simulation is also included in the figure to emphasize the overshoot that actually occurs in the 0.5gm device. The velocity attained by the electrons is seen to exceed its assumed saturated value of 9.8 x 106 cm/sec near the drain. This overshoot is produced primarily by increase in electron mobility due to Tn lagging E, and becomes more significant as the channel length decreases. Note that the velocity overshoot can be beneficial to the performance of the device since it increases the drain current and reduces the sourcedrain transition time, but perhaps more significantly tends to reduce T, for a given bias condition. Fig. 2.8 shows the twodimensional impactionization generation rate simulated with the SEB/SHD and full HD models at VGS = VDS = 5.0 V. The 31 6000 * HD 5000 O SHD 4000 3000 2000 1000 0.0 0.2 0.4 0.6 0.8 1.0 1.2 Distance ( gm ) Fig. 2.6 Electron temperature distribution along the channel. The location and magnitude of the peak electron temperature show close agreement for both models. 32 2.0 I * HD SHD u 1.5 U DD S1.0 0.5 0.0 0.0 0.2 0.4 0.6 0.8 1.0 1.2 Distance ( tm ) Fig. 2.7 Average electron velocity distribution. Velocity overshoot is observed near the drain. 33 So"' ce 30 U 25, U < 20e 15 0 10 5 0> 1 0.8 0.6 1 ) 0.4 0 0.8 ".2 0.4 (6V. )0 0 0.2 Fig. 2.8 (a) Generation rate by the impactionization process simulated with the HD model of FLOODS. 34 30 C) 25 20. S15 10 , 0.8 0.6 1 S0.8 00..0.2 0.4 0 0 0 0.2 riss Fig. 2.8 (b) Generation rate by the impactionization process simulated with the SEB/SHD model of FLOODS. 35 values used for parameters Wth and ~oI in (214) are 4.0 eV and 7 x 1015 s1, respectively. As seen in (214), the impactionization rate is an increasing function of the carrier energy, and hence its distribution follows the temperature distribution. The SHD and HD predictions agree quite well, thereby verifying the validity of the SHD option. These twodimensional generation rates were calculated via the postprocessing option, which is adequate for the common case of weak impact ionization. As we mentioned in the previous section however, (214), inserted in (22), could be an integral part of the system. The substrate current of the MOSFET is calculated by the integration of the generation rate over the simulation mesh. For SHD and HD simulations, the substrate currents are 3.4x108 A/gm and 3.0x108 A/jm, respectively, further confirming the validity of the SHD option. To demonstrate the real utility of the SHD option, we discuss the computational efficiency of the new model. Since the SHD simulation solves five nonlinear partial differential equations, i.e., Poisson's equation, the electron and hole continuity equations (with the carriervelocity expressions), and the energybalance equations (SEBs) for electrons and holes, its Jacobian has the same order as that of the full HD model. But, as we described in Section 2.3, the SHD model assumes fixed velocity, or energyrelaxation length distributions in the evaluation of the SEB equations. This simplification removes the n and p dependences of carrier energies and simplifies the Jacobian accordingly; the linearity of the system is thus enhanced. This enhanced linearity of the SEB/SHD model relative to the HD model results in reduced numerical iteration counts for a given bias, and hence gives faster run times. Fig. 2.9 shows the average iteration counts for the SHD and full HD DC MOSFET simulations done by sweeping VDS from 0 V to 5.0 V in equal voltage 36 40 H HD ee SHD 30 20 = 20 Ii 10 0.0 0.2 0.4 0.6 0.8 Voltage Step (AVDS) Fig. 2.9 Average iteration count for the SHD and HD MOSFET simulations for each specified biases reached by successively increasing VDS from 0 V to 5.0 V in equal voltage steps at VGS=5.0 V. As the voltage step increases, computation efficiency of the SHD simulation becomes prominent. 37 steps AVDs, for different AVDS. As can be seen in the figure, for AVDS = 0.2V, the SHD and HD simulations need 100 (=25x4) and 125 (=25x5) iterations, respectively, to obtain electron temperature distributions at VDS = 5.0 V; while for AVDS = 0.5V, the SHD and HD simulations need 50 (=10x5) and 150 (=10x15) iterations, respectively. Note that as the simulation voltage step AVDs increases, the relative computational efficiency of the SEB/SHD model increases; that is, for the small voltage step, the simulation time of the SHD model is just 20% shorter than that of the full HD model, whereas for the large voltage step, the SHD simulation time is reduced by a factor of three relative to the HD simulation. This improvement means that the SHD option can enable much faster, yet reliable characterization of scaled devices at biases where the hotcarrier effects are prominent. 2.4.2 Scaled SOI MOSFET Unlike a bulk MOSFET, there is no substrate contact in the thinfilm (commonly floatingbody) SOI MOSFET, and thus majority carriers generated near the drain by impactionization are injected into the film body where they accumulate. This charge modulation forwardbiases the source junction and thus activates the parasitic bipolar transistor [Fos90]. Therefore it is very difficult to characterize the impactionization process explicitly. To clarify this mechanism, more physical simulations including the carrier energy balance equation are required. In this section, SHD simulations are performed to investigate the hotcarrier effect in a deepsubmicron fullydepleted (FD) SOI MOSFET. Fig. 2.10 shows the simulated twodimensional nchannel device structure. The device has a silicon film thickness of 62 nm, a front gate oxide thickness of 10 nm, a back oxide thickness of 0.38 m, and a channel length of 38 G n+ Poly n+ 1.1x1017 cm3 n+ Back Oxide Si tox,=10nm, tsi = 62 nm, tbox=0.38 pm Fig. 2.10 The simulated nchannel SOI MOSFET structure.The S/D junction depth is 0.3 gm, and the channel length is 0.3 uim. 39 0.3 pm. The doping concentration in the body is 1.1 x 1017 cm3 and the doping in the source and drain regions is 5 x 1019 cm3. To focus on the hotcarrier effects, the lattice temperature is set to 300 K, and selfheating is ignored. Fig. 2.11 shows the electron velocity distribution simulated with the SEB/SHD option, and compares it with the DD simulation result. The velocity is defined as the average across the channel. Significant nonlocal effect is reflected near the drain where the attained electron velocity is much higher than that predicted by the DD model. The nonphysical electron velocity spike that commonly occurs in the conventional fullHD simulation is eliminated in the SHD option. Finally, Fig. 2.12 shows the twodimensional carrier generation rate due to the impactionization process, simulated with the SEB/SHD models. Most impact ionization is seen to occur near the drain as we expected. 2.5 Model Extension The SEB/SHD model can be easily extended to threedimensional case, for which excessive, usually prohibitive CPU time is required for solving the HD model due to coupling to the third dimension. In the general threedimensional simulator, the box method, analogous to the twodimensional one, is commonly used to discretize the carrier transport equations, which are integrated over a small volume enclosing each node [Tma92]. The discretized threedimensional SEB equation has the same form as (212) with minor changes; length and area in (212) are replaced by area and volume, respectively. This volume can be divided by several triangular prisms. In this case, the current density is defined along each side of each triangular prism, where then the carrier velocity, and thus the energyrelaxation length distributions can be evaluated accordingly from the preceding DD analysis. This capability of extending to the threedimensional case, potentially with even more 40 2.0 (0 SHD 0 DD 1.5 0 > 1.0 0.5 0.0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 Distance( m ) Fig. 2.11 Average electron velocity distribution for the SOI MOSFET 41 30 o 25 20 15 5 0.06 0.04 0.6 0.5 S 002 0.2 0.3 0 0.0 0 0 0. ace Fig. 2.12 Generation rate by the impactionization process simulated with the SEB/SHD model of FLOODS. 42 enhancement in relative computational efficiency, obviously increases the utility of the SEB/SHD model. So far attention was given to DC, or steadystate device simulation. We next discuss extending the SEB/SHD model for transient simulation, and exploiting the benefit afforded with respect to the time discretization. For the general case, this extension is not straightforward; n and p can be important in the timedependent energybalance equation, which tends to prevent its simplification. However, the extension can be done effectively for special cases. For example, when T,(dn/dt) < n VTn + ) dT = q v E (215) 5 Ir dt 5k n which can be expressed as dT dt= Fn(E, T,). (216) Now the timederivative in (216) must be discretized, considering stability and local truncation error in the selection of the numerical integration scheme. Firstorder backward difference formula (BDF1), secondorder backward difference formula (BDF2), trapezoidal rule (TR), and BDF2TR composite methods have been developed for this purpose [Ban86], and any of these schemes can be used for transient simulation via the SEB/SHD model. For example with BDF1, (216) can be discretized as follows: Tk k1 Atk = Fn(Ek, T) (217) 43 where At = tk tk . The use of (217) at each mesh point follows the spatial discretization of (216) as indicated by (212), and the transient SEB/SHD model is then implemented analogously to the DC case as illustrated in Fig. 2.3. Note that relative to HD, SHD simulation times can be much shorter since larger time steps Atk for transient simulation can be used as implied by the inherent reduced iteration count in Fig. 2.9. Based on the DC simulations we have done with FLOODS, we project an orderofmagnitude improvement in run time for general twodimensional transient MOSFET simulation afforded by the SEB/SHD option; and even more for threedimensional numerical simulation because of the inherent higher degree of nonlinearity in the HD system. 2.6 Summary An efficient multidimensional SEB/SHD model was developed based on assumed fixed energyrelaxation length distributions, derived from the preprocess DD simulation, in the SEB equations for the SHD simulation. The coupling between electric field and carrier energies was properly taken into account, in contrast to postprocessing analytical models and decoupled numerical simulations. However, the coupling between energy and carrier densities, which is typically weak, is implicitly neglected, thereby enhancing the linearity of the system. The new twodimensional SEB/SHD DC model was implemented in FLOODS as a pragmatic option without substantial sourcecode modification. All physical models such as impactionization and mobility are used without change. Model verification was accomplished through submicron bulkSi and SOI MOSFETs simulations. The SHD simulation provides reliable detailed 44 twodimensional representations for carrier temperature, carrier velocity, and impactionization rate. The accuracy of the model is comparable to the full HD model in FLOODS. The most noteworthy feature of the new model is in its computational efficiency, which is attributed to the enhanced linearity of the model transport equations. Substantively reduced iteration counts for larger voltage steps, in DC simulation for example, underlies the enhanced efficiency, which in fact can enable pragmatic threedimensional device simulation. The use of the SEB/SHD model in transient simulation, where quite significant improvements in computational efficiency can result due to large time steps enabled by the enhanced linearity, is not straightforward generally, but follows analogously for special cases. In such cases, and for general DC simulation, the new SEB/SHD option can be easily implemented in any robust two or threedimensional numerical HD simulator. CHAPTER 3 IMPLEMENTATION OF TEMPERATURE DEPENDENCE IN MMSPICE HBT/BJT MODEL 3.1 Introduction MMSPICE was first developed to incorporate the novel concept of seminumerical mixedmode device/circuit simulation [Jeo89], and was verified, through subsequent enhancement [Jin92, Hon94], to be a predictive and reliable IC TCAD tool. Since the chargebased BJT model in MMSPICE is physical and has parameters that depend directly on both device structure and technology, the model capability can be extended straightforwardly to account for new phenomena. The main enhancement of MMSPICE3.0 [Hon94] is a physical but compact model for the SiGebase heterojunction bipolar transistor (HBT), which has been widely investigated in the past few years [Cra93, Gha91, Gib88] as a viable technology due to excellent currentdriving capability and highspeed performance. Although this new version of MMSPICE simulated HBT/BJT devices and circuits successfully at room temperature, it was not applicable to temperaturedependent simulation. Recent interest in lowtemperature operation of HBT/BJTs, however, has emphasized the need for an expanded version of the simulator capable of predicting transistor behavior over a wide temperature range. This is the motivation of this chapter. Bipolar transistors are generally not considered for operation at low temperature (e.g., 77K) because of reduced current gain, increased base resistance, and degraded transit time [Dum81]. The dominant reason for the 45 46 current gain reduction is larger bandgap narrowing in the heavily doped emitter than in the base. In thinbase Si transistors, the base current is typically limited by hole injection into the emitter, and therefore its temperature dependence is largely determined by the effective bandgap of the emitter, whereas the temperature dependence of the collector current is mainly determined by the bandgap in the base. The difference in bandgap between emitter and base causes the current gain to degrade exponentially with decreasing temperature. In addition to reduced current gain, Si BJTs suffer from reduced base conductivity at low temperature because of carrier freezeout which degrades switching performance. The recent advances of bipolar technology such as sustained device scaling [Com91], lowtemperature epitaxialbase process [Gan90], and SiGebased HBT [Pat89], change this situation dramatically. For the advanced Si bipolar transistors, the base doping level must be high to prevent punchthrough as the base width is thinned to minimize the base transit time. The high base doping level combined with polysilicon emitter contact leads to improvement in current gain at low temperature [Woo88]. Also, by incorporating Ge into the base, further improvement in current gain can be obtained since Ge introduces additional bandgap narrowing in the base. In fact, for the recently developed SiGe HBT having larger bandgap narrowing in the base than in the emitter, the current gain actually increases as temperature decreases. Furthermore, carrier freezeout in the base is diminished by using a thin base heavily doped above the Mott transition. This is more prominent in epitaxialbase technology which is capable of achieving significantly thinner base width than an ionimplanted base process. This chapter describes a detailed investigation of the temperature dependence for these advanced Si and SiGebase HBT/BJTs, modeling of the 47 temperature dependence, and implementation of the modeling into MMSPICE3.0. Although most concerns are given to the lowtemperature operation, the device operation above room temperature is also examined for hightemperature applications. For verification and physical insight, MMSPICE simulations in temperature are performed and compared with measurements based on recently developed HBT/BJT technologies. 3.2 Si and SiGebase BJT Characteristics in Temperature Prior to defining the temperature dependence for the MMSPICE HBT/ BJT model, which is facilitated by the physical nature of the model parameters [Hon94], it is useful to discuss the general temperature dependences of familiar transistor characteristics. This discussion will provide good perspective for the MMSPICE implementation described later. 3.2.1 Collector Current The integral relation for the current flow through the quasineutral base region of a homojunction bipolar transistor was first derived by Moll and Ross [Mol56]. In the case of a SiGe heterojunction bipolar transistor with a nonuniform energy bandgap in the base region, the MollRoss current density expression can be generalized in the forwardactive mode as follows [Kro85]: q qVBE/kT Jco = e (31) P2D dx Dnnib 48 2 where nib, the effective intrinsic carrier concentration in the base, is generally positiondependent and given by 2 2 AEgb, HD + AEgb, Ge nib = nioexp kT (32) where AEgb, HD is the apparent bandgap narrowing associated with the heavily doped base region [Rei79] and AEgb, Ge is the additional bandgap reduction caused by germanium. Substituting (32) into (31), we can obtain an analytic expression for the collector saturation current density under low injection for a device with a linearly graded Ge profile across the base: q 2n l'exp( '+ AEgb(0) Jco = (33) WBNAO expil' 1 = Jco'(T) f (', T) (34) where Dn is the average electron diffusion coefficient in the base, WB is the quasineutral base width, NAO is the ionized peak dopant concentration, and 4' = l1+AEgb(grade)/kT, where AEgb(grade) = AEgb(WB) AEgb(O) and AEgb(X) = AEgb, HD(X) + AEgb, Ge(X). A onedimensional SiGebase HBT structure is assumed in the derivation of (33), in which the doping profile can be represented with an exponential function, NA(x) = NAO exp(Tix/WBM) where WBM is the metallurgical base width. In (33), Jco is dependent on temperature both explicitly and implicitly through the several physical parameters. The most dominant term determining the temperature dependence of 2 the collector current is nio. As will be seen in (320), it is an exponential 49 function of temperature, and is divided by four for every 80C decrease near room temperature. This Jco degrades significantly with cooling. Also, the freezeout of majority carriers plays an important role in determining the temperature dependence of the collector current at sufficiently low temperature. This effect can be understood qualitatively as follows. The number of ionized acceptors in the base depends exponentially on the difference between the Fermi level and acceptor ionization energy. Because the Fermi level moves closer to the band edge as the temperature is reduced, it becomes increasing favorable for the acceptors to remain neutral rather than become ionized. This neutralization of the acceptors removes holes from the valence band; the integral in (31), or the effective Gummel number in the base, is reduced so that the collector current increases. Further the reduced hole concentration lowers the base conductivity. Another important parameter related to the temperature dependence of the collector current is Dn in (33), which is defined by the product of the thermal voltage and mobility under low injection. It is slightly reduced with decreasing temperature since the electron mobility does not increase fast enough to overcome the decrease of thermal voltage. We will investigate the temperature dependences of these parameters in more detail in section 3.3. Let us next examine the temperature dependence of f(7', T) in (33), which represents the multiplication factor of the collector current corresponding to both the bandgap narrowing and builtin field: For the advanced Si BJT with an ionimplanted base doping profile, rj' becomes ir AEgb(O)/kT since AEgb(WB) is much smaller than AEgb(O). Thus f(71', T) can be simplified to f(r', T) = exp(n). (35) expl' 1 50 Although (35) is exact, it is more insightful to consider what this equation becomes in two extreme limiting cases. If T >> AEgb(0)/lk, ir' reduces to r, which means that f(Tl', T) is independent of the bandgap narrowing in the base. In such a case the collector current is only dependent on temperature by Jco'(T). Whereas if T < f(Tj', T) = 'exp (A (O (36) Since AEgb(grade) has a positive value for a steep Ge gradient across the base, AEgb(O) is reduced as the emitterbase bias rises, which significantly decreases the collector current. This nonideal behavior of the collector current is physical and verified in the literature [Cra93]. Note that the transconductance, which can be approximated by qIc/kT, increases for the same collector current with cooling due to the reduction in the thermal voltage. This improved transconductance at low temperature is of great practical importance because it allows more efficient driving of capacitive loads [Cre90]. 51 3.2.2 Base Current The collector current given in (33) to first order depends only on the base profile. Unfortunately, the same cannot be said of base current, which depends on emitter contact and emitterbase junction leakage as well as emitter profile. So, temperature dependence of the base current is more complicated than that of the collector current. Under normal emitterbase bias, the predominant mechanism defining the base current is the injection of holes from the base into the emitter which can be expressed as [Fos81] 2 qno 1 WH u t (Age qVBE/kT J = + 1 + exp (37) ND Sp DpoJ TA) kT where Sp is the surface recombination velocity, Dpo is the hole diffusivity in the emitter, TA is the Auger recombination lifetime, and AEge is the effective bandgap reduction in the emitter due to heavy doping. The saturation current 2 is dependent on temperature mainly through both nio and AEge at all temperatures. The nonideal base current at room temperature under lowforward bias is in general caused by the ShockleyReadHall (SRH) recombination within the depletion layer through bandgap traps. As will be discussed in section 3.3, its temperature dependence is largely determined by the intrinsic carrier concentration. However, as the temperature is reduced, the base current has a nonideality greater than two expected in SRH recombination. Various mechanisms such as fieldassisted PooleFrenkel SRH recombination, bandtoband tunneling, and trapassisted tunneling are responsible for this large nonideal factor. The PooleFrenkel effect, which is described by the enhancement of the 52 carrier emission rate due to the large junction electric field, can be written as [Hac88] Iqp1/2E 1/2 IPF = ISEO exp/ 2E1/2 exp (qVBE 1 (38) where E is the baseemitter junction electric field and 3 is q/ni. Since E is a positive function of Obi VBE, IPF has an exponential dependence on VBE weaker than exp(qVBE/2kT), which implies a nonideality larger than two. Emitterbase tunneling is a major limitation to sustained profile scaling in bipolar devices. This tunneling constraint, though problematic for roomtemperature designs, is much more severe for LNT designs since tunneling current is only weakly temperaturedependent and thus does not decrease as rapidly with cooling as diffusion or recombination currents. Since silicon is an indirect bandgap material, any observed tunneling current under forward bias is either phononassisted or trapassisted (excess current) [Hac86]. Recent studies on heavily doped junctions with doping levels between 1018 and 1019 cm3 showed that the excess tunneling dominates the base current at low forward bias[Li88]. 3.2.3 Current Gain Based on the analyses in section 3.2.1 and 3.2.2, the temperature dependence of the current gain is mainly determined by the difference of the bandgap narrowing in the emitter and base: P(T) o exp( kT ge (39) The current gain degrades exponentially with decreasing temperature for conventional homojunction BJTs since the emitter is far more heavily doped 53 than the base, and thus the exponent in (39) is negative. For the advanced Si bipolar transistors, the base doping level should be large as the base width is thinned to minimize the base transit time. The high base doping level combined with polysilicon emitter contact leads to the improvement in current gain at low temperature. Also, by incorporating Ge into the base, further improvement in current gain can be obtained since Ge introduces additional bandgap narrowing in the base. In fact, for the recently developed SiGe HBT having larger bandgap narrowing in the base than in the emitter, the current gain increases as temperature decreases [Cra93]. 3.2.4. Output Conductance To design a highperformance device for analog circuit applications, high Early voltage is most desirable since it allows large smallsignal output impedance. The Early voltage VA is formulated directly from the variation of collector current with collectorbase bias, based on the generalized MollRoss relation [Kro851: DIc p(WB) w B p dx] dWB IC VcB D(WB)nb(WB)o Dnnib dVc VA (310) which can be simplified to VA = qDn(WB)n(WB)fW dx. (311) CCB 0 Dnnib Substituting (32) into (311), a closedform expression for the Early voltage can be obtained: 54 qWBNAO 1 exp(q') (AEgb(grade) VA CCB r' K T (312) In the deviation of (312), a constant diffusivity in the base is assumed. As can be seen in the equation, Ge grading in the base increases AEgb(grade) and thus enhances the Early voltage. The enhancement is more prominent at low temperature since it is dependent exponentially on l/kT. Also, since the collectorbase junction depletion capacitance CCB is reduced with decreasing temperature, VA can be further improved. 3.2.5 Base Transit Time The electron transit time through the base is a key parameter in determining the frequency response of the HBT/BJT. It can be written, based on the MollRoss extension, as [Kro85] S= dy dz (313) =o fP [z Dnnib J Substituting (32) into (313), we can obtain a closedform expression for a device with a triangular Ge profile: TB W 1 (1e)). (314) While this equation is valid under the assumption of low injection and Boltzmann statistics, interesting physical insight can be gained by comparing the difference between Si and SiGe devices at various temperatures. Let us first investigate the Si transistor with an ionimplanted base doping profile. If AEgb(0) < 55 through Dn, that is, 1 /TBDn is constant as illustrated by the curve 'c' in Fig. 3.1. The plots in the figure are normalized at 300K. Even though the minoritycarrier mobility increases with cooling due to a reduction in phonon scattering, the diffusivity degrades at low temperature since the mobility increase is not large enough to offset the decrease of the thermal voltage. Therefore, the base transit time increases with cooling in this case. Whereas if AEgb(0)>>ikT, B can be approximated as 2 AEgb(O) WB kT kT kT DBDn AEgb(O) AEgb(O)e In this case, 1/TBD decreases with cooling for the temperature ranges of interest due to the retarding fields induced from nonuniform bandgap narrowing, as can be seen in curve 'd' of Fig. 3.1. In a Gebased bipolar transistor, 'r' is much larger than one, and thus TB can be simplified to TB = 1 (316) In such a case, the base transit time decreases with cooling because the drift field induced by the Ge grading can compensate for the degradation of the diffusivity in the base. This case is illustrated in Fig. 3.1 by curves 'a' and 'b'. 3.2.6 Emitter Transit Time The temperature dependence of the emitter transit time is defined by the bandgap narrowing difference between emitter and base [Suz91]: TE = TEO exp (317) 56 3 2.5 2 a 1 \b TBDn 1.5 [cm2 1 c 0.5 0 100 150 200 250 300 T[K] Fig. 3.1 Calculated 1/tBDn as a function of temperature with various graded bandgap reductions due to both Ge incorporation and heavily doped impurity. Curves a, b, c and d represent the cases AEgb(grade) = 100, 50, 0, and 50 meV, respectively. This curves are calculated for Tl=3. 57 where tEO is a weak function of temperature involving the emitter junction depth, the polysilicon thickness, and the base width. In a conventional BJT, the emitter transit time in (317) increases drastically with cooling due to the bandgap narrowing in the heavily doped emitter. In a highly graded SiGebase HBT, the emitter transit may be reduced because of larger bandgap reduction in the base than in the emitter and this reduction is most significant at low temperature. This strong temperature dependence in TE may be important for future profile scaling, since in ultrathin base bipolar transistors the emitter transit time is comparable to the base transit time. 3.2.7 Cutoff Frequency The cutoff frequency fT is defined as the frequency at which the commonemitter shortcircuit current gain is unity. This is related to the physical structure of the transistor through the total forward transit time TF by fT = 1/(2n7TF) where loosely 'F can be expressed as the summation of the various delay components: F = EB + E + B + C + BC, scl (318) where TEB and 'BC,scl are the emitterbase and basecollector depletion layer charging time, respectively, and TC is the quasineutral collector time constant, which is given by rCC = RC CBC. The main delay component at low collector current is 'EB, which can be written as kT TEB = C (CEB + CBC) (319) where CEB and CBC are the emitterbase and basecollector depletionlayer capacitances, respectively. These capacitances decrease monotonically from 58 300K to 77K by roughly 20% due to the increases in the builtin voltage for both Si and SiGe devices. Also, as can be seen in (319), TEB is proportional to the temperature and thus favorably influenced by cooling. The last delay term in (319) is the basecollector depletionlayer transit time, tB ,scl = WBC, scl/2vd where WBC,scl is the basecollector space charge layer width and ud is the saturation drift velocity in this layer, which increases by roughly 30% between 300K and 77K. From the analysis of the temperature dependence of each delay components in (319), we can estimate the temperature dependence of the cutoff frequency. The peak cutoff frequency is in general dominated at all temperatures by the base transit time tB. As stated in section 3.2.5, its temperature dependence is complicated and dependent on the technology used. 3.3 MMSPICE Implementation The HBT/BJT model in MMSPICE3.0 [Hon94] is amenable to straightforward implementation of temperature dependence because of its physical nature. Key model parameters are physical and have well characterized temperature dependences: Table 3.1 lists the parameters, which underlie the coupled physicsbased regional analyses that constitute the model [Hon94]. We first discuss their temperature dependences, and then describe the implementation of these dependences in the MMSPICE. 3.3.1 Temperature Dependences of HBT/BJT Model Parameters a) Intrinsic Carrier Concentration The product of the equilibrium hole and electron concentrations in nondegenerate semiconductors is a temperaturedependent constant equal to the square of the intrinsic carrier concentration: 59 Table 3.1 MMSPICE3.0 HBT/BJT Model Parameters Name Description Units Typical UNEPI Electron mobility in epi collector cm2*V1*s1 1.0e3 NEPI Epi doping density cm3 1.0e16 WEPI Epicollector width m 5.0e7 WBM Metallurgical base width m 2.0e7 NAO Extrapolated peak base doping cm3 1.0e18 ETA Base doping gradient factor 3.0 AI Preexponential coeff. for impact ionization rate cm1 7.03e5 BI Exponential coeff. for impact ionization rate VOcm1 1.23e6 VS Electron saturated drift velocity cm*s1 1.0e7 JEO Emitter saturation current density Aom2 1.0e8 JSEO BE SCR saturation current density A*m2 1.0e4 NEB BE SCR emission coefficient 2.0 WSEO Zero bias BE SCR width m 5.0e8 PE BE junction potential barrier V 1.0 ME BE junction exponential factor 0.4 PC BC junction potential barrier V 0.8 DNB Electron diffusivity in base for low injection cm2es1 10.0 PS CS junction potential barrier V 0.6 MS CS junction exponential factor 0.4 CJS Zerobias CS junction capacitance Fom2 1.0e4 CRBI Intrinsic base resistance coefficient V*s*cm2 2.0e4 TB Carrier lifetime in base s 1.0e7 TC Carrier lifetime in collector s 1.0e6 CIF Forward current coupling coefficient 1.0 CIR Reverse current coupling coefficient 1.0 RC Extrinsic collector resistance Q 50.0 RB Extrinsic base resistance Q 100.0 RE Extrinsic emitter resistance 0 10.0 FB Basecharge partition factor 0.5 FC Collectorcharge partition factor 0.5 TE Emitter hole transit time s 1.0e9 UPBASE Hole mobility in intrinsic base cm2*V1*s"1 230. NEBP Peripheral BE SCR emission coefficient 2.0 JEOP Peripheral BE SCR saturation current density A/m 2.0e11 LBE Width of extrinsic base region m 1.5e6 JCOP Substrate saturation current density A/m 1.5e15 DNBH Electron diffusivity in base for high injection cm29s1 10.0 DEGE Bandgap reduction in base at emitter junction eV 0.0 DEGC Bandgap reduction in base at collector junction eV 0.0 WEM Emitter junction depth m 0.0 DEGA Average bandgap reduction in emitter eV 0.0 60 2 = (27kT 3 3/2 3/2 Eg/kT noPo = nio =4 h2 . ) mnmPe (320) which depends on energy bandgap and effective masses for electrons and holes as well as temperature. If we ignore the relatively weak temperature dependences of the effective masses, the temperature dependence of nio in the HBT/ BJT model of MMSPICE is determined by 3 nio(T2) (T2 2 Eg(T2) Eg(T,) (321) ni(T) T1 2kT2 + 2kT, where T1 and T2 are considered as follows. If the nominal temperature TNOM in the .OPTIONS card of a MMSPICE input file is not specified, then T1=TNOM=TREF=270C. Ifa new TNOMnew is specified, then Tl=TNOM=270C and T2=TNOMnew. These two cases are valid if only one temperature request is made. If more than one temperature is requested using .TEMP card, then T1=TNOM=270C and T2 is the temperature specified in the .TEMP card. Since the energy bandgap for silicon increases with cooling, a correction of Eg with temperature is necessary for exact calculation of nio, and is given by [Sze81] aT2 Eg(T) = E(O) T+ (322) with a=7.02 x 104, P=1108 and Eg(0)=1.16eV. As the doping concentration increases above 1017 cm3 the silicon band structure changes due to the appearance of bandtails and an impurities band, the merging of the impurity band with the valence or conduction bands, and modification of the density of states [Swi86]. The net effect of these changes is 61 an increase in the equilibrium carrier concentration product n oP = nie over 2 nio which can be formulated by introducing a phenomenological model parameter AEg, HD as can be seen in (32). This apparent bandgap shrinkage parameter is a strong function of doping concentration, but its temperature dependence is negligible. (b) Base Current Components The backinjected hole current is the dominant base current component under medium emitterbase bias condition, and was given by equation (37). This current contains several parameters which depend on temperature. From the equation, the temperature dependence of this emitter saturation current density JEO in MMSPICE can be modeled by JEO(T2) ioT2) 2 (Aga T2 y(T1, T2) nT exp (323) JEO(T1) io 1( ) kT2 1 where y(T 1, T2) represents temperature dependences of hole diffusivity in the emitter and Auger recombination lifetime, and A~ga is the effective bandgap reduction in the emitter which is defined by both the bandgap narrowing due to heavy doping and the potential barrier of the polysilicon emitter contact. The temperature dependence of JEO is largely dominated by the strong dependence on temperature of the intrinsic carrier concentration. Also, since the emitter saturation current density is exponentially dependent on the ratio of AEga to kT, an exact model for AEga is required. However, it is one of the most difficult parameters to evaluate due to its variation with emitter doping profile and emitter contact technology used. In MMSPICE3.0, we introduce a parameter, DEGA, to represent AEga, which is the only additional parameter had to be added for the temperaturedependent simulation. 62 As mentioned in section 3.2.3, since the BE SCR recombination current, the dominant component of the base current at low forward bias, is dependent on temperature through the intrinsic carrier concentration, its temperature dependence can be modeled in MMSPICE by JSEO(T2) ni(T2) = .(324) JSEO(TI) ni(T1) In advanced scaleddown bipolar transistors, there are many nonideal base current components which can be significant under some condition of bias, geometry and temperature. One of the nonideal base currents is the peripheral leakage current IBP which can be written by [Jin92] IBP = PE JEOP .expT E 1 (325) where PE is the length of the baseemitter perimeter, JEOP is the peripheral saturation current density and rIEBP is the peripheral emission coefficient. The variation with temperature of the JEOP is given by [Cha91] JEOP = A. exp() (326) where the activation energy, Ea, is 0.64 eV which is close to half the silicon bandgap energy, and thus IBP is identified by a generationrecombination current along the baseemitter junction perimeter. The ideality factor is found to vary from 1.75 at 140k to 1.45 at 300K. Other nonideal components of the base current such as fieldassisted PooleFrenkel SRH recombination, bandtoband tunneling, and trapassisted tunneling, which can be dominant at low temperature, are not modeled in current version of MMSPICE. 63 (c) Mobilities and Diffusivities Accurate models for carrier mobilities are required if predictive simulation results are to be obtained. In MMSPICE, the models for the mobilities are classified, based on electric field, doping level, and temperature, as follows. The carrier drift velocity is proportional to the electric field when the field intensity is less than a critical value, where the proportionality constant is the mobility, and the velocity is saturated at vs above the critical electric field. At low electric field, the majoritycarrier mobility varies more rapidly with temperature for lower doping concentrations, where the strongly temperaturedependent lattice scattering mechanism is dominant. For very high doping concentration where ionized impurity scattering dominates, the mobility is almost temperatureindependent. Therefore an adequate model for carrier mobility must account for each scattering mechanism by reflecting the variation in mobility with lattice temperature and impurity concentration. The empirical mobility model introduced by Arora [Aro82] is used as a basis in our temperaturedependence modeling. From the model by Arora, the majorityelectron mobility is expressed as (T) 88 57 7.4 x 108 T2.33 n, maj(T) = 88 + (327) nm 30NEPI I T 0.146 1 + NEPI 0.88 1.2 x 1017 (T/300)2.4). 300 where NEPI is the donor, or epicollector concentration. Recent experimental results have shown that, starting at a doping concentration of 1017 cm3, the minoritycarrier mobility in silicon exceeds the majoritycarrier mobility, up to a factor of three at a concentration of 1020 cm3 [Swi86]. So different models are needed for minority and majoritycarrier mobilities. Empirical models for the doping dependence of minorityelectron mobility have been defined [Swi861 64 and extended to incorporate the temperature dependence. From this model, the temperature dependence of the minorityelectron mobility is written as ,mi,(T)= 232T 0.57 6.97 x 108 T2.33 (T) = 232 + (328) 0NA T 0146 1 + .88 9.22 x 10 16 (T/300)2.4 300 where NA is the acceptor, or spatial average of the intrinsic base concentration, and is calculated by NAO exp (ri/2). Fig. 3.2 shows the electron mobility as a function of temperature for various impurity concentrations. Solid lines represent the majorityelectron mobilities, which increase with cooling at lower doping levels due to reduced lattice scattering and become flat for the higher doping concentrations. Dashed lines reflect model calculations for the minorityelectron mobilities which have trends similar to with the majorityelectron mobility, but have almost a factor of two higher values for heavily doped silicon. Similarly, the majorityhole mobility in the intrinsic base region can be expressed as [Aro82] p, mj(T) = 54 T 0.57 1.36 x 108. T2.33329) 1 + 17 2 .88(2.35 x 107 (T/300)2.4 300 Note that the mobility expressions in (327), (328), and (329) are not directly used to define actual mobilities in MMSPICE, but the expressions implicitely reveals the accounting for the temperature dependences of the mobilities. For example, from (327), the temperature dependence of the majorityelectron mobility in the epicollector region is modeled in MMSPICE by ( Tl "0.57T2"2.33 2 0.57 n, maj(T2) = (tn, maj(T) 88(T0. I T 2. (F(TI T2)+ 88 300 T(330300 (330) 65 104 " 3 N = 1016 cm3 N = 10 cm3 Eln 10 cm2Vs ] N =1018 cm3 10 50 100 150 200 250 300 350 400 T[K] Fig. 3.2 Majority (solid) and minority (dashed) electron mobilities as a function of temperature for various impurity concentrations. 66 where I+ NEPI 10.146 1 + .0.88/ 1 1.2 x 1017 (T,/300)24 300 F(T1, T2) = x I ) (331) NEPI .088T2 0.146 1.2 x 1017. (T2/300)24 0 Similar expressions for the temperature dependences of the minorityelectron mobility and the majorityhole mobility in the base can be obtained, based on (3.28) and (3.29), respectively. At high electric field, the drift velocity is not proportional to the electric field any more, but become saturation. Several empirical models for the saturation velocity have been presented. After comparing the previous published models, Selberherr's model was chosen because it provides the best fit with the fewest fitting parameters. Based on the model, the drift saturation velocity for electron can be written by vs = 1.38 1017 tanh( (332) From above equation, the temperature dependence of the saturation velocity in MMSPICE is determined by tanh(175/T2) VS(T2) = Vy(TI). (333) Vtanh(175/T 1) As can be seen in equation (32) and (310), the physical parameter determining the collector current and the base transit time is the electron diffusion coefficient in the base instead of electron mobility. So it is necessary to consider the temperature dependence of the diffusion coefficient. The electron diffusion coefficient is generally defined by 67 Dn = 2kTF(l/2 EkT Fn 1/2(EcEFn (334) _qkT ) Y334 where F is the Fermi integral. For nondegenerate case, the above equation can be reduced to well known Einstein's relation, which is assumed to be valid in our model. (d) Junction Potentials The builtin potential across a emitterbase junction under zero bias can be expressed by OE = kTIn (335) q nibnie where nib and nie are the effective intrinsic carrier concentration in both emitter and base, respectively. Therefore, the effect of temperature on the builtin potential OE in MMSPICE is modeled as follows: T2 k T T 3 I AE 2T2 E(2 OE(T ) In ) Eg(T)Eg(T (336) where AEg is the spatial average of the effective bandgap narrowing in the space charge region of the baseemitter junction. Generally, OE increases with decreasing temperature due to significant reduction of the intrinsic carrier concentration. Similar expressions can be used for both basecollector and collectorsubstrate builtin potentials. (e) Recombination Lifetimes Recombination lifetimes in the base and in the collector can be largely determined by the ShockleyReadHall (SRH) recombination process since the 68 impurity concentrations are not so high enough to be influenced by Auger recombination process. The SRH lifetime for electron is well defined as In = 1/((nVthNt), where on is the capture cross section, Vth is the thermal velocity of the electron, and Nt is the trap density. The temperature dependence of SRH recombination lifetime is generally controlled by the temperature dependence of on and Vth, and the result yields [Hen77] SRH(T) = AT .exp (337) where A is a constant determined by the doping concentration and AE is the activation energy for the minoritycarrier capture process and its value is typically 0.1 eV. Consequently, the temperature dependence of the SRH recombination lifetime can be modeled in MMSPICE as TSRH(T2) = TSRH(TI) exp 1160( (338) (T 0( I (338) 3.3.2 Model Implementation in MMSPICE3.0 Unlike the GummelPoon BJT model for which additional empirical parameters are needed to account for transistor behaviors with temperature, the HBT/BJT model in MMSPICE does not require those parameters to describe the temperaturedependent simulation because of the physical nature of the model. In fact, only one additional parameter DEGA, which is defined in section 3.3.1, is introduced to complete the temperature dependence of the model. To recognize the new parameter through the MMSPICE input file, eight related subroutines are modified: ADDELT, ALTER, FIND, READIN, SPICE, TMPUPD, MODCHK, and QBBJT. 69 The default simulation temperature in MMSPICE is 270C, for which userspecified model parameter must apply. The TNOM option can be used to change the temperature corresponding to the given parameter values. The actual simulation temperature is defined by the .TEMP control card if not defaulted. Temperature appears explicitly and implicitly through physical parameters in the HBT/BJT model equations defined in section 3.3.1. We encoded the temperature dependences, based on the temperature specified on the .TEMP card, of the following HBT/BJT model parameters into TMPUPD: electron mobility in epi (UNEPI), carrier saturated drift velocity (VS), emitter saturation current density (JEO), BE SCR saturation current density (JSEO), BE equilibrium SCR width (WSEO), BE junction potential (PE), BC junction potential (PC), electron diffusivity in base (DNB), substrate junction potential (PS), zerobias collectorsubstrate junction capacitance (CJS), intrinsic base resistance coefficient (CRBI), carrier lifetime in base (TB) and collector (TC), peripheral SCR emission coefficient (NEBP), peripheral saturation current density (JEOP), hole mobility in intrinsic base (UPBASE), and substrate saturation current density (JCOP). 3.4 Verification and Discussion The temperature dependence of the HBT/BJT model in MMSPICE3.0 has been verified based on the most advanced devices including two SiGebase HBTs from IBM and one lesser scaled Sibase BJT from Motorola. The first SiGe device, SiGe1 [Cre93], was fabricated using the lowtemperature epitaxy process which have great potential for achieving the verythin, abrupt base profiles required for continued vertical scaling of silicon bipolar transistors. The device has a base width of 56 nm, a peak base doping concentration of 70 3.7x1018 cm3 and a emitter area of 0.6x4.3 pm2. The Ge is linearly graded across the base with 0% at the emitter junction and 8.4% at the collector junction. The second device, SiGe2 [Cra93], was fabricated in a selfaligned structure with a reduced thermalcycle emitter process using phosphorus. The triangular Ge profile is grades 0% at the emitterjunction to approximately 15% at the collector junction. The base width is 35nm and the intrinsic base sheet resistance 16 KQ/square. The emitter is heavily doped and contacted by polysilicon. The effective emitter area is 0.6x4.4 ptm2. The third device, MOSAIC3 [Zde87], is a selfaligned BJT with polysilicon emitter and oxideisolation. The device has a peak base doping concentration of 1x1018 cm3, a base width of 150 nm, and a collector doping concentration of 2x1016 cm3. The effective emitter area of the simulated device is 1.2x3.2 11m2 Though most of the HBT/BJT model parameters in MMSPICE are physical and can be determined directly from the technology using SUMM [Gre90], a few parameters such as diffusivity and emitter saturation current density cannot be determined accurately due to uncertainties in physical mechanisms and process variations. Hence, some parameter tuning is necessary to enhance the accuracy. Fig. 3.3 shows the measured (dashed lines) and simulated (solid lines) collector currents with temperature for the model parameters based on SiGel. The simulation data are nearly the same as the measurement data with minor parameter tuning at room temperature. Furthermore, as the temperature decreases, the simulations predict, without additional tuning, the temperaturedependent behaviors of the collector current over a wide temperature range quite well, except for extremely low temperature at which the carrier freezeout effect is prominent. Since the current version of MMSPICE does not include the proper model for the carrier freezeout effect, a discrepancy is found between the simulation and measure 71 100 VCB = 0 V 101 WB=56nm NB=3.7x1018cm3 2 Ge(WB)=8.4% 10.3 103/ / / 4 / 104 Ic 105 [A] 6 106 310 10.1 250 7 10 182K 133 10 84 ' 10 1010 0.4 0.6 0.8 1 1.2 VBE [V ] Fig. 3.3 Measured (dashed) and simulated (solid) collector current as a function of emitterbase voltage for SiGel. This figure shows an improved transconductance with decreasing temperature. 72 ment data at 84K. In the real device, there is considerable incomplete ionization of dopant impurities at that temperature, which yields the reduction of the hole density in the base. Consequently the effective Gummel number decreases with cooling and the collector current is larger than expected values in the simulation, as can be seen in Fig. 3.3. However, our temperature model is generally accurate enough to estimate the temperature dependence of the collector current. Note that this favorable result is obtained without introducing any new parameter for the collector current, which shows the predictive and physical nature of MMSPICE. This is clearly different from other simulators like SPICE 2G.6 and HSPICE in which additional parameters have to be introduced for the simulation of the collectorcurrent variation with temperature because of the empirical properties of their model parameters. Fig. 3.3 also reveals that the transconductance, the slope of the collector current versus the emitterbase voltage curve, increases with cooling, as discussed in section 3.2.1. Fig. 3.4 illustrates the base current as a function of emitterbase voltage for various temperatures. Under moderate injection levels, the backinjected base current by hole into the emitter is dominant component at all temperaqVBE/kT tures, which can be verified by its dependence on e This term is also strongly dependent on the total bandgap reduction through the minoritycarrier transport in the emitter, which can be explained as follows. The hole concentration increases due to the bandgap narrowing in heavily doped emitter. Since the total bandgap reduction is very sensitive to the process and the technology, it is impossible to use a hardwired value in the simulator, and thus we conclude that a new parameter, DEGA must be introduced. In fact, DEGA is the only parameter to be added for the temperaturedependent simulation of MMSPICE. Under lower injection levels, SRH recombination in the spacecharge region of the emitterbase junction is a dominant mechanism 73 10 102 WB=56nm VCB = 0 V NB=3.7x010cm3 Ge(WB)=8.4% 104 5 10 IB [A] 106 310K 7 107 250 / 182K 8 10 13 / 19 ,', 84K 10 0.4 0.6 0.8 1 1.2 VBE [V ] Fig. 3.4 Measured (dashed) and simulated (solid base current as a function of emitterbase voltage for SiGe at a variety of temperatures. 74 responsible for the base current near room temperature. However, other parasitic base leakage mechanisms such as tunneling and PF recombination become significant with decreasing temperature. Since the current version of MMSPICE does not include models for these mechanisms, there are some discrepancies at low temperature, as shown in the figure. Under high injection, since the base current is limited by several factors like the lateral spreading of electrons, internal and external resistors, heterojunction barrier effect, it is difficult to predict accurately the temperature dependence of the base current. However, the simulation results are reasonably good. Fig. 3.5 shows output currentvoltage characteristics for SiGe1 simulated with MMSPICE3.0. For this device, since the total bandgap reduction is larger in the emitter than in the base, the collector current is reduced with cooling for the same base current, as the collector current of a conventional BJT is. Therefore, the simulations predict that the commonemitter breakdown BVCEO, which is determined by the condition P(M 1) = 1, increases with cooling due to reduced current gain, as can be seen in the figure. But the measurement [Cre93a] indicates that BVCEO is almost the same, 3.9V, for different temperatures. This can be accounted for by the variation of the multiplication factor (M1) with temperature which is assumed to be constant in our model. In fact, at low temperature, the electron impact ionization rate increases due to reduced phonon scattering, which causes the increase in M1. This increase in the multiplication factor with cooling offsets the reduction of P and yields temperatureinsensitive BVCEO. Therefore, to predict the variation of the commonemitter breakdown voltage with temperature, a accurate temperature model for M1 is necessary, which is left for future work. Also, the output conductance decreases with cooling for this device as can be seen in the figure, which is very useful in future lowtemperature analog applications. 75 x 103 0.9 WB=56nm NB=3.7x1018cm3 0.8 Ge(WB)=8.4% 0.7 0.6 C 0.5 [A] 0.4 182K 0.3 85K 0.2 0.1 IB = 10 pA 01 2 3 4 VCE [V ] Fig. 3.5 Simulated IC versus VCE characteristics for SiGe] with IB = 10 gLA. The Early voltage, and hence output impedance, increases with cooling. 76 As discussed in section 3.2.7, the total transit time, hence cutoff frequency, is dependent on temperature through the combination of several delay components. At low collector current region, fT is mainly limited by the emitterbase depletion layer charging time which is favorable with decreasing temperature due to the reduction of both the capacitance and the thermal voltage. Therefore, fT increases with cooling at that region, as shown in Fig. 3.6. Also, the peak cutoff frequency, which is generally determined by the base transit time, rises slightly with decreasing temperature due to reduced base transit time as described in section 3.2.5. Fig. 3.7 shows the simulated current gain versus the collector current at various temperatures for SiGe2, which was fabricated using a more aggressive technology. That is, the device has a much shorter base width and much larger germanium dose compared with SiGel. As can be seen in the figure, the peak current gain rises with decreasing temperature since the bandgap narrowing induced by both germanium and dopant impurities in the base overcome the total bandgap reduction in the emitter. This is in clear contrast to the behaviors in temperature of a conventional BJT or SiGe1. So far a verification has focused on lowtemperature operation for SiGebase HBTs, partly because of the growing interest in the lowtemperature operation of SiGebase HBTs. However, for the conventional BJTs, it is necessary to look at the transistor behavior near and above room temperature where these devices typically operate. Fig. 3.8 shows the measured and simulated current gain for MOSAIC3 at and above room temperature. As expected from the analysis in section 3.2.3 and verified with measurement data, the simulated current gain rises with increasing temperature. Under high injection levels, however, some discrepancies are seen, which are more severe at higher temperatures. The limited verification is further demonstrated by the 77 70 WB=56nm VCB = 1 V 60 NB=3.7x108cm3 Ge(WB)=8.4% , 50 85K / 40 309K [GHz ] 30 20 10 0.01 2 100 10 10 Ic [A] Fig. 3.6 Measured (dashed) and simulated (solid) transistor cutoff frequency as a function of collector current at 309K and 85K for SiGe'. 78 800 WB=35nm NB=4.0xl018cm3 VCB = 0 V Ge(WB)=15% 63K 600 192K 234K S400 300K 200 08 108 106 104 102 Ic [A] Fig. 3.7 Simulated current gain versus collector current with various temperatures for SiGe2 at VCB = 0 V. This figure shows that the current gain rises with cooling due to larger effective bandgap narrowing in the base than in the emitter. 79 AT = 25 K WB=150nm VCB = 1 V NB=1.0x1018cm3 1.5 398 K 298 K 0.5 10.5 104 103 102 Ic [A] Fig. 3.8 Measured (dashed lines) and simulated (solid lines) forward current gain versus collector current over temperatures for MOSAIC3. 80 emitterbase voltage versus temperature curves for various collector currents illustrated in Fig. 3.9. We believe the highVBE (highcurrent) discrepancies are related to the temperature dependences of both internal and external series resistances as well as diffusivity under highinjection conditions. 3.5 Summary A detailed analysis of temperature dependences of Si and SiGebase HBT/BJTs was presented and implemented into the existing bipolar transistor model in MMSPICE3.0 with one additional parameter. It was shown that key device characteristics have different temperature dependences according to bandgap narrowing and doping gradient in the base. For example, the peak current gain increases with decreasing temperature when the bandgap narrowing induced by both germanium and dopant impurities in the base overcomes the total bandgap reduction in the emitter. The cutoff frequency, which is generally determined by the base transit time, rises slightly with decreasing temperature due to reduced base transit time. This is in clear contrast to the behaviors in temperature of a conventional BJT. The MMSPICEsimulated and measured DC and AC characteristics are in good agreement for both the SiGe HBT and Si BJT, over a wide temperature range, which proves the practical utility of the temperature model. 81 1.1 WB=150nm NB=1.0x1018cm3 0.9 IE= 6 mA 0.8 VBE 0.7 1 mA [V] 0.6 0.1 mA lO pA 0.5 1pA 0.4 0.3 300 320 340 360 380 400 420 T [K] Fig. 3.9 Measured (dashed) and simulated (solid) emitterbase voltage versus temperature with various emitter currents for MOSAIC3. CHAPTER 4 PHYSICAL SOI MOSFET CHARGE MODEL FOR LOWVOLTAGE ANALOG AND MIXEDSIGNAL CIRCUIT APPLICATIONS 4.1 Introduction MOSFETs built on SOI have yielded superior performance in high frequency analog and mixedsignal applications mainly because of the inherent device isolation afforded by the underlying oxide, which reduces parasitic source/drain capacitance and allows use of highresistivity substrates for lowloss transmission lines [Co196] [Hur96]. The thinfilm nature of the SOI MOSFET however can underlie mechanisms that complicate circuit simulation and portend equivocation in design [Fos97]. Mainly because of floatingbody effects, empirical compact models like those used for designing bulkSi CMOS circuits are not adequate for reliable SOI circuit simulation. SOISPICE has been developed mainly for digital applications and verified as a reliable simulation tool for lowvoltage and highspeed digital circuit design [Fos97]. Its physical, processbased models can be effectively used to characterize the fully depleted (FD) and nonfully depleted (NFD) SOI MOSFETs. However to apply the models to analog circuit design, model enhancement is necessary because, unlike digital circuits, where performance depends on relatively few parameters such as threshold voltage and drive current, the performance of the analog circuit is quite sensitive to the smallsignal parameters as well as largesignal characteristics of the SOI MOSFET. Therefore the SOI analog model should not only meet common requirements 82 83 for digital design such as reasonable IV characteristics, delay prediction, and charge conservation, but also should give accurate values for all smallsignal quantities such as transconductance and transcapacitances; and all of these parameters should be continuous with respect to any terminal voltage. Furthermore, floatingbody effects, which introduce additional difficulties in SOI analog modeling due to frequency dependences, must be properly accounted for. In this chapter, we present scalable charge models of the FD and NFD SOI MOSFETs for analog applications of SOISPICE (Ver. 4.4) focusing on the continuity of capacitances. The modeling, implemented directly without additional parameters, enables reasonable AC simulation and reduces the simulation runtime relative to that of previous versions of the SOISPICE models. A unified floatingbody model is described based on the charge models, which predicts well the frequencydependent floatingbody effect as well as DC and transient floatingbody effects. 4.2 SOI MOSFET Charge Models As illustrated by the network representation in Fig. 4.1, SOISPICE FD and NFD SOI models are chargebased with five terminals and they allow the floatingbody option. The parasitic bipolar current and charge are intrinsically coupled to the MOS modeling [Kri96]. The carrier generation/recombination currents are properly accounted for. Since the current and charge models have been developed based on the regional approach, the derivatives of drain current and charges are not continuous at some boundaries. Therefore these models can not be used for analog circuit simulation without upgrade. In this section the charge models are expanded and improved so that all of the charges 84 Gf SdQGf ICH dt IBJT dQD RS+RLDS VLDD RD+RLDD SR S ds D D dt VdQs dt iGi IRGt(VBS) B' IRGt(VBD) dQGb R dt RB B Gb Fig. 4.1 Network representation of the chargebased SOI (nchannel) MOSFET models in SOISPICE. 85 and their first derivatives are continuous through all regions of operations. The current model is also modified accordingly. 4.2.1 NFD/SOI The intrinsic quasistatic frontgate charge can be modeled using Gauss's law as QGf = WCofL (VGfS Ofs Vsf (Y))dy (41) where W is the channel width, Cof (= ,ox/toxf) is the frontgate capacitance per unit area, <, is the frontgate workfunction difference, and Wf(Y) is the frontsurface potential. Fig. 4.2 shows MEDICIsimulated potential distributions along the surface at different drain and gate biases for a 0.2 pm NFD/SOI MOSFET. As can be seen in Fig. 4.2, in weak inversion, the VDS dependence of the surface potential is relatively small and restricted to near the drain; the draininduced barrier lowering is hardly noticeable. In such case, (41) can be simplified as QGfW = CofWLe(VGfs ms sf) (42) where Le is an effective channel length, and ysf is the average surface potential calculated by [Suh95] 1 r V Qb EstbVDs = f + B2C + cVYS + EtbVDS (43) oFB Cof(1 + cx)L In (43), a = Cb/Cof where Cb = s/tb, and Qb is the body charge density. The effective channel length, Le=LyDYS, is defined in weak inversion as the portion of the channel where carrier diffusion is predominant and where the 86 3.0 VGfS=O.OV (Weak inversion) 2.5    VGfs=.0OV (Strong inversion) VDS=1.5V 2.0 1.5 1.0 0.5 0.0 0.1 0.2 0.3 0.4 y (im) Fig. 4.2 MEDICIsimulated potential distributions along the surface at different drain and gate biases for a 0.2 gm NFD/SOI MOSFET. 87 longitudinal electric field is small relative to the fields closer to the source and drain junctions [Suh95]; YD and ys account for a weakinversion channellength modulation in the drain and source sides, respectively. To characterize YD, the twodimensional Gauss's law was applied near the drain, yielding a VDSindependent [Gre93, Suh95] yD =icln 2 (44) Estb where l = 2Co + ) This VDsindependence implies yD = ys for any VDS 2Cof(1 + cX) because the symmetry of the MOSFET device means yD = ys for VDS = 0. However we found some inaccuracy in (44) due to an improperly defined boundary condition, AEy(L YD) = 0, used in its derivation. To rectify this problem, instead of rigorously modeling Le, which would add excessive nonlinearity, we decided to modify the original expression for Le only slightly by using an empirical parameter (( 1): Le L 2yD (45) Note this simple fix is justified by MEDICI simulations, which suggest yD depends only weakly on VDS. However, the expression for Le in (45) becomes invalid when VGfS goes to the accumulation region since (44) was derived assuming only depletion charges under the gate, which is reasonable in the weakinversion region but not in the accumulation region. In such case, Le goes back to L and ,sf approaches VBS as the depletion region disappears. Therefore the frontgate charge is defined in accumulation as QGfA = CofWL(VGfs ms VBS) (46) 88 To ensure a continuous transition of QGf between the accumulation and weakinversion regions, a smoothing function is introduced: QGf Q ln(1 +eA(QGfW QGfA)) QGf = Qef A (47) where A = A'/(CofWL) and the coefficient A' is a hardwired parameter that controls the stiffness of the transition. Note the smoothing function guarantees the continuity of the charge and its derivatives as a function of VGfS. In weak inversion, the channel charges are neglected. Thus, aside from the parasitic bipolar components [Kri96], the drain and source charges are only due to bodycharge modulation by yD and ys, respectively, and defined as QD = qtbNBODYWyD (48) Qs = qtbNBODYWYS (49) where tb is the body thickness, and NBODY is the body doping density. The above charge expressions should be modified in the accumulation region in accord with (46), but since the inaccuracies are small, these modifications were not done. Next consider the stronginversion case. Much of the recent analog modeling effort for bulkSi CMOS has been aimed at removing the discontinuity in smallsignal conductance at the transition between the linear and saturation regions [Cha95]. For this purpose, the following smoothing function has been commonly used [Par91] 1 +eB(I VDS/VDSAT)) F(VDS, VDSAT) = 1 (410) In(1 + eB) 89 where VDSAT is the drain saturation voltage and B is a constant which controls the stiffness of the smoothing function. Then an effective drainsource voltage VDSX can be defined as VDSX = F VDSAT (411) which smoothly varies between VDS and VDSAT. The same smoothing function can be applied to the SOI modeling. In such a case, the frontgate charge for strong inversion is expressed as lin sat (412) QGf = GfS + GfS(412) lin sat where QGfS is the gate charge component between y=O to y=LAL and QGfS is the gate charge component in the saturation region from y=LAL to y=L; AL is the modulated channel length in strong inversion [Suh95]. These charges, with (411), are given by lin f L)C of VDSX VDSX(1+s)(1+a ) GfS = W(L L) VGfS m s fS 2 12 j() (413) of Gf WLC S W mssfS VDSX 2eff oI c 1 (414) where s=4effVDSX/(2vsatL), u=Qc(O)/[Cof(l+o)VDsx], Qc(O) is the channel charge at y=O, and 'VsfS is the surface potential in the stronginversion region. The AL can be obtained by solving a differential equation based on the twodimensional Gauss's law near the drain. The solution [Suh95] of the differential equation, with (411) is 90 VDS VDSx 2 sinh + 2 .[cosh 1 (415) 4 eff .lc 7, From (415), AL can be calculated straightforwardly with a given VDSX. Note lin that when AL=0, QGfS = GfS' which is consistent with the conventional definition of the charge. Fig. 43 shows VDsx and (LAL)/L as a function of VDS for a 0.2 lm nchannel NFD MOSFET. As can be seen in the figure, smooth transitions are achieved for both VDSX and the modulated channel length, in contrast to the corresponding results of the previous version of the model (in SOISPICE4.3), also shown in the figure. The channelcharge component of the stronginversion drain and source charges are computed from the original formalism [Suh95], but with smooth VDSx and AL. First total channel charge is calculated by lin sat QCH = CH+ QCH (416) where lin (2z (z 1) (417) QlCH = W(L AL)Co(1 + o)VDSX 2z 1 +(u z) (417) and QsaH = WALQc(L AL) (418) where QC(LAL) is the channel charge at y=LAL. The drain charge lin sat QD( = QD + QD ) is then computed from conventional charge partitioning method [Vee88] 91 1.2 I I I 0.2 gm NFD SOI MOSFET (VGfs=1.5V) 1.0 (L AL)/L 0.8 VDSAT 0.6 0.4 VDSX 0.2 0.0 0.0 0.2 0.4 0.6 0.8 1.0 1.2 VDS (V) Fig. 4.3 Effective drain voltage and normalized channel length for the NFD device. The dashed curves reflect the previous version (in SOISPICE4.3) of the model. 92 lin 42( 1) 4 (z 1)5 (u z) QD =WLCf(+ )VDSX[ 2z 1 15 (2z 1)2 (419) and QD = W AL c(Le). (420) Then the source charge is given by Qs = QCH QD. (421) So far we described charge modeling in the weak and stronginversion regions. Next consider charges in moderate inversion region. Previously, a linear interpolation was used for the charge calculation in this region. However, it shows discontinuities at the boundaries, VTW and VTS, which is not acceptable for analog circuit design, especially for lowvoltage applications where operating bias points might be in moderate inversion. In the new charge model, numerical interpolation as used in BSIM3 [Hua94] is used to describe the moderate inversion charges. For the interpolation, four boundary values ( QTW, QTW', QTS, QTS' ) are needed to define the charges; QTW is the charge at VGfS=VTW and QTW' is the first derivative evaluated there, and QTS and QTS' are the charge and its first derivative atVGfs = VTS. Note here that Q represents QGf, QD, or Qs. The interpolation is efficiently implemented in SOISPICE4.4 so that the additional calculations associated with the new numerical scheme can be minimized. Hence, the charges in moderate inversion can be written as QMOD = (1 t)2QTW + 2(1 T)TQp + T2QTS (422) 93 where (VTw Vp) + J(VTW VP)2 (VTW 2Vp + VTS)(VTW VGfS) (423) VTW 2Vp + VTS The intersection points, Vp and Qp in (422) and (423) and illustrated in Fig. 4.4, are calculated as follows: VP = (VTsQTs' VTWQTw' QTS + QTW) (424) QTS' QTW' and QP = QTW + (VGfs VTW)QTW'. (425) Henceforth, continuity of the charges and transcapacitances at the moderateinversion boundaries is ensured. Gauss's law gives the backgate charge under the intrinsic body region as QGb = WLCob(VGbS ms VBS) (426) where QB = (Qs + QD + QGf + QGb) (427) The terminal charging/discharging currents indicated in Fig. 4.1 are then characterized using the quasistatic approximation: dQ. QidV dt = aV3 dt (428) 