PRAGMATIC AND RELIABLE DEVICE/CIRCUIT SIMULATION FOR DESIGN IN ADVANCED SILICON-BASED TECHNOLOGIES
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
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 Chen-Chi 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. Ghy-Boong Hong and Mr. Ming-Yeh 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, Meng-Hsueh Chiang, and Jonathan Scott Brodsky.
I also thank my parents, father-in-law, Dr. Jong Duk Lee, mother-in-law, 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.
TABLE OF CONTENTS
ACKNOW LEDGM ENTS .................................................... ...................... ii
ABSTRACT ............................................................. vi
1 INTRODUCTION ..................................... ................................... 1
2 SIMPLIFIED ENERGY BALANCE MODEL FOR PRAGMATIC
MULTI-DIMENSIONAL DEVICE SIMULATION ...................... 8
2.1 Introduction .............................................................................. 8
2.2 Hot-Carrier Transport Theory ...................................... 11
2.2.1 Carrier Transport Equations ...................................... 11
2.2.2 Non-Local 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 Hot-Carrier 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 SiGe-base 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
3.3.1 Temperature Dependences of HBT/BJT
Model Parameters..................... ............ 58
3.3.2 Model Implementation in MMSPICE-3.0.................. 68
3.4 Verification and Discussion ................................................. 69
3.5 Sum m ary ..................................... ......................... 80
4 PHYSICAL SOI MOSFET CHARGE MODEL FOR LOWVOLTAGE MIXED-SIGNAL 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 Floating-Body 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 10-bit 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
AN AUGMENTED DRIFT-DIFFUSION MODEL FOR
EFFICIENT NUMERICAL DEVICE SIMULATION ................... 150
REFERENCES ........................................ 165
BIOGRAPHICAL SKETCH ..................................... 172
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 SILICON-BASED TECHNOLOGIES
Chairman: Dr. Jerry G. Fossum
Major Department: Electrical and Computer Engineering
This dissertation describes pragmatic application-specific mixed-mode device/circuit simulation approaches for design in contemporary Si-based 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 energy-relaxation length is estimated from a pre-process drift-diffusion simulation using the carrier-velocity 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 two-dimensional SHD and full-HD DC simulations of submicron bulk-Si 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 semi-numerical approach is considered for mixed-mode device/circuit simulation where physical compact models are developed based on device structures. The approach is first applied to advanced SiGe-base HBTs through MMSPICE to investigate their temperature dependences. MMSPICE-simulated and measured DC and
AC characteristics are shown to be in good agreement for SiGe HBTs, over a wide temperature range, thus proving the practical utility of the semi-numercial approach. The approach is next applied to submicron fully depleted (FD) and non-fully depleted (NFD) SOI MOSFETs via SOISPICE to analyze floating-body effects on circuit performance. New continuous charge models for scaled devices, having continuous derivatives in all regions of operation as well, are presented. The floating-body 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 floating-body 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 simulation-based design of a 10-bit current-steering digital-to-analog converter (DAC) is described to demonstrate a novel approach to ensure kink-free operation of floating-body NFD/SOI analog circuits.
Low-power and high-speed integrated-circuit (IC) technology is the key to future portable and wireless communication systems [Ish97]. Scaling of conventional complementary-metal-oxide-semiconductor (CMOS) and bipolar devices has been the common way to advance the technology. However, as the devices are scaled towards tenth-micron and smaller dimensions, a significant tradeoff between performance and reliability has to be made. For example, scaling of the MOS device to deep-submicron regimes (e.g., 0.1 gm) imposes tremendous challenges to device designers mainly due to hot-carrier effects (HCEs) and short-channel 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 threshold-voltage 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 hot-carrier-related reliability problems can be improved with reduced power supply voltage, but then SCEs become a dominant factor limiting device performance in low-voltage/low-power 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
performance due to decreased mobility and increased junction capacitance [Yan92]. Use of thin-film Silicon-On-Insulator (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 bulk-Si technology due to the reduced source/drain parasitic capacitance, simpler device isolation, absence of latch-up, increased radiation hardness, feasibility of resistors and capacitors free of parasitic junction effects, and less noise coupling afforded by high-resistivity substrates [Yeh95, Suh95]. In fact, various circuits utilizing deep-submicron 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, giga-bit DRAM and SRAM, and microwave RF circuits.
Also, the recent advances in bipolar technology, such as sustained device scaling, low-temperature epitaxial-base process [Gan90], and Si/SiGe heterojunction [Kon96, Jos96] have significantly extended the utility of the Si-based bipolar junction transistor (BJT). The graded SiGe-base 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 Si-based 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 technology-CAD (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
compact device models like BSIM3 (for MOSFETs) and Gummel-Poon (for BJTs). Numerically intensive device simulation, based on the process simulation, is done for device characterization, and then model-parameters are optimized for the circuit simulation through rigorous parameter-extraction 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 turn-around time and loss of correlation between predicted circuit performance and device structures and process flow. For example, the thin-film nature of the SOI MOSFET can underlie mechanisms that complicate circuit simulation and portend equivocation in design [Fos97]. Mainly because of floating-body effects, which in fact obtain even in body-tied structures because of unavoidable high resistance, empirical compact models like those used for designing bulk-Si CMOS circuits may not be useful for SOI. Not only are such models deficient, but they furthermore cannot be easily calibrated using common parameter-optimization techniques because of equivocal data acquisition implied by SOI device selfheating in DC measurements and the floating-body effects in pulse measurements designed to avoid the self-heating [Fos97].
The efficiency and reliability of the IC TCAD system could be greatly enhanced by a pragmatic mixed-mode device/circuit simulator based on physical, compact models which have process-based parameters that relate directly to the device structure. With this objective, MMSPICE (for bulk-Si and SiGebase bipolar and CMOS) [Jeo89, Jin92, Hon94] and SOISPICE (for SOI CMOS) [Vee89, Suh95, Yeh95] were created. Here, we further extend the process-based device models for general circuit applications to demonstrate the utility of this semi-numerical mixed-mode simulator. Further, the resulting enhanced models are used in various device/circuit simulations in actual technologies to
demonstrate viability and to reveal unique simulation and design issues for scaled SOI CMOS and SiGe-base HBT circuits.
In Chapter 2, a purely numerical way to improve the efficiency of the conventional TCAD system is examined. To pragmatically account for non-local carrier heating and hot-carrier effects such as velocity overshoot and impact ionization in multi-dimensional numerical device simulation, a new simplified energy-balance (SEB) model is developed and implemented in FLOODS [Lia94a] as a hierarchical option. In the SEB model, the energy-relaxation length is estimated from a pre-process drift-diffusion simulation using the carrier-velocity 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 two-dimensional SHD and full-HD DC simulations of submicron bulk-Si and SOI MOSFETs. The SHD simulations yield detailed distributions of carrier temperature, carrier velocity, and impact-ionization rate, which agree well with the full-HD 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 full-HD, 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 three-dimensional SHD device simulation as well.
In Chapter 3, a detailed analysis of temperature dependences in Si and SiGe-base HBT/BJTs is presented and implemented in MMSPICE to demonstrate the utility of semi-numerical coupled mixed-mode device/circuit simulation for bipolar ICs, and to gain physical insight regarding the performance of advanced bipolar transistors in temperature. Conventional Si bipolar
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 low-temperature operation, the device operation above room temperature is also examined for high-temperature applications. For verification, MMSPICE simulations in temperature are performed and compared with measurements based on recently developed HBT/BJT technologies.
In Chapter 4, the semi-numerical approach to mixed-mode simulation is applied to submicron fully depleted (FD) and non-fully depleted (NFD) SOI MOSFETs via SOISPICE to analyze the floating-body effects on circuit performance. SOISPICE was previously developed mainly for digital applications and verified as a reliable simulation tool for low-voltage and high-speed digital circuit design [Fos97]. Its physical, process-based models can be effectively used to characterize FD and NFD SOI MOSFETs. However, to apply the
models to analog and mixed-signal 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 small-signal parameters as well as large-signal 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 quasi-static 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 hybrid-FD/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 floating-body model is described based on the charge models, and is shown to predict well the frequency-dependent floating-body effect as well as DC and transient floatingbody effects.
In Chapter 5, the upgraded SOI MOSFET models are calibrated and used in circuit-simulation 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 next-generation 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 floating-body 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 floating-body effects. The dynamic floating-body effect associated with the sense amplifier can be solved using the FD device because it effectively decouples the body voltage from the
channel current. However, like the NFD device, transient leakage current is observed due to the capacitance coupling. Next a 10-bit current-steering DAC design is presented. FD or body-tied NFD SOI devices have been used almost exclusively for analog/mixed-signal 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 threshold-voltage 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 floating-body NFD/SOI MOSFET has not been seriously considered as an option for analog circuits since it suffers from the kink effect. Recently, NFD/SOI circuit-design approaches to avoid the kink effect, i.e., high-frequency 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 mixed-signal circuits. Here a 10-bit current-steering DAC based on floating-body NFD, bodytied NFD, and FD SOI MOSFET technologies is presented, and a novel approach to ensure general kink-free operation of floating-body NFD/SOI analog circuits is proposed. SOISPICE, with its device models upgraded and refined for high-frequency 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 one-dimensional hot-carrier transport model is developed and implemented in FLOODS to account for non-local carrier heating and velocity overshoot for device simulation. The model is verified through simulation of several advanced device structures.
SIMPLIFIED ENERGY-BALANCE MODEL FOR PRAGMATIC
MULTI-DIMENSIONAL DEVICE SIMULATION
In contemporary deep-submicron 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 drift-diffusion (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 non-local 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
band [Bor90], non-Maxwellian energy distribution [Che92], and higher-order moments of the BTE (e.g., heat flux) [Bor91]. Also, it suffers from relatively much longer simulation run-time 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 drift-diffusion (ADD) model (Che91], in which non-local effects are accounted for by introducing additional analytic field-dependent terms in drift-diffusion current equations, has been proposed as a way of efficiently extending the utility of drift-diffusion 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 higher-order derivatives of the electric field are significant, which is not uncommon in small devices. In the Appendix, we present an enhanced ADD model in one-dimensional 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 energy-balance equation has been recently proposed for post-processing use. The post-processing technique, in which electron energy is calculated from the electric field distribution obtained from a DD simulation, has been primarily used to predict the energy-dependent impact-ionization rate [Cra90, Slo91, Sou93, Ago93]. Slotboom et al. [Slo91] simplified the energy-balance 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.
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 energy-dependent two-dimensional impact-ionization, or substrate current model for the simulation of submicron MOSFET's, in which the simplified one-dimensional energy-balance equation proposed by Cook [Coo82] is applied to many current contours in order to generate a two-dimensional representation of carrier energy. However the average energy predicted by this post-processing technique is not correct in highly scaled devices since it decouples the electric field from the energy. For example, in deep-submicron MOSFETs, velocity overshoot influences the electric-field 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 impact-ionization current because of the lower electric field near the drain. The significant electric field-energy 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 energy-balance 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 base-emitter 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
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 energy-balance 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 non-local effects with reasonable accuracy. The new SEB model is formulated for practical multi-dimensional device simulation and is implemented in the two-dimensional tool FLOODS (FLorida Objected-Orient Device Simulator) [Lia94a] as a pragmatic option. The physical models underlying the SEB formalism are also discussed. Nonlocal effects in the submicron bulk-silicon 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 three-dimensional and transient device simulation, is addressed as well.
2.2 Hot-Carrier Transport Theory
2.2.1 Carrier Transport Equations
The energy distribution of hot carriers plays a crucial role in predicting the long-term reliability of advanced devices. For example, in a deep-submicron 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
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 semi-classical dynamics, under the assumptions that carriers behave as classical particles and collisions are instantaneous [Lun90]:
af +* (Vrf)+Fe (pf) = afc (2-1)
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 second-order ones that relate to the average carrier density, velocity, and energy, respectively, are usually sufficient to characterize semiconductor devices. The zeroth-order moment equation is the familiar carrier continuity equation. In the steady-state condition, this equation for electrons is given by
Ve(n ,) + (Rn Gn) = 0 (2-2)
where n is the electron concentration, and (Rn-Gn) is the recombination/ generation rate; un now represents the average electron velocity. Note that (22) physically describes the conservation of the electron density.
The first-order moment equation is usually referred to as the momentumbalance equation. In the steady state, this equation for electrons becomes
kB kTnVn (kBn) (2-3
4n -- n E n + V (2-3)
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 (2-3) and assume the fielddependent mobility model, the momentum-balance equation reduces to the well known drift-diffusion 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 drift-diffusion model is not adequate for investigating modern scaled devices.
A better approach is to include the temperature gradient term in the momentum-balance equation and to use the energy-dependent mobility model, as is done in (2-3). 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 steady-state condition, is given by
S* Vwn = (n E) V(nkgTnn + Q) n (2-4) n 'r."
where wn is the average electron energy, Q is the heat flux, and Tw is the energy relaxation time. In (2-4), wn and wo are defined as
3 1 2
Wn = kBTn + -m n (2-5)
wo = kBTL (2-6)
where TL is the lattice temperature and m* is the effective electron mass. The left-hand side of (2-4) 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 right-hand 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 higher-order moments. The zeroth-order moment equation includes the first-order moment, the average velocity. The first-order moment equation includes the second-order moment, the carrier average energy, or temperature. And the second-order moment equation includes the heat flux, the third-order 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 heat-flux term, which represents the contribution of the third moment to the energybalance equation, is commonly assumed as the phenomenological Fourier law,
Q = -KVTn (2-7)
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 moment-equation approach is far simpler for characterizing the devices than the direct-solution 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 higher-order 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 Non-Local 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 hot-carrier 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 non-local carrier heating can be easily explained quantitatively by using a simplified one-dimensional energy-balance equation proposed by Slotboom [Slo91]. The equation can be derived from (2-4) by assuming that the effect of the heat flux and diffusion are negligible and v, = U1sat:
dwn W O -( 3
+ = E (2-8)
dx kw 5
where the energy relaxation length X which is defined as 5t, sat/3, is set to a constant. If we iteratively solve (2-8) for the carrier energy, we obtain
3 d (3 d d
W -(WwEE- XW-x( E- wwE+(...
3 3 2 dE 3 3 d2E (3)
= E- 5 wx -+ 5Xw dx2 +0 (2-9)
As can be seen in (2-9), 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 right-hand-side term in (2-9)) 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 carrier-transport models involves a trade-off 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, full-energy 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 high-field carrier transport in devices. It is also difficult to
Monte Carlo Full Energy Band Structure] Scattering
Boltzmann Transport Eq.
Non-Parabolic Band Hydrodynamic
Carrier Continuity Eq. Degeneracy
Momentum Balance Eq.
Energy Balance Eq. Energy Relaxation Time
Computational Simplified Energy Balance Efficiency
Carrier Continuity Eq.
Fig. 2.1 The hierarchy of numerical simulation.
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 non-local transport and hot-carrier 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 high-energy electrons, the energy-band structure becomes non-parabolic and the energy distribution is not Maxwellian. These modifications introduce additional terms into the hydrodynamic models since they basically depend on higher-order 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 energy-balance 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 one-dimensional 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 energy-dependent impactionization rate. Therefore, in the next sections, we develop a new coupled twodimensional simplified energy-balance model and demonstrate its utility.
2.3 Simplified Energy-Balance Model
2.3.1 Model Development and Its Implementation in FLOODS
We here propose a new simplified energy-balance (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.2-gm device, the drift energy can be nearly 15% of the total kinetic energy [Ste93]. We next assume that the heat-flux term Q is negligible in determining the electron temperature. This term in (2-4) 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 non-physical 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 heat-flux term is frequently and justifiably neglected in the energy-balance model [Go188, Slo91].
Now focusing on (2-4), we assume that V*(nn) 0, which, as (2-2) shows, is equivalent to neglecting carrier recombination/generation in regions where carriers are heated. Then (2-4) simplifies to
(VTn TL) 2qn E. (2-10)
n n 5 rw 5k n
Expressing the electric field on right-hand side by the gradient of the conduction band edge EC, which renders (2-10) applicable for homostructures as well as heterostructures, we rewrite the steady-state SEB equation as
V Tn + + "-E~ ] (TnTL) = 0. (2-11)
In (2-11), we have implicitly assumed that the effect of diffusion is negligible since the electrons are heated mainly by drifting.
For numerical device simulation, (2-11) can be spatially discretized using the box method. For the two-dimensional case on a triangular mesh structure as illustrated in Fig. 2.2, the discretization gives
lij( E + Tnwn, ij+A(T,- To)i = 0 (2-12)
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 energy-relaxation length defined on triangle edge (i-j) of the grid, which is perpendicular to lij. The energy-relaxation length, first defined by Slotboom [Slo91] with constant uD = usat, is defined here as
1wn = 3Twn7n (2-13)
where un is spatially varying. In the first term of (2-12), 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 energy-balance equation greatly increases the nonlinearity of the entire system, and thus the computational time during
Fig. 2.2 A two-dimensional structure with triangular
simulation. To reduce this computational burden, Slotboom assumed a constant velocity throughout the device and applied the model to predict the energy-dependent impact-ionization rate along a single current path in such devices as BJT and MOSFET. However the model does not provide detailed information about the two-dimensional distribution of the electron energy within the devices. Recently, Agostinelli and Tasch [Ago94] demonstrated an energy-dependent two-dimensional substrate current model for the simulation of submicron MOSFETs in which the simplified one-dimensional energy balance equation proposed by Cook [Coo82] was applied to many current contours in order to generate a two-dimensional representation of electron energy, and thus impact-ionization rate. However, as we noted in the introduction, these models are only available as post-processing in device simulation, and thus exaggerate the electron energy, and hence overly predict the impactionization rates.
In this section, we present an energy-balance 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 high-level 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 energy-balance equation. In MEDICI [Tma92], a decoupled approach for the
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 energy-balance equations for the carrier energies. The carrier energies are used to re-evaluate I, n, and p. This procedure is repeated until the system converges. However due to the decoupled nature between the energy-balance equation and DD equation, convergence of the HD simulation in MEDICI can be slow for highly scaled devices. Alternatively, FLOODS has chosen the full Newton-like approach with the DD simulation as a pre-processor. 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 post-processing analytical methods for estimating carrier temperatures, e.g., the one proposed by Slotboom [Slo91]. According to Slotboom's model, the electron energy-relaxation 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 one-dimensional devices, it is not physical and it is not sufficiently accurate for use in general multi-dimensional device simulation.
DD Simulation (W, n,p)
which define on each edge
SHD Simulation; Poisson's Eq. (W) Carrier Continuity Eqs. (n, p)
SEB Eqs. (Tn Tp)
Set New Bias
Fig. 2.3 The flow chart of the SEB/SHD model for a specific
Alternatively in our SEB model, the energy-relaxation length in (2-13) is estimated from a pre-process drift-diffusion 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 (2-11) for electrons, are solved simultaneously. The carrier energies are thus directly coupled to and defined by the electric field distribution; the pre-defined energy-relaxation 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 object-oriented 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 energy-balance equation is replaced by the SEB equation, i.e., (2-12) with Xwn in (2-13) 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, impact-ionization rates, and energy-relaxation times have been used without change. The energy dependence of the carrier mobility is incorporated through Klaassen's model [Kla92].
2.3.2 Hot-Carrier Effects
In deep-submicron 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 hot-carrier effects since the effects are fundamentally described by the carrier energy, not by the electric field. In this section prominent non-local effects caused by hot-carriers such as velocity overshoot and impact-ionization 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 non-local 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 impact-ionization rate exist in the literature [Oga65, Bar91, Sou93]. In FLOODS, the Scholl-Quade 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
Gn = -jexp -erfc (2-14)
0o 7 w ;
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 (2-14), 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 (2-2). 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 impact-ionization as well as other hot-carrier effects directly in the coupled numerical iteration process. The SEB/SHD model is not merely decoupled post-processing.
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 n-channel 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 cm-3. An LDD with 0.1 gim depth and 1.0 x 1018 cm-3 doping density is part of the device structure.
Fig. 2.5 illustrates the electric field distributions along the channel simulated with the DD and full-HD models of FLOODS at VDS=5.0V and VGS=5.0V. The electric field and field gradient near the drain predicted by the
S D n+ poly n+ n- 0.5m n+
1.4 x 1017 cm-3
Fig. 2.4 The simulated MOSFET structure. The gate
oxide thickness is 15 nm, the S/D junction depth
is 0.3 gm.
350 -O DD
6 250 S 150
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.
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 energy-relaxation time (tw = 0.35ps). Note that the strong non-equilibrium behaviors are predicted near the drain end of the channel by both models. The hot-carrier 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.5-gm 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 source-drain transition time, but perhaps more significantly tends to reduce T, for a given bias condition.
Fig. 2.8 shows the two-dimensional impact-ionization generation rate simulated with the SEB/SHD and full HD models at VGS = VDS = 5.0 V. The
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.
SHD u 1.5
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.
< 20e 15
) 0.4 0 0.8 ".2 0.4 (6V. )0 0 0.2
Fig. 2.8 (a) Generation rate by the impact-ionization process
simulated with the HD model of FLOODS.
C) 25 20.
00..0.2 0.4 0 0 0 0.2 riss
Fig. 2.8 (b) Generation rate by the impact-ionization process simulated with the SEB/SHD model of FLOODS.
values used for parameters Wth and ~oI in (2-14) are 4.0 eV and 7 x 1015 s-1, respectively. As seen in (2-14), the impact-ionization 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 two-dimensional generation rates were calculated via the post-processing option, which is adequate for the common case of weak impact ionization. As we mentioned in the previous section however, (2-14), inserted in (2-2), 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.4x10-8 A/gm and 3.0x10-8 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 carrier-velocity expressions), and the energy-balance 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 energy-relaxation 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
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
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 hot-carrier effects are prominent.
2.4.2 Scaled SOI MOSFET
Unlike a bulk MOSFET, there is no substrate contact in the thin-film (commonly floating-body) SOI MOSFET, and thus majority carriers generated near the drain by impact-ionization are injected into the film body where they accumulate. This charge modulation forward-biases the source junction and thus activates the parasitic bipolar transistor [Fos90]. Therefore it is very difficult to characterize the impact-ionization 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 hot-carrier effect in a deep-submicron fully-depleted (FD) SOI MOSFET.
Fig. 2.10 shows the simulated two-dimensional n-channel 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
n+ 1.1x1017 cm-3 n+ Back Oxide Si
tox,=10nm, tsi = 62 nm, tbox=0.38 pm
Fig. 2.10 The simulated n-channel SOI MOSFET structure.The S/D junction depth is 0.3 gm, and the channel length is 0.3 uim.
0.3 pm. The doping concentration in the body is 1.1 x 1017 cm-3 and the doping in the source and drain regions is 5 x 1019 cm-3. To focus on the hot-carrier effects, the lattice temperature is set to 300 K, and self-heating 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 non-local 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 full-HD simulation is eliminated in the SHD option. Finally, Fig. 2.12 shows the two-dimensional carrier generation rate due to the impact-ionization 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 three-dimensional 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 three-dimensional simulator, the box method, analogous to the two-dimensional one, is commonly used to discretize the carrier transport equations, which are integrated over a small volume enclosing each node [Tma92]. The discretized three-dimensional SEB equation has the same form as (2-12) with minor changes; length and area in (2-12) 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 energy-relaxation length distributions can be evaluated accordingly from the preceding DD analysis. This capability of extending to the three-dimensional case, potentially with even more
(-0 SHD 0 DD
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
0.5 S 002 0.2 0.3 0
0.0 0 0 0. ace
Fig. 2.12 Generation rate by the impact-ionization process simulated with the SEB/SHD model of FLOODS.
enhancement in relative computational efficiency, obviously increases the utility of the SEB/SHD model.
So far attention was given to DC, or steady-state 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 time-dependent energy-balance 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 (2-15)
5 Ir dt 5k n
which can be expressed as
dt= Fn(E, T,). (2-16)
Now the time-derivative in (2-16) must be discretized, considering stability and local truncation error in the selection of the numerical integration scheme. First-order backward difference formula (BDF1), second-order backward difference formula (BDF2), trapezoidal rule (TR), and BDF2-TR 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, (2-16) can be discretized as follows:
Atk = Fn(Ek, T) (2-17)
where At = tk tk -. The use of (2-17) at each mesh point follows the spatial discretization of (2-16) as indicated by (2-12), 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 order-of-magnitude improvement in run time for general two-dimensional transient MOSFET simulation afforded by the SEB/SHD option; and even more for three-dimensional numerical simulation because of the inherent higher degree of nonlinearity in the HD system.
An efficient multi-dimensional SEB/SHD model was developed based on assumed fixed energy-relaxation 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 post-processing analytical models and de-coupled 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 two-dimensional SEB/SHD DC model was implemented in FLOODS as a pragmatic option without substantial source-code modification. All physical models such as impact-ionization and mobility are used without change. Model verification was accomplished through submicron bulk-Si and SOI MOSFETs simulations. The SHD simulation provides reliable detailed
two-dimensional representations for carrier temperature, carrier velocity, and impact-ionization 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 three-dimensional 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.
IMPLEMENTATION OF TEMPERATURE DEPENDENCE IN MMSPICE HBT/BJT MODEL
MMSPICE was first developed to incorporate the novel concept of seminumerical mixed-mode device/circuit simulation [Jeo89], and was verified, through subsequent enhancement [Jin92, Hon94], to be a predictive and reliable IC TCAD tool. Since the charge-based 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 MMSPICE-3.0 [Hon94] is a physical but compact model for the SiGe-base heterojunction bipolar transistor (HBT), which has been widely investigated in the past few years [Cra93, Gha91, Gib88] as a viable technology due to excellent current-driving capability and high-speed performance. Although this new version of MMSPICE simulated HBT/BJT devices and circuits successfully at room temperature, it was not applicable to temperature-dependent simulation. Recent interest in low-temperature 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
current gain reduction is larger bandgap narrowing in the heavily doped emitter than in the base. In thin-base 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 freeze-out which degrades switching performance.
The recent advances of bipolar technology such as sustained device scaling [Com91], low-temperature epitaxial-base 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 freeze-out in the base is diminished by using a thin base heavily doped above the Mott transition. This is more prominent in epitaxial-base technology which is capable of achieving significantly thinner base width than an ion-implanted base process.
This chapter describes a detailed investigation of the temperature dependence for these advanced Si- and SiGe-base HBT/BJTs, modeling of the
temperature dependence, and implementation of the modeling into MMSPICE3.0. Although most concerns are given to the low-temperature operation, the device operation above room temperature is also examined for high-temperature 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 SiGe-base 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 quasi-neutral 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 Moll-Ross current density expression can be generalized in the forward-active mode as follows [Kro85]:
Jco = e (3-1)
where nib, the effective intrinsic carrier concentration in the base, is generally position-dependent and given by
2 2 AEgb, HD + AEgb, Ge
nib = nioexp kT (3-2)
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 (3-2) into (3-1), 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 = (3-3)
WBNAO expil'- 1
= Jco'(T) f (', T) (3-4)
where Dn is the average electron diffusion coefficient in the base, WB is the quasi-neutral 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 one-dimensional SiGe-base HBT structure is assumed in the derivation of (3-3), 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 (3-3), Jco is dependent on temperature both explicitly and implicitly through the several physical parameters. The most dominant term determining the temperature dependence of
the collector current is nio. As will be seen in (3-20), it is an exponential
function of temperature, and is divided by four for every 80C decrease near room temperature. This Jco degrades significantly with cooling.
Also, the freeze-out 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 (3-1), 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 (3-3), 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 (3-3), which represents the multiplication factor of the collector current corresponding to both the bandgap narrowing and built-in field: For the advanced Si BJT with an ion-implanted 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). (3-5)
Although (3-5) 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 (3-6)
Since AEgb(grade) has a positive value for a steep Ge gradient across the base, AEgb(O) is reduced as the emitter-base 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].
3.2.2 Base Current
The collector current given in (3-3) to first order depends only on the base profile. Unfortunately, the same cannot be said of base current, which depends on emitter contact and emitter-base 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 emitter-base 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]
qno 1 WH u t (Age qVBE/kT
J = + 1 + exp (3-7)
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
is dependent on temperature mainly through both nio and AEge at all temperatures.
The non-ideal base current at room temperature under low-forward bias is in general caused by the Shockley-Read-Hall (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 field-assisted Poole-Frenkel SRH recombination, band-to-band tunneling, and trap-assisted tunneling are responsible for this large non-ideal factor. The Poole-Frenkel effect, which is described by the enhancement of the
carrier emission rate due to the large junction electric field, can be written as [Hac88]
IP-F = ISEO exp/ 2E1/2 exp (qVBE 1 (3-8)
where E is the base-emitter junction electric field and 3 is q/ni. Since E is a positive function of Obi VBE, IP-F has an exponential dependence on VBE weaker than exp(qVBE/2kT), which implies a nonideality larger than two.
Emitter-base 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 temperature-dependent and thus does not decrease as rapidly with cooling as diffusion or recombination currents. Since silicon is an indirect band-gap material, any observed tunneling current under forward bias is either phonon-assisted or trap-assisted (excess current) [Hac86]. Recent studies on heavily doped junctions with doping levels between 1018 and 1019 cm-3 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 (3-9)
The current gain degrades exponentially with decreasing temperature for conventional homojunction BJTs since the emitter is far more heavily doped
than the base, and thus the exponent in (3-9) 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 high-performance device for analog circuit applications, high Early voltage is most desirable since it allows large small-signal output impedance. The Early voltage VA is formulated directly from the variation of collector current with collector-base bias, based on the generalized Moll-Ross relation [Kro851:
DIc p(WB) w B p dx]- dWB IC VcB D(WB)nb(WB)o Dnnib dVc VA (3-10) which can be simplified to
VA = qDn(WB)n(WB)fW dx. (3-11) CCB 0 Dnnib
Substituting (3-2) into (3-11), a closed-form expression for the Early voltage can be obtained:
qWBNAO 1 exp(-q') (AEgb(grade)
VA CCB r' K T (3-12)
In the deviation of (3-12), 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 collector-base 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 Moll-Ross extension, as [Kro85]
S= dy dz (3-13)
=o fP [z Dnnib J
Substituting (3-2) into (3-13), we can obtain a closed-form expression for a device with a triangular Ge profile:
TB W-- 1 (1-e-)). (3-14)
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 ion-implanted base doping profile. If AEgb(0) <
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
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 Ge-based bipolar transistor, 'r' is much larger than one, and thus TB can be simplified to
TB = 1 (3-16)
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 (3-17)
1 \b TBDn 1.5 [cm-2
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.
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 (3-17) increases drastically with cooling due to the bandgap narrowing in the heavily doped emitter. In a highly graded SiGe-base 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 ultra-thin base bipolar transistors the emitter transit time is comparable to the base transit time.
3.2.7 Cutoff Frequency
The cut-off frequency fT is defined as the frequency at which the common-emitter short-circuit 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 (3-18)
where TEB and 'BC,scl are the emitter-base and base-collector depletion layer charging time, respectively, and TC is the quasi-neutral 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
TEB = C (CEB + CBC) (3-19)
where CEB and CBC are the emitter-base and base-collector depletion-layer capacitances, respectively. These capacitances decrease monotonically from
300K to 77K by roughly 20% due to the increases in the built-in voltage for both Si and SiGe devices. Also, as can be seen in (3-19), TEB is proportional to the temperature and thus favorably influenced by cooling. The last delay term in (3-19) is the base-collector depletion-layer transit time, tB ,scl = WBC, scl/2vd where WBC,scl is the base-collector 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 (3-19), 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 MMSPICE-3.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 physics-based 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 temperature-dependent constant equal to the square of the intrinsic carrier concentration:
Table 3.1 MMSPICE-3.0 HBT/BJT Model Parameters
Name Description Units Typical
UNEPI Electron mobility in epi collector cm2*V-1*s-1 1.0e3 NEPI Epi doping density cm-3 1.0e16 WEPI Epi-collector width m 5.0e-7 WBM Metallurgical base width m 2.0e-7 NAO Extrapolated peak base doping cm-3 1.0e18 ETA Base doping gradient factor 3.0 AI Pre-exponential coeff. for impact ionization rate cm-1 7.03e5 BI Exponential coeff. for impact ionization rate VOcm-1 1.23e6 VS Electron saturated drift velocity cm*s-1 1.0e7 JEO Emitter saturation current density Aom-2 1.0e-8 JSEO B-E SCR saturation current density A*m2 1.0e-4 NEB B-E SCR emission coefficient 2.0 WSEO Zero bias B-E SCR width m 5.0e-8 PE B-E junction potential barrier V 1.0 ME B-E junction exponential factor 0.4 PC B-C junction potential barrier V 0.8 DNB Electron diffusivity in base for low injection cm2es-1 10.0 PS C-S junction potential barrier V 0.6 MS C-S junction exponential factor 0.4 CJS Zero-bias C-S junction capacitance Fom-2 1.0e-4 CRBI Intrinsic base resistance coefficient V*s*cm-2 2.0e-4 TB Carrier lifetime in base s 1.0e-7 TC Carrier lifetime in collector s 1.0e-6 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 Base-charge partition factor 0.5 FC Collector-charge partition factor 0.5 TE Emitter hole transit time s 1.0e-9 UPBASE Hole mobility in intrinsic base cm2*V1*s"1 230. NEBP Peripheral B-E SCR emission coefficient 2.0 JEOP Peripheral B-E SCR saturation current density A/m 2.0e-11 LBE Width of extrinsic base region m 1.5e-6 JCOP Substrate saturation current density A/m 1.5e-15 DNBH Electron diffusivity in base for high injection cm29s-1 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
2 = (27kT 3 3/2 3/2 -Eg/kT
noPo = nio =4 h2 -. ) mnmPe (3-20)
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
nio(T2) (T2 2 Eg(T2) Eg(T,) (3-21)
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]
Eg(T) = E(O) T+ (3-22)
with a=7.02 x 10-4, P=1108 and Eg(0)=1.16eV.
As the doping concentration increases above 1017 cm-3 the silicon band structure changes due to the appearance of band-tails 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
an increase in the equilibrium carrier concentration product n oP = nie over
nio which can be formulated by introducing a phenomenological model parameter AEg, HD as can be seen in (3-2). This apparent bandgap shrinkage parameter is a strong function of doping concentration, but its temperature dependence is negligible.
(b) Base Current Components
The back-injected hole current is the dominant base current component under medium emitter-base bias condition, and was given by equation (3-7). 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 (3-23)
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 MMSPICE-3.0, we introduce a parameter, DEGA, to represent AEga, which is the only additional parameter had to be added for the temperature-dependent simulation.
As mentioned in section 3.2.3, since the B-E 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
In advanced scaled-down 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 (3-25)
where PE is the length of the base-emitter 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(-) (3-26)
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 generation-recombination current along the base-emitter 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 field-assisted Poole-Frenkel SRH recombination, band-toband tunneling, and trap-assisted tunneling, which can be dominant at low temperature, are not modeled in current version of MMSPICE.
(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 majority-carrier mobility varies more rapidly with temperature for lower doping concentrations, where the strongly temperature-dependent lattice scattering mechanism is dominant. For very high doping concentration where ionized impurity scattering dominates, the mobility is almost temperature-independent. 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 temperature-dependence modeling. From the model by Arora, the majority-electron mobility is expressed as
(T) 88 57 7.4 x 108 T-2.33
n, maj(T) = 88 + (3-27) 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 epi-collector concentration. Recent experimental results have shown that, starting at a doping concentration of 1017 cm-3, the minority-carrier mobility in silicon exceeds the majority-carrier mobility, up to a factor of three at a concentration of 1020 cm-3 [Swi86]. So different models are needed for minority- and majority-carrier mobilities. Empirical models for the doping dependence of minority-electron mobility have been defined [Swi861
and extended to incorporate the temperature dependence. From this model, the temperature dependence of the minority-electron mobility is written as
,mi,(T)= 232T 0.57 6.97 x 108 T-2.33
(T) = 232 + (3-28) 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 majority-electron 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 minority-electron mobilities which have trends similar to with the majority-electron mobility, but have almost a factor of two higher values for heavily doped silicon. Similarly, the majority-hole mobility in the intrinsic base region can be expressed as [Aro82]
p, mj(T) = 54 T -0.57 1.36 x 108. T-2.333-29)
1 + 17 2 .88(2.35 x 107 (T/300)2.4 300
Note that the mobility expressions in (3-27), (3-28), and (3-29) 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 (3-27), the temperature dependence of the majority-electron mobility in the epi-collector 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(3-30300
N = 1016 cm3
N = 10 cm-3
Eln 10 cm2V-s ]
N =1018 cm-3
50 100 150 200 250 300 350 400
Fig. 3.2 Majority (solid) and minority (dashed) electron mobilities
as a function of temperature for various impurity
I+ NEPI 1-0.146
1 + .0.88/ --1
1.2 x 1017 (T,/300)24 300
F(T1, T2) = x --I ) (3-31) NEPI .088T2 0.146
1.2 x 1017. (T2/300)24 0
Similar expressions for the temperature dependences of the minority-electron mobility and the majority-hole 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( (3-32)
From above equation, the temperature dependence of the saturation velocity in MMSPICE is determined by
VS(T2) = Vy(TI). (3-33) Vtanh(175/T 1)
As can be seen in equation (3-2) and (3-10), 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
Dn = 2kTF(l/2 EkT Fn -1/2(EcEFn (3-34)
_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 built-in potential across a emitter-base junction under zero bias can be expressed by
OE = kTIn (3-35)
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
where AEg is the spatial average of the effective bandgap narrowing in the space charge region of the base-emitter junction. Generally, OE increases with decreasing temperature due to significant reduction of the intrinsic carrier concentration. Similar expressions can be used for both base-collector and collector-substrate built-in potentials.
(e) Recombination Lifetimes
Recombination lifetimes in the base and in the collector can be largely determined by the Shockley-Read-Hall (SRH) recombination process since the
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 (3-37)
where A is a constant determined by the doping concentration and AE is the activation energy for the minority-carrier 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( (3-38) (T 0( I (3-38)
3.3.2 Model Implementation in MMSPICE-3.0
Unlike the Gummel-Poon 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 temperature-dependent 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.
The default simulation temperature in MMSPICE is 270C, for which user-specified 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), B-E SCR saturation current density (JSEO), B-E equilibrium SCR width (WSEO), B-E junction potential (PE), B-C junction potential (PC), electron diffusivity in base (DNB), substrate junction potential
(PS), zero-bias collector-substrate 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 MMSPICE-3.0 has been verified based on the most advanced devices including two SiGe-base HBTs from IBM and one lesser scaled Si-base BJT from Motorola. The first SiGe device, SiGe1 [Cre93], was fabricated using the low-temperature epitaxy process which have great potential for achieving the very-thin, 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
3.7x1018 cm-3 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 self-aligned structure with a reduced thermal-cycle 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, MOSAIC-3 [Zde87], is a self-aligned BJT with polysilicon emitter and oxideisolation. The device has a peak base doping concentration of 1x1018 cm-3, a base width of 150 nm, and a collector doping concentration of 2x1016 cm-3. 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 temperature-dependent behaviors of the collector current over a wide temperature range quite well, except for extremely low temperature at which the carrier freeze-out effect is prominent. Since the current version of MMSPICE does not include the proper model for the carrier freeze-out effect, a discrepancy is found between the simulation and measure-
VCB = 0 V 101 WB=56nm
10-3/ / /
Ic 10-5 [A]
10 182K 133
10 84 '
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 emitter-base voltage for SiGel. This figure
shows an improved transconductance with decreasing
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 collector-current 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 emitter-base voltage curve, increases with cooling, as discussed in section 3.2.1.
Fig. 3.4 illustrates the base current as a function of emitter-base voltage for various temperatures. Under moderate injection levels, the back-injected 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 hard-wired 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 temperature-dependent simulation of MMSPICE. Under lower- injection levels, SRH recombination in the space-charge region of the emitter-base junction is a dominant mechanism
10-2 WB=56nm VCB = 0 V
107 250 / 182K
10 13 /
19 ,', 84K
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 emitter-base voltage for SiGe at a variety of
responsible for the base current near room temperature. However, other parasitic base leakage mechanisms such as tunneling and P-F 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 current-voltage characteristics for SiGe1 simulated with MMSPICE-3.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 common-emitter 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 (M-1) 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 M-1. This increase in the multiplication factor with cooling offsets the reduction of P and yields temperature-insensitive BVCEO. Therefore, to predict the variation of the common-emitter breakdown voltage with temperature, a accurate temperature model for M-1 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 low-temperature analog applications.
0.4 182K 0.3 85K
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.
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 emitter-base 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 low-temperature operation for SiGebase HBTs, partly because of the growing interest in the low-temperature operation of SiGe-base 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 MOSAIC-3 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
WB=56nm VCB = 1 V
40 309K [GHz ] 30
100 10 10
Fig. 3.6 Measured (dashed) and simulated (solid) transistor cutoff frequency as a function of collector current at 309K and 85K
NB=4.0xl018cm-3 VCB = 0 V
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.
AT = 25 K WB=150nm
VCB = 1 V NB=1.0x1018cm-3
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
emitter-base voltage versus temperature curves for various collector currents illustrated in Fig. 3.9. We believe the high-VBE (high-current) discrepancies are related to the temperature dependences of both internal and external series resistances as well as diffusivity under high-injection conditions.
A detailed analysis of temperature dependences of Si- and SiGe-base HBT/BJTs was presented and implemented into the existing bipolar transistor model in MMSPICE-3.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 MMSPICE-simulated 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.
IE= 6 mA
VBE 0.7 1 mA [V]
0.6- 0.1 mA lO pA
300 320 340 360 380 400 420
Fig. 3.9 Measured (dashed) and simulated (solid) emitter-base voltage versus temperature with various emitter currents for
PHYSICAL SOI MOSFET CHARGE MODEL FOR LOW-VOLTAGE
ANALOG AND MIXED-SIGNAL CIRCUIT APPLICATIONS
MOSFETs built on SOI have yielded superior performance in high frequency analog and mixed-signal applications mainly because of the inherent device isolation afforded by the underlying oxide, which reduces parasitic source/drain capacitance and allows use of high-resistivity substrates for lowloss transmission lines [Co196] [Hur96]. The thin-film nature of the SOI MOSFET however can underlie mechanisms that complicate circuit simulation and portend equivocation in design [Fos97]. Mainly because of floating-body effects, empirical compact models like those used for designing bulk-Si 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 low-voltage and high-speed digital circuit design [Fos97]. Its physical, process-based models can be effectively used to characterize the fully depleted (FD) and non-fully 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 large-signal characteristics of the SOI MOSFET. Therefore the SOI analog model should not only meet common requirements
for digital design such as reasonable I-V characteristics, delay prediction, and charge conservation, but also should give accurate values for all small-signal quantities such as transconductance and transcapacitances; and all of these parameters should be continuous with respect to any terminal voltage. Furthermore, floating-body 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 run-time relative to that of previous versions of the SOISPICE models. A unified floating-body model is described based on the charge models, which predicts well the frequency-dependent floating-body effect as well as DC and transient floating-body effects.
4.2 SOI MOSFET Charge Models
As illustrated by the network representation in Fig. 4.1, SOISPICE FD and NFD SOI models are charge-based with five terminals and they allow the floating-body 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
SdQGf ICH dt
RS+RLDS VLDD RD+RLDD S-R S ds D D
IRGt(VBS) B' -IRGt(VBD)
Fig. 4.1 Network representation of the charge-based SOI (nchannel) MOSFET models in SOISPICE.
and their first derivatives are continuous through all regions of operations. The current model is also modified accordingly.
The intrinsic quasi-static front-gate charge can be modeled using Gauss's law as
QGf = WCofL (VGfS- Ofs- Vsf (Y))dy (4-1)
where W is the channel width, Cof (= ,ox/toxf) is the front-gate capacitance per unit area, <, is the front-gate work-function difference, and Wf(Y) is the front-surface potential. Fig. 4.2 shows MEDICI-simulated 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 drain-induced barrier lowering is hardly noticeable. In such case, (4-1) can be simplified as
QGfW = CofWLe(VGfs ms sf) (4-2)
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 (4-3) oFB Cof(1 + cx)L
In (4-3), a = Cb/Cof where Cb = -s/tb, and Qb is the body charge density. The effective channel length, Le=L-yD-YS, is defined in weak inversion as the portion of the channel where carrier diffusion is predominant and where the
VGfS=O.OV (Weak inversion) 2.5 - - -- VGfs=.0OV (Strong inversion)
0.0 0.1 0.2 0.3 0.4 y (im)
Fig. 4.2 MEDICI-simulated potential distributions along the
surface at different drain and gate biases for a 0.2 gm
longitudinal electric field is small relative to the fields closer to the source and drain junctions [Suh95]; YD and ys account for a weak-inversion channellength modulation in the drain and source sides, respectively.
To characterize YD, the two-dimensional Gauss's law was applied near the drain, yielding a VDS-independent [Gre93, Suh95]
yD =icln 2 (4-4)
where l = 2Co + ) This VDs-independence 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 (4-4) 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 (4-5)
Note this simple fix is justified by MEDICI simulations, which suggest yD depends only weakly on VDS.
However, the expression for Le in (4-5) becomes invalid when VGfS goes to the accumulation region since (4-4) was derived assuming only depletion charges under the gate, which is reasonable in the weak-inversion 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 front-gate charge is defined in accumulation as
QGfA = CofWL(VGfs ms VBS) (4-6)
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 (4-7)
where A = A'/(CofWL) and the coefficient A' is a hard-wired 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 body-charge modulation by yD and ys, respectively, and defined as
QD = qtbNBODYWyD (4-8)
Qs = qtbNBODYWYS (4-9)
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 (4-6), but since the inaccuracies are small, these modifications were not done.
Next consider the strong-inversion case. Much of the recent analog modeling effort for bulk-Si CMOS has been aimed at removing the discontinuity in small-signal 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 (4-10) In(1 + eB)
where VDSAT is the drain saturation voltage and B is a constant which controls the stiffness of the smoothing function. Then an effective drain-source voltage VDSX can be defined as
VDSX = F VDSAT (4-11)
which smoothly varies between VDS and VDSAT.
The same smoothing function can be applied to the SOI modeling. In such a case, the front-gate charge for strong inversion is expressed as
lin sat (4-12)
QGf = GfS + GfS(4-12)
where QGfS is the gate charge component between y=O to y=L-AL and QGfS is the gate charge component in the saturation region from y=L-AL to y=L; AL is the modulated channel length in strong inversion [Suh95]. These charges, with (4-11), 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() (4-13) of
Gf WLC S W ms-sfS VDSX 2eff oI c- 1 (4-14)
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 strong-inversion region. The AL can be obtained by solving a differential equation based on the two-dimensional Gauss's law near the drain. The solution [Suh95] of the differential equation, with (4-11) is
VDS VDSx- 2 sinh + 2 .[cosh-- 1 (4-15)
4 eff .lc 7,
From (4-15), 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. 4-3 shows VDsx and (L-AL)/L as a function of VDS for a 0.2 lm n-channel 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 SOISPICE-4.3), also shown in the figure.
The channel-charge component of the strong-inversion drain and source charges are computed from the original formalism [Suh95], but with smooth VDSx and AL. First total channel charge is calculated by
QCH = CH+ QCH (4-16)
lin (2z (z 1) (4-17)
QlCH = -W(L AL)Co(1 + o)VDSX- 2z -1 +(u -z) (4-17)
QsaH = WALQc(L -AL) (4-18)
where QC(L-AL) is the channel charge at y=L-AL. The drain charge
QD( = QD + QD ) is then computed from conventional charge partitioning method [Vee88]
1.2 I I I
0.2 gm NFD SOI MOSFET (VGfs=1.5V)
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 SOISPICE-4.3) of the model.
lin 42( 1) 4 (z- 1)5 (u z)
QD =-WLCf(+ )VDSX[ 2z -1 15 (2z 1)2 (4-19)
QD = W AL c(Le). (4-20)
Then the source charge is given by
Qs = QCH- QD. (4-21)
So far we described charge modeling in the weak- and strong-inversion 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 low-voltage 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 SOISPICE-4.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 (4-22)
(VTw Vp) + J(VTW VP)2 (VTW 2Vp + VTS)(VTW VGfS) (4-23) VTW 2Vp + VTS
The intersection points, Vp and Qp in (4-22) and (4-23) and illustrated in Fig.
4.4, are calculated as follows:
VP = (VTsQTs' VTWQTw' QTS + QTW) (4-24) QTS'- QTW'
QP = QTW + (VGfs VTW)QTW'. (4-25) Henceforth, continuity of the charges and transcapacitances at the moderateinversion boundaries is ensured.
Gauss's law gives the back-gate charge under the intrinsic body region as
QGb = WLCob(VGbS ms VBS) (4-26) where
QB = -(Qs + QD + QGf + QGb) (4-27) The terminal charging/discharging currents indicated in Fig. 4.1 are then characterized using the quasi-static approximation:
dt = aV3 dt (4-28)