Интерны | Megalo Box | فيلم Bioscopewala 2018 مترجم مشاهدة وتحميل فيلم الدراما الهندي Bioscopewala 2018 مترجم بجودة 1080P WEB-DL اون لاين وتحميل مباشر احدث...

Hydrodynamic Simulation of Semiconductor Devices

Hydrodynamic Simulation of Semiconductor Devices

PII: Prog. Quant. Electr. Vol. 21, No. 4, pp. 293±360, 1998 # 1998 Elsevier Science Ltd. All rights reserved Printed in Great Britain S0079-6727(97)0...

3MB Sizes 0 Downloads 91 Views

Recommend Documents

Simulation of semiconductor devices with a local numerical approach
A numerical solution of the Drift-Diffusion Model for simulation of semiconductor devices based on the local meshless nu

Monte Carlo simulation of charge transport in semiconductor devices
Recent applications of Monte Carlo device simulation are reviewed, focussing mainly on GaAs structures. The results show

Polymeric semiconductor devices
We have constructed various semiconductor device structures with polyacetylene prepared by the Durham precursor route as

Fast semiconductor thermoelectric devices
The parameters of a thermoelectric radiation detector on a Bi2Te3 basis were calculated. Maximum specific detectivity re

Reliability of compound semiconductor devices
This paper reviews the reliability of III–V semiconductor devices with particular attention to the failure mechanisms ty


Prog. Quant. Electr. Vol. 21, No. 4, pp. 293±360, 1998 # 1998 Elsevier Science Ltd. All rights reserved Printed in Great Britain S0079-6727(97)00013-X 0079-6727/98 $19.00 + 0.00

HYDRODYNAMIC SIMULATION OF SEMICONDUCTOR DEVICES ARLYNN W. SMITH* and KEVIN F. BRENNAN ITT Night Vision, 7635 Plantation Road, N.W., Roanoke, VA 24019-3257, U.S.A. and School of Electrical and Computer Engineering, 777 Atlantic Dr., Georgia Tech, Atlanta, GA 30332-0250, U.S.A. AbstractÐIn this paper, we present an introduction to hydrodynamic-based simulation of semiconductor devices. A very general approach is given to illustrate a mathematical method of the solution. A general di€erential equation which is characteristic of the relevant semiconductor equations is ®rst presented and discussed. Solution of this general di€erential equation forms the basis of the mathematical model employed. A discretization equation based on the general di€erential equation is next determined. From use of Newton's method, the system of equations arising from the discretization equation applied to the nodes in the simulation domain, can then be solved. This general technique is next applied to the speci®c problem of the hydrodynamic device simulation equations. The basic equations and their auxiliary relationships governing the hydrodynamic simulation are developed. The boundary conditions and material constituent relations are presented which complete the solution. Finally, we illustrate the method with several computational examples. The workings of silicon and GaAs based ballistic diodes are examined using the drift±di€usion and hydrodynamic simulations. Additionally, we examine a GaAs/AlGaAs heterostructure bipolar transistor. Both narrow and wide base devices are simulated to illustrate the importance and relevance of the hydrodynamic method in a realistic and important semiconductor device. It is shown that the hydrodynamic simulation recovers more of the important physics governing these devices than the standard drift±di€usion method. # 1998 Elsevier Science Ltd. All rights reserved

CONTENTS 1. Introduction 2. The Boltzmann equation and its moments 3. A general di€erential equation 3.1. Creation of the discretization equations 3.2. Pro®le assumption 3.2.1. Derivation of the pro®le with constant coecients 3.2.2. General pro®le equation 3.3. Flux equation formulation 3.3.1. General ¯ux equation 3.3.2. General ¯ux equation for non-zero source terms 3.4. General discretization equation 3.4.1. General discretization equation with no source term 3.4.2. Two-dimensional discretization equation with source terms 3.5. Solution of the algebraic equationsÐNewton's method 3.5.1. Single variable Newton's method 3.5.2. Multi-variable Newton's method 3.6. General Jacobian equation 3.6.1. General partial di€erential of the one-dimensional discretization equation without source terms 3.6.2. General partial di€erential of the two-dimensional source term augmented discretization equation 3.7. Numerical issues 3.8. Node position and general boundary conditions 4. Hydrodynamic model 4.1. Background and history 4.2. Derivation of the transport model 4.2.1. Particle ¯ux equations 4.2.2. Energy ¯ux equations 4.3. System of equations for the hydrodynamic model 4.4. Boundary conditions 4.4.1. Thermal boundary conditions

*E-mails: [email protected]; [email protected] 293

295 298 301 303 304 304 305 305 306 306 307 307 308 309 309 310 312 313 315 317 319 320 320 322 323 325 328 332 334


A. W. Smith and K. F. Brennan

4.5. Electrical boundary conditions 4.6. Advanced physics 4.7. Material parameters 5. Application of the model 5.1. Silicon ballistic diode 5.2. Gallium arsenide ballistic diode 5.3. Heterojunction bipolar transistor, HBT, simulation 6. Conclusions References

NOMENCLATURE f fo kx,ky,kz x,y,z t ~u; v~g F~ext h ~ r ~x r; ~k r S(k,k') f SS G r u P L fo,fl D d n p c Te Th Tl S Dx Dy Z DZ Vkn Fluxi,j a Band Y e W m(k) t FI m K w q Nd Na C Cp ni J R Eg Rsurf F~

distribution function equilibrium distribution function momentum dimensions real space dimensions time velocity external forces Planck's constant divided by 2p real space gradient momentum space gradient scattering matrix unknown variable source and sink terms di€usion coecient density function drift velocity (convection) Peclet number [(g V L)/G] distance between two nodes unknown values at two nodes di€usivity distance between node and face electron concentration hole concentration electrostatic potential electron temperature hole temperature lattice temperature source term at node length of control volume in x direction length of control volume in y direction Bernoulli function [P/(exp (P)ÿ)] derivative of Bernoulli function unknown variable k at node n ¯ux of unknown variable between nodes i and j non-parabolicity factor variable for expected value material permittivity energy momentum dependent e€ective mass relaxation time „1 1 ej Fermi integral of ith order: ˆ G…j‡1† 0 …1‡exp…eÿp†† de Mobility Boltzmann's constant material anity electronic charge donor doping concentration acceptor doping concentration net doping concentration heat capacity intrinsic carrier concentration carrier current recombination rate bandgap recombination velocity ®eld induced force

334 335 336 337 337 342 345 358 359

Hydrodynamic simulation of semiconductor devices


1. I N T R O D U C T I O N

The emergence of relatively low cost, powerful computers has led to a fundamental restructuring in engineering and science. Traditionally, these disciplines have been founded on the careful interplay between theory and experiment. Though theory and experiment will always remain essential ingredients to the scienti®c method, the advent of numerical computing adds an additional means of probing and engineering the natural world. Numerical computing greatly extends the reach of theory and experiment by providing a powerful means of examining the behavior of non-linear systems as well as by providing a new means of observing e€ects which cannot be experimentally isolated. As a result, modern science and engineering are becoming increasingly dependent on numerical simulation and modeling to complement theoretical and experimental activities. Simulation and modeling has become important in almost all areas of engineering. As experimental, capital, and safety costs rise and computational costs continue to diminish, computer simulation of complicated processes and devices becomes far more attractive. Additionally, computer simulation provides a new means of examining processes and events which cannot ordinarily be examined through direct experimentation; the dynamics behind a large nuclear explosion for example. Accurate simulations enable an engineer to test new designs and processes rapidly, and often without incurring large capital costs. A predictive simulation enables the selection of an optimal design from among many possible choices, and provides a deeper understanding of the workings of each design. All computer simulation codes are based on two fundamental aspects. Every simulator must translate the problem into a mathematical model which accurately describes the physical situation. If the mathematics do not adequately describe the characteristics of the system then, no matter how complex the computer program is, incomplete or misleading results will be obtained. In addition, if an important property of the system is neglected then the model will not re¯ect changes in the system due to this property. Secondly, each simulator must employ a numerical technique, which when applied to the problem, provides a computationally ecient solution. Obviously, a simulation program must have an accurate mathematical model as its basis, for if the mathematical model is ¯awed then no matter how exotic the numerical technique is, a useful solution will not be obtained. The need for a robust numerical technique is just as important. A set of mathematically concise equations is useless for simulations unless a convergent and accurate numerical method is available for their computer implementation. The need for good mathematical models and numerical methods is not the only issue of importance in simulations. There is also the need for good engineering knowledge and understanding. A model misapplied can still produce results. It is up to the engineer to correctly interpret the data. A key concept in modeling is not to abandon common sense to the simulation program. The results should always be examined to ensure their compatibility with a physical understanding of the problem. In addition, solutions can be sensitive to the degree of re®nement of the numerical simulation mesh. Di€erent simulations can result in di€erent solutions for identical structures. Resolution of such con¯icting results often requires further mesh re®nement within the simulators. In the speci®c engineering problem of simulating semiconductor devices, many di€erent approaches are employed. All of these techniques are fundamentally dependent upon the solution of the basic carrier transport equations along with the Poisson equation. In the classical and semiclassical regimes, the transport is described by the Boltzmann equation. For the vast majority of semiconductor devices presently used, the Boltzmann equation provides a highly accurate description of the transport physics. Although quantum e€ects can be important in some semiconductor devices, the essential quantum physics, such as


A. W. Smith and K. F. Brennan

spatial quantization e€ects, can often be satisfactorily incorporated within a semiclassical simulator. Owing to the fact that the Boltzmann equation is an integro-di€erential equation, it is dicult to solve, particularly in most situations of interest in semiconductor device design. Analytical solutions of the Boltzmann equation can be obtained but only under extraordinarily limiting circumstances and through the application of many simplifying assumptions. As a result, numerical solutions provide the only practical means of solving the Boltzmann equation in general. Numerical solutions of the Boltzmann equation are typically provided by either its direct solution through microscopic simulation or by the method of moments. The microscopic simulation technique used for solving the Boltzmann equation is the Monte Carlo procedure. The Monte Carlo method consists of the direct simulation of the trajectory of one or more carriers in a semiconductor material or device structure subject to the action of applied electric ®elds, device geometry, and carrier, impurity or lattice scattering mechanisms. The carrier free ¯ight time and the scattering agents are selected stochastically in accordance with some given probabilities describing the microscopic physics. The Monte Carlo method, which draws its name from the stochastic aspect of the technique, provides a rigorous solution to the Boltzmann equation limited only by the extent to which the underlying physics of the problem is included. As a result, the Monte Carlo method a€ords a complete solution of the Boltzmann equation and provides a detailed description of the microscopic transport processes in a real device. Unfortunately, the Monte Carlo method is extremely computationally expensive. Additionally, there are many device applications that are not well suited to direct microscopic analysis due to both the practical and physical limitations of the model. Speci®cally, the Monte Carlo method is not well suited to studying device structures with large spatial dimensions, device structures which operate over many orders of magnitude in timescales, or devices in which thermal e€ects are important. Finally, direct comparison to experimental measurements is often frustratingly dicult to accomplish with the Monte Carlo method due to the diculty of including nonidealities in its formulation. For these reasons, other techniques are often employed in semiconductor device simulation. Old techniques such as expansion of the distribution function in spherical harmonics using orthogonal polynomials have recently undergone a renaissance and have become popular. The speed of modern computers makes it possible to discretize the energy space which is necessary for the orthogonal polynomials. The resurgence of this approach comes from the realization that the solution provides the distribution function. However, many orders of the expansion may be necessary to describe the distribution function and it has not been proven that this solution technique is any more accurate than a moment method. The most prevalent semiconductor device simulation tool is based on the coupled solution of the carrier drift±di€usion equations and the Poisson equation. The drift±di€usion equations are obtained by taking the ®rst two moments of the Boltzmann equation (the moment generating procedure is described below). The particular advantage of the drift± di€usion formulation is that it is computationally fast and requires relatively little parametrization. For these reasons, it has been used extensively in semiconductor device simulation. However, the drift±di€usion solution is highly approximate since by truncating the solution of the Boltzmann equation, it neglects many important physical e€ects such as thermal gradients, non-stationary transport e€ects, and energy dependent phenomena. These e€ects, though of little importance in large device structures and low power devices, are of increasing importance in state-of-the-art submicron feature length semiconductor devices. Therefore, usage of standard drift±di€usion simulators are becoming suspect in designing emerging device structures.

Hydrodynamic simulation of semiconductor devices


There is some need then to provide a simulation tool which is intermediate between the fast but over simpli®ed drift±di€usion simulators and the complete, but computationally demanding Monte Carlo simulators. One of the most promising approaches is the hydrodynamic simulation technique. The hydrodynamic solution is obtained by taking the next two higher moments of the Boltzmann equation over that for the drift±di€usion solution. As a result, the hydrodynamic simulation provides a practical compromise between these two methods. The hydrodynamic simulation includes more physical phenomena than the drift±di€usion model thereby extending its range of validity and accuracy, but with much more computational eciency than the Monte Carlo method. In addition, it operates at the macroscopic level thereby enabling simulation of large spatial structures, operation over many orders of timescales and inclusion of nonidealities. For these reasons, hydrodynamic simulations of semiconductor devices are becoming a standard method for designing advanced devices. The numerical methods and equations used in the simulation have now reached a level of maturity which makes their application at the engineering level possible. A great deal of engineering and mathematical e€ort has been concentrated on the solution of partial di€erential equations of the type employed to describe ¯ow processes used in device simulation. This type of problem appears in mechanical, chemical, and aeronautical engineering as well as electrical engineering. The interested reader should not ignore the wealth of work performed in these other engineering ®elds where simulation has been applied far longer than in electrical engineering. In this paper, we will discuss the basic formulation and implementation of the hydrodynamic technique for semiconductor device analysis. This work will concentrate on the two-dimensional simulation of the hydrodynamic model using the rectangular control volume formulation (also called ®nite boxes). It will not cover ®nite element methods, however a triangular control volume approach is also available with characteristics similar to those presented here(1). Each method has its advantages and disadvantages, some of the less obvious examples will be discussed in this text. This paper describes the basic numerical algorithms and methods as well as the necessary semiconductor equations for the simulation of advanced devices. After reading the paper the reader should have obtained enough knowledge to correctly apply a simulation program or create their own basic program if they so desire. This paper will not provide a comprehensive survey of the methods disclosed within the literature, nor will full mathematical proofs be given. Unlike other presentations of the hydrodynamic model for semiconductor devices this paper will approach the problem from a di€erent perspective. First, some background material on the Boltzmann equation and its moments will be presented. Next, a general di€erential equation is presented and discussed which will form the basis of the mathematical model employed. If the equations of the hydrodynamic model can be manipulated into this general form, only one numerical method would need to be programmed. A discretization equation based on the general di€erential equation is next introduced. The algorithms used to implement the discretization are then derived. The partial di€erential equations for the hydrodynamic model are then presented. A discussion of how the terms in the model can be written in the form of the general partial di€erential equation is given. Powerful enhancements of the basic model are included in the original form of the equations from the start instead of being added as an after thought. Finally some general observations are discussed.


A. W. Smith and K. F. Brennan

2. T H E B O L T Z M A N N E Q U A T I O N A N D I T S M O M E N T S

The Boltzmann equation provides the fundamental description of transport phenomena in classical and semi-classical systems. The Boltzmann equation provides a description for how the non-equilibrium distribution function evolves both spatially and temporally. In a semiclassical or classical system, the non-equilibrium distribution function can be expressed as a function of both spatial and momentum coordinates. This description is in general unsatisfactory since it violates the Heisenberg Uncertainty relationship due to the fact that it assumes that the spatial and momentum coordinates can be speci®ed simultaneously. However, for many situations present in semiconductor devices, neglect of the Heisenberg relationship is an acceptable approximation since quantum mechanical e€ects are negligible. The Boltzmann transport equation, BTE, is a valid representation of the time dependent motion of charge carriers and energy, provided the following assumptions hold: 1. The number of carriers in the device is large enough that a statistical distribution function applies(2). 2. Quantum mechanical phenomena are included phenomenologically in the e€ective mass and scattering coecient terms(2). 3. There is little carrier±carrier interaction(3). 4. The scattering probability is independent of external forces and collisions are instantaneous(3). 5. Forces applied to the electron wavepacket do not vary appreciably over the wavelength of the packet, i.e. 10 nm(4). In most devices these assumptions can be justi®ed. Only in quantum structures where quantized, distinct energy states arise will these assumptions begin to fail. The reason for this inability to model quantum phenomena is due to the semi-classical nature of the BTE; both position and momentum of the particles are assumed to be known at all times. If the regions with quantum mechanical e€ects can be considered separate and modelled by some other method (Monte Carlo techniques), then the bulk and quantum systems may be coupled as done by Widiger for HEMT's(5). The distribution function is de®ned as f (x,k,t), where x describes the spatial position, k the wavevector and t the time. Recall that k is simply related to the momentum, p, for free particles by Planck's constant. The distribution function describes the state of the system in phase space and represents the probability of ®nding a particle at a speci®ed position in phase space at a given time. Phase space is de®ned as the collective space consisting of three dimensions in real and momentum space. Therefore, a complete description of phase space requires the speci®cation of three real space coordinates, three momentum coordinates and time. The Boltzmann equation describes the evolution of the distribution function in these coordinates. The Boltzmann equation can be derived by recognizing that the probability distribution function, f (x,k,t), must be conserved. In other words, f (x,k,t) must obey a continuity equation. In the special case where there are no scatterings, the total derivative of f (x,k,t) must be zero. Hence, df ˆ0 dt


The above equation states that the total time rate of change of the probability distribution function must be zero; there is no loss of probability and hence no loss of particles; the probability ¯ux is conserved. However, in semiconductors, carriers undergo scatterings between themselves and their surroundings, in particular impurities and phonons. In

Hydrodynamic simulation of semiconductor devices


addition, carriers can move between bands through generation and recombination events. Therefore, in general the distribution changes due to scatterings and generation and recombination events. Lumping all of these e€ects together into one term, the Boltzmann equation can be written as,   df @f ˆ …2† dt @ t scatterings The above equation is called the Boltzmann equation. It can be recast however into a more convenient form by expanding out the total time derivative on the left hand side. Since f is a function of x, k, and t, the total time derivative can be expanded as, df @f @f dx @ f dy @ f dz df dkx @ f dky @f dkz ‡ ‡ ˆ ‡ ‡ ‡ ‡ dt @t @ x dt @ y dt @z dt @kx dt @ ky dt @kz dt


The above can be simpli®ed by recognizing that ~ ext dx ^ dy ^ dz ^ dkx ^ dky ^ dkz ^ eF i‡ kˆ i ‡ j ‡ k ˆ ~u and j‡ dt dt dt h dt dt dt


The force is assumed to only be external to the system. Further recognizing that, @f ^ @ f ^ @f ^ ~ @f ^ @f ^ @ f ^ ~ i‡ k ˆ rk f i ‡ j ‡ k ˆ rx f and j‡ @x @y @z @ kx @ky @kz


the Boltzmann equation can be written as,   ~ @f ~ x f  ~u ‡ r ~ k f  Fext ÿ @ f ‡r ˆ0 h @t @t scatterings


where the last term in the left hand side above includes the change in the distribution ~ ext is the function due to scatterings between particles and their surroundings, and F external force term. The scattering term in the Boltzmann equation can be de®ned in the following way,   … @f ˆ ‰ f …x,k,t†…1 ÿ f …x,k 0 ,t††S…k,k 0 † ÿ f …x,k 0 ,t†…1 ÿ f …x,k,t††S…k 0 ,k†Šdk 0 …7† @ t scatterings where (1 ÿ f(x,k,t)) is the probability of a vacancy existing in the state denoted by (x,k) at the time t. S(k,k') is the transition rate between the initial state k and the ®nal state k'. The ®rst set of terms then represents the scattering rate from the initial state x, k into the ®nal state x, k' at the time t assuming a coupling expressed by S(k,k') taking into account that the initial state must be occupied and the ®nal state empty in order that the transition can occur. Similarly, the second term represents the scattering rate from the state x, k' into the state x, k at time t. Therefore, the net rate at which the state x, k changes with time due to scatterings is given by the di€erence between the outgoing and incoming scattering rates. The identity of the term S(k,k') depends upon the physical mechanism which causes the scattering interaction. In semiconductors, the most important scattering mechanisms are carrier-phonons, carrier-impurities, and carrier±carrier interactions. Additionally, carriers can be generated or recombined through radiative, thermal or Auger interactions. These mechanisms also contribute to the collective scattering term. The other terms in the Boltzmann equation can be understood as follows. The ®rst term on the right hand side is the time rate of change of the probability distribution function. The


A. W. Smith and K. F. Brennan

second term on the right hand side represents the ¯ux of the distribution in the spatial dimension while the third term represents the ¯ux of the distribution in momentum space. As discussed in the introduction, the solution of the Boltzmann equation is in general quite dicult owing primarily to the fact that it is an integro-di€erential equation. Some analytical solutions of the Boltzmann equation have been provided(6±8), however in general numerical techniques must be employed for most problems. As mentioned above, the Boltzmann equation can be solved directly by microscopically simulating the trajectories of one or more carriers in phase space subject to the interactions with all scattering agents, applied potentials and boundaries. The most successful such microscopic approach is the Monte Carlo method. However, owing to the diculties encountered in the use of Monte Carlo simulations for large device structures, structures which operate over long timescales and structures which are dominated by imperfections, other techniques are typically employed. Most commonly used semiconductor device simulators use either the drift±di€usion or hydrodynamic solutions of the BTE(9±11). In either case, the basic equations which govern these approaches are based on the method of moments solution of the BTE. Moments of the BTE are obtained by multiplying Eqn. (6) by the quantity to be conserved and integrating over all k space. In this manner, the expected value of the parameter is determined. To completely describe the system, an in®nite set of moments must be calculated because each moment will contain a moment of the next higher order. Therefore, at some point an approximation must be made for the highest moment. In the derivation to follow, the ®rst moment is the continuity equation (the quantity to be conserved is the particle number). It has been assumed that the highest moment is the ®fth moment which can be approximated as an energy squared term. The second moment provides the expected value of the momentum, the third moment the expected value of the energy, and the fourth moment provides the expected value of energy ¯ux. Below, the resulting equations for each of these moments will be presented. The primary diculty in solving the BTE through the method of moments is the treatment of the scattering term. Generally, the exact nature of each scattering mechanism must be determined and an averaging over that mechanism must be made. However, a simpli®cation can be made to the scattering term which greatly reduces the mathematical complexity, yet retains much of the important physics. The following observations can be made about the scattering term which can be utilized in its simpli®cation: 1. The scattering term must vanish in equilibrium (scattering still occurs but the term must evaluate to zero). 2. When an external stimuli disturbs the distribution, scattering ultimately returns the system to equilibrium once the action of the external stimuli is removed. In other words, if the system is initially driven into non-equilibrium by some external driving force, then when the external forces are all removed, the system must return to equilibrium. The scattering mechanisms are responsible for restoring the system back to equilibrium. (The only exception to this rule is a superconductor, which does not relax back to equilibrium when all external driving forces are removed). The simplest form for the scattering term that is consistent with these two criteria is:   @f ÿ… f ÿ fo † …8† ˆ @t c t The quantity t is the relaxation time, and fo and f are the equilibrium and non-equilibrium distribution functions respectively. Notice that the complicated integral expression for the scattering term has been replaced by a much simpler expression. The basic assumption

Hydrodynamic simulation of semiconductor devices


made in Eqn. (8) is called the relaxation time approximation. The above equation implies that the system will relax to the equilibrium distribution, fo, in a characteristic time speci®ed by the relaxation time, t. The relaxation time physically represents the mean time, in which, through the action of scatterings, the system relaxes from an initial nonequilibrium state into a ®nal equilibrium state. Substitution of this formulation into Eqn. (6) produces the ®nal form of the Boltzmann transport equation that will be used in the derivations presented in this chapter,   @f ~ ‡ F~  r ~ k f ˆ ÿ… f ÿ fo † …9† ‡ ~u  rf @t t ~ ~ cÿrW) for electrons(12). To simplify where the ®eld induced general force F~ is given by (ÿre the derivation, the crystal structure is assumed to be cubic with ellipsoidal energy surfaces. This reduces the e€ective mass tensor to a vector and the temperature tensor into a scalar(2). In all the cases to be described in this chapter, moments can be generated for holes by a similar derivation. If the semiconductor has satellite valleys then moment equations for each valley must be formulated. The equations for di€erent valleys would then be coupled by the scattering terms. Unlike the drift±di€usion model which takes only the ®rst two moments, the hydrodynamic model is created by taking the ®rst four moments of the Boltzmann transport equation (BTE). Therefore, the drift±di€usion model will require an approximation for the third moment to form a closed system of equations. Since the hydrodynamic model contains the ®rst four moments, the ®fth moment closes the system, it should be obvious that by using suitable approximations, the drift±di€usion model may be recovered. Before, we present the derivation of the drift±di€usion and hydrodynamic models using the moment method, we will present a more general formalism for numerically solving any general ¯ow equation. Based on this formalism, one can develop a numerical solution to the Boltzmann equation to any moment order. It should be noted that one of the most important problems in solving the BTE through the use of higher order moment equations is that the resulting system of equations, particularly for solutions involving more than the ®rst two moments, i.e., those solutions that go beyond the drift±di€usion formulation, are numerically often very dicult to solve. This is due to non-linearities in the equations as well as the fact that the carrier concentrations can change by many orders of magnitude throughout the simulation in either time or space(13, 14). The solution of even the drift±di€usion equations coupled with the Poisson equation requires a special formulation to ensure numerical stability. The ®rst such formulation in electrical engineering literature was presented by Scharfetter and Gummel(13). Their approach has been widely used in most semiconductor device simulations because it allows for relatively large changes in the potential and carrier concentrations between neighboring mesh points without causing numerical instability. In the following sections, we present a generalization of the technique developed by Scharfetter and Gummel(13), which can be used in solving higher moment solutions of the BTE. It should be further noted however that the method presented herein is not limited solely to solutions of the BTE. It is suciently general that it can be applied to any system in which the ¯ow of some quantity is to be studied. 3. A G E N E R A L D I F F E R E N T I A L E Q U A T I O N

In this chapter, we present the details for formulating a numerical solution of the basic semiconductor transport equations based on a more general formalism. The particular advantage of this scheme is that it can be applied to any general ¯ux equation, therefore its


A. W. Smith and K. F. Brennan

extension to any of the moment equations useful in semiconductor device modeling is immediate, though not necessarily simple. The approach taken is often referred to as the control volume technique(1). The control volume technique is somewhat intermediate between the more commonly used ®nite di€erence and ®nite element methods, however the ®nite di€erence scheme can be recovered from the control volume formulation under some assumptions. In general the e€ect which is to be modelled is the ¯ow of an unknown quantity within some volume. What is desired is a knowledge of how much of the quantity is contained in a volume in space as a function of time. The physical quantity can ¯ow out of the volume into a surrounding volume or be generated or lost within the interaction volume itself. This description applies to electrons\holes ¯owing and recombining in a semiconductor, as well as chemicals ¯owing and reacting in a vessel. A general conservation equation which possesses this structure for the unknown variable f can be written as, …time rate of change in f† ‡ …Change in flux of f† ˆ …Source ÿ Sink of f†


which can be expressed mathematically in the following form @…rf† ~ ÿÿÿÿ4 @ …rf† ~ ~ ˆ SS ‡ r  Flux ˆ SS or ‡ r  …r~ uf ÿ Grf† @t @t


In the second form of the equation, the expression for the ¯ux has been substituted. In special cases, the above equations can contain terms which may be non-linear or even zero. The numerical method should apply equally as well and as smoothly to these cases as well as the full equation. The ¯ux of the unknown parameter f includes both a convection ~ (drift) (r~uf) and di€usion component (Grf). The parameter r is the density of the unknown parameter (f). The ¯ow velocity of f is u, and G is the di€usion coecient. The source and sink terms (generation\recombination) are depicted phenomenologically by the symbol SS. The general di€erential is a vector equation in terms of the ¯ow velocity and the di€usion term. In the above context, the word di€usion is meant in a generalized sense and can represent any di€usion ¯ux due to the gradient of the general variable f. The conservation equation alone does not specify the entire problem. As is always the case in solving di€erential equations, the boundary conditions must be speci®ed. Boundary conditions must be speci®ed for the unknown parameters or the ¯ux at the edge of the simulation region forming the boundary with the external environment. In the hydrodynamic simulation six of these equations are required, one for each of the variables. General analytical solutions to these equations may not be available, and if they are available these solutions may contain in®nite series or other special functions. What is required is a method of solving the equation exactly over a small volume where the parameters are not changing rapidly. By solving the di€erential equation in this small volume with information from the surrounding volumes as boundary information the total solution over the entire simulation domain can be constructed. The outer boundary conditions of the simulation domain occur only at the outer small volumes of the domain. Our approach can be outlined as follows. Given the general equation that we seek to solve, the ®rst step in its numerical solution is to select a discretization mesh on which the equation is to be solved. Next a pro®le for the parameter of interest is chosen. The pro®le is the interpolation scheme for the value of the unknown de®ned between any two points in the discretization mesh. From the pro®le a general ¯ux equation can be generated. A similar result is found for the general ¯ux equation as compared to that of the Scharfetter± Gummel formulation. The general ¯ux equation is then used to create a discretization equation for the general di€erential equation. This converts the di€erential equation into an algebraic system of equations which relate the value of the unknown at one node to the

Hydrodynamic simulation of semiconductor devices


values of the unknown at the surrounding nodes. Once the di€erential equation has been written in terms of a system of algebraic equations these equations are solved using Newton's method. The system of equations are often quite extensive. For example, for the hydrodynamic model discussed below, six equations are needed. If the simplest stencil is employed, where only the four nearest neighbor points to each node are used to discretize the equations, then the number of linear algebraic equations in the system is 30 n  30 n for a two-dimensional solution where n is the number of nodes in the domain. To solve such a large system various numerical techniques can be employed. We will present a short discussion of some of the options one can select. Below we discuss in some detail each of the steps in the solution. We begin with a short discussion of the discretization scheme. 3.1. Creation of the discretization equations Consider the volume shown schematically in Fig. 1A. If the parameters which describe the general unknown are to be found, then the di€erential equations must be solved over the entire domain. Now consider the same volume but with it broken down into small discrete blocks, as shown in Fig. 1B. With suitable approximations, the di€erential equations can be solved on these discrete blocks with boundary conditions from the surrounding elements. For instance, if it is acceptable to ®nd the values of the parameters at discrete points (at the center of the blocks) in the volume, then a suitable interpolation scheme may be constructed to estimate the parameter between the nodal points. The values of the parameters at the nodal points can be found by solving the di€erential equations on an area encompassed around the node, using information from the surrounding nodes as boundary data. This leads to sets of algebraic equations which relate the values of the unknown parameter of the general di€erential equation at one node to the values of the parameter at the surrounding nodes. Thus the system of di€erential equations can be transformed to a system of algebraic equations (discretized equations) which may be solved by numerical techniques. With the advent of large digital computers and numerical methods, these tools hold the promise that mathematical models may be created to determine all of the parameters in the simulation domain for the most complex physical problems.

FIG. 1. Device simulation domain (A) in physical space; (B) discretized domain.


A. W. Smith and K. F. Brennan

3.2. Pro®le assumption The ®rst step in creating a useful program for solving the general di€erential equation is to understand how the forces acting on the unknown quantity can a€ect the pro®le of the unknown between two discrete points. The pro®le is de®ned as how the variation of the unknown between any two discrete mesh points can be modeled. In other words, speci®cation of the pro®le can be interpreted as de®ning the interpolation scheme of the parameter between the two nodal points. The pro®le is determined by the boundary conditions, the form of the di€erential equation, and the forces which are acting on the unknown quantity. The simple case of a one-dimensional steady state ¯ow with constant forcing terms (r, u, and G constant), and no source terms will be examined ®rst, then the more complicated general case is examined. 3.2.1. Derivation of the pro®le with constant coecients. The pro®le for the unknown parameter (f) in the case of constant coecients (r, u, G) will now be presented assuming the steady state situation with no source or sink terms (SS = 0). The node cluster is shown in Fig. 2. First, the general di€erential Eqn. (11), after integrating once, is rewritten in the following form, ruf ÿ G

@f ˆ Flux @x


Setting the Flux = 0 (for the zero ¯ow condition) and re-arranging the equation one obtains, ruf ˆ G

@f @x


The solution for the zero ¯ow pro®le will be solved using this equation. It will then be assumed that the non-equilibrium pro®le does not deviate substantially from the equilibrium pro®le. Separating the variables results in the following equation, ru df dx ˆ G f


At this point, the above di€erential equation can be readily solved with the following assumptions and by specifying certain boundary conditions as,

FIG. 2. One-dimensional node cluster for the de®nition of the ¯ux terms.

Hydrodynamic simulation of semiconductor devices

i. ii. iii. iv. v.


ru is constant G is constant at x = 0 f = f0 at x = L f = fL de®ne P = ruL/G and therefore (P/L)dx = df/f

P is called the Peclet number and will be used throughout this chapter. It provides a measure of the strength of convection (drift) compared to the amount of di€usion in the system. With the assumptions and conditions speci®ed above, the resulting pro®le can be obtained by solving for the parameter from x = 0 to x and from x = 0 to x = L. The pro®le is given as the ratio of the di€erence in the parameter over these two ranges of x. Performing the integrations, the pro®le becomes(1),   3 P x ÿ1 exp 6 7 f ÿ fo L 7  ˆ6  5 P fL ÿ fo 4 L ÿ1 exp L 2


Notice that at x = 0 the value of f is fo and at x = L the value is fL as speci®ed by the boundary conditions. The above analysis is suitable for constant coecients for both the di€usion and convection (drift) components. However, in the context of the hydrodynamic equations this assumption may no longer be valid. Therefore, a need exists for a more general pro®le which will now be presented. 3.2.2. General pro®le equation. If the di€usion coecient and convection components are integrable over the discretization volume then these parameters may be substituted into the derivation above. The general pro®le equation using these non-constant parameters from the general di€erential equation (Eqn. (11)) is then 2

3  ru d x ÿ 17 6 exp f ÿ fo G 6 7 ˆ 6 … 0L  7 5 fL ÿ fo 4 ru dx ÿ 1 exp 0 G … x


A comparison of Eqn. (16) with Eqn. (15) con®rms that these two equations will result in the same type of pro®les. This equation gives the one-dimensional interpolated pro®le of the unknown variable between two locations for the case of no source terms. 3.3. Flux equation formulation Recall the general di€erential equation previously described (Eqn. (11)) @…rf† ~ ÿÿÿÿ4 ‡ r  Flux …f† ˆ SS @t


From inspection of Eqn. (17) it is obvious that the ¯ux of the unknown is required, not the pro®le of the unknown. Now given that an equation for the pro®le of the unknown variable f between two nodes has been formulated, the next step in the analysis is creation of a ¯ux equation between two nodes. This is accomplished by substituting the exact pro®le for the unknown and the derivative of the pro®le into the appropriate ¯ux equation and simplifying the result.


A. W. Smith and K. F. Brennan

3.3.1. General ¯ux equation. Both Eqns (15) and (16) can be substituted into the current ¯ow Eqn. (12) and can be written in the following form   fi exp…P † ÿ fj …18† Flux…f†i, j ˆ D exp…P † ÿ 1 which is the ¯ux of the unknown variable between two nodes i and j and D is the di€usion conductance (D = G/L). This general form is similar to the Scharfetter±Gummel formulation(13). In other engineering literature(7) this form is called the exact formulation or exponential form. The fact that a simple solution exists for the ¯ux term is important because it provides a convenient general formulation which can be applied to all six of the unknowns (n, c, p, Te, Tl, and Th) that appear within the hydrodynamic model. 3.3.2. General ¯ux equation for non-zero source terms. The above discussion for the ¯ux equations applies only for di€erential equations in which the source term is set to zero (SS = 0). Clearly, in the general case of systems of equations, the source terms in all of the equations of the model will not necessarily be zero. Since the goal of this paper is to generate a general two-dimensional model, derivation of a ¯ux equation for the twodimensional non-zero source term case is required. The general ¯ux equation for this case is given by 2 0 13 P   ÿ 1 Bexp…P † ÿ 1 6 fi exp…P † ÿ fj d C7 7 ‡ C ‡ Si Li B …19† Flux…f†i,j ˆ DY 6 @ 4D P exp…P † ÿ 1 Li A5 The spatial parameters are de®ned in Fig. 3. The ¯ux is de®ned for ¯ow between any of the outer nodes and the central node separated by distance Li. The parameters D and Y are due to the second spatial direction. For the one-dimensional case the terms Y and D both become unity (1). The source term associated with the node face is Si, the harmonic average

FIG. 3. Two-dimensional node cluster for the de®nition of the ¯ux terms with and without source terms in the formulation.

Hydrodynamic simulation of semiconductor devices

of the values of the source terms at the node points. Si is given by  ÿ1 1 ÿ fi fi ‡ Si ˆ Sp S1



where SI is the source term of the surrounding node (I $ [E, W, N, S]). Comparison of Eqns (18) and (19) reveals that they have the same dependence on D; the factor multiplied by D in both cases is identical. Notice that when the source term is included as an additional term, the original term in the zero source formulation is unchanged. Examination of Eqn. (19) reveals that the ¯ux in this case will be position dependent; i.e. the ¯ux is linear in position due to the source term component. Notice that unlike Eqn. (18), this formulation will require more extensive programming complexity in order to be incorporated into existing models. Nevertheless, some authors have included the non-zero source term in their models. A one-dimensional method which accounts for non-zero source terms is contained in PC-1D(15), and is also derived by a di€erent method in Miller(16, 17). Comparison of the one-dimensional ¯ux equations from this work and those in references(15±17) shows good agreement. At this point, a discussion of grids and discretization schemes is in order. The ¯ux equations presented so far are general in nature, applying to ¯ow between any two points given the parameters at the points. In rectangular control volumes the ¯ow is always perpendicular to the faces of the volume. In this way the formulation is exact for ¯ow in one direction. In triangular control volumes the ¯ux equations are again for ¯ow perpendicular to the triangle sides. However, it will not give the exact solution for simple one-dimensional ¯ow due to the error of the triangles not aligning to the ¯ow ®eld. This is not to say that the rectangular volumes are always superior. For example, for ¯ow which is at 45 degrees to a rectangular volume an error is introduced due to the ¯ow angle, which is referred to as false di€usion. If triangular control volumes are used instead, false di€usion will not occur. Several attempts have been made to minimize false di€usion for a rectangular control volume(18). In general, the following tradeo€ exists. The rectangular control volume is numerically exact for one-dimensional ¯ows but introduces false di€usion for angular ¯ows, while triangular volumes minimize the false di€usion e€ect but are not numerically exact for one-dimensional ¯ows. 3.4. General discretization equation At this point the general ¯ux equation will be used to create the ®nal discretization equation for the general di€erential equation. This converts the general di€erential equation into an algebraic system of equations which relate the value of the unknown at one node to the values of the unknown at the surrounding nodes. The use of the general discretization equation facilitates programming due to the fact that only one general equation is programmed, di€erent coecients are calculated depending upon which of the six unknowns is to be solved. 3.4.1. General discretization equation with no source term. The steady state general di€erential equation in two dimensions with no source terms can be written as, (see Fig. 3 for the associated node cluster), ÿÿÿÿ 4

~  Flux ˆ 0 r @ Fluxx @ Fluxy ‡ ˆ0 @x @y


Substituting the values of the ¯uxes at the node faces into this equation and applying a central di€erence for the partials yields the following equation


A. W. Smith and K. F. Brennan

Fluxe ÿ Fluxw Fluxn ÿ Fluxs ‡ ˆ0 Dx Dy


The general discretization equation for the case of no source term is then obtained by substituting Eqn. (18) into Eqn. (22) and gives the ®nal form ÿÿÿÿ 4

~  Flux †DxDy ˆ 0 …r ÿ ÿ  ÿ ÿ   ÿ  ÿ  1Dy De fE Z ÿ Pe ÿ fP Z Pe ÿ Dw fP Z ÿ Pw ÿ fW Z Pw ÿ ÿ  ÿ ÿ   ÿ   ÿ ‡Dx Dn fN Z ÿ Pn ÿ fP Z Pn ÿ Ds fP Z ÿ Ps ÿ fP Z Ps


Where Z(P) is the Bernoulli function {P/(exp(P) ÿ 1)}, P is the Peclet number associated with the node face, and D is the di€usivity associated with the node face. Eqn. (23) relates the value of the unknown parameter f at the point P to the values of the parameter f at the four surrounding nodes E, W, N, and S. 3.4.2. Two-Dimensional discretization equation with source terms. The two-dimensional discretization equation including the source terms, derived in the same manner as that above, is given by    ÿÿÿÿ 4  ÿ ~  Flux ÿ SP DxDy1DyY Fluxo ÿ Fluxo ‡ Se Le C‡ …Pe † ‡ de r e w Le      ds dn ‡ DxX Fluxon ÿ Fluxos ‡ Sn Ln C‡ …Pn † ‡ ÿ Sw Lw Cÿ …Pw † ÿ Lw Ln   ds …24† ÿ SP DxDy…X ‡ Y ÿ XY † ˆ 0 ÿ Ss Ls Cÿ …Ps † ÿ Ls where the following de®nitions have been made:  ÿ Fluxoi ˆ Di fi Z…ÿPi † ÿ fj Z…Pi †


SP is the source term at node P. Se,w,n,s are the harmonic averages of the source terms at the node faces and   1 0  dn ds Ln C‡ …Pn † ‡ ÿ Ls Cÿ …Ps † ÿ B Ln Ls C C Yˆ1ÿB @ A Dy    1 de dw Le C‡ …Pe † ‡ ÿ Lw Cÿ …Pw † ÿ B Le Lw C C Xˆ1ÿB @ A Dx 0


where the terms C+ and Cÿ are de®ned by C‡ …Pi † ˆ Cÿ …Pi † ˆ

Z…Pi † ÿ 1 Pi

Z…ÿPi † ÿ 1 Pi


Hydrodynamic simulation of semiconductor devices


Notice that as the source term (S) goes to zero the two forms (with and without the source terms) of the discretization equation approach the same value, as is expected. As implied in reference(17), for the one dimensional case, this form should provide a stable numerical discretization scheme for regions of very strong source or sink terms or in regions of very strong convection (drift) with an associated source term. However, inclusion of the source term can lead to numerical instabilities if not properly damped. We have found that inclusion of the source term in a twodimensional simulation often produces numerical instabilities particularly if the source terms vary greatly between adjacent volumes.

3.5. Solution of the algebraic equationsÐNewton's method Now that a general algebraic equation has been derived the next step is the development of a technique for solving the equations. The mathematical models for the discretization equations of the unknown variables were introduced in the previous sections. The next step is to use these equations to ®nd the values of the unknown variable at the node locations in the discretized domain. The equations encountered in most applications are highly non-linear. Newton's method is a commonly used technique for ®nding the roots of a non-linear equation. Newton's method for one variable exhibits quadratic convergence to the root if the initial guess is relatively close to the real root. In solving a system of multiple equations which are coupled through some term, a multi-variable Newton's method is required. In the multi-variable formulation the convergence rate is not guaranteed to be the same as the single variable method. It is hoped that by using a previous solution of the parameters as the initial guess super-linear convergence can be achieved. Applying the multi-variable Newton's method to the system of equations for the model creates a linearized Jacobian matrix. The elements of the Jacobian are created by a general partial di€erential of the general discretization equation. In this and the following sections the partial di€erential and its associated partial derivatives are introduced. First we present the details of the single variable Newton's method. The method of solution for the linear set of equations required for the multivariable Newton's method is discussed next. In Section 3.6 the numerical approximations for the implementation of certain functions in the discretization equation and the general Jacobian are discussed. 3.5.1. Single variable Newton's method. Before discussing the multi-variable Newton's method, the single variable Newton's method is discussed because it forms the basis for the multi-variable method and is used in other sections of the hydrodynamic model. Newton's method is a powerful technique to solve for the roots of a non-linear equation. The method relies on a Taylor series expansion of the function. For instance, the problem may be described as(19) find x * such that F…x * † ˆ 0


Using a Taylor series expansion of the function F about some point x in the neighborhood of x* and de®ning F'(x) as the ®rst derivative of F(x) obtains, F…x * † ˆ 0 ˆ F…x† ‡ F 0 …x†Dx where Dx ˆ x * ÿ x



A. W. Smith and K. F. Brennan

Solving this equation for the value of x* gives x* ˆ x ÿ

F…x† F 0 …x†


Given an initial guess of x, Eqn. (30) provides a better estimate of x, namely x*, closer to the root of the equation F(x). By repeatedly substituting the value of x* in for x the equation may be iterated until convergence to the root is obtained (F(x*) = 0), provided that the derivative exists and is not equal to zero. This same concept may be applied to equations which contain many variables as discussed below. 3.5.2. Multi-Variable Newton's method. The multi-dimensional Newton method starts from the same point as the one-dimensional case, except now a system of equations (fi, 1 Ri Rn) which are functions of several variables (xi, 1 Ri Rn) must be solved. De®ne the parameter Zi as the solution (or a better estimate of the solution) and x io as the initial guess. First, the Taylor series expansion is written, neglecting second order terms as(19), @fi …x 1o ,:::,x no † @x 1 1 n 1 n @ fi …x o ,:::,x o † n n @ fi …x o ,:::,x o † ‡ …Z2 ÿ x 2o † ::: ‡ …Z ÿ x † o @x2 @x n

0 ˆ fi …Z1 ,:::,Zn †1fi …x 1o ,:::,x no † ‡ …Z1 ÿ x 1o †


A similar equation is obtained for all values of i (i $ [1, . . . n] to cover all of the parameters in the model). Consequently, a system of n equations in n unknowns is obtained. If we specialize the system to six equations in six unknowns, which is the situation for the hydrodynamic model, the resulting system of equations can be expressed in a vector matrix form as, 3 2 13 2 ÿf1 …x o † Dx 7 7 6 6 7 7 6 6 7 7 6 6 7 7 6 6 7 7 6 36 2 7 @f1 @ f1 @f1 @f1 @ f1 @f1 6 2 7 6 6 Dx 7 6 ÿf2 …x o † 7 7 7 6 6 @x 1 @ x 2 @x 3 @x 4 @ x 5 @x 6 76 7 7 6 76 6 7 6 7 7 6 @f 6 7 7 6 6 2 @ f2 @f2 @f2 @ f2 @f2 76 7 6 7 7 6 1 6 @ x 2 @x 3 @x 4 @ x 5 @x 6 76 7 7 6 6 @x 7 76 3 7 6 6 6 @f3 @ f3 @f3 @f3 @ f3 @f3 76 Dx 7 6 ÿf3 …x o † 7 7 7 6 76 6 7 7 6 6 @x 1 @ x 2 @x 3 @x 4 @ x 5 @x 6 76 7 7ˆ6 76 6 …32† 7 6 7 7 6 @f 6 7 7 6 6 4 @ f4 @f4 @f4 @ f4 @f4 76 7 7 6 76 6 1 @ x 2 @x 3 @x 4 @ x 5 @x 6 76 7 7 6 6 @x 76 Dx 4 7 6 ÿf4 …x o † 7 6 7 7 6 6 @f5 @ f5 @f5 @f5 @ f5 @f5 76 7 7 6 76 6 7 7 6 6 @x 1 @ x 2 @x 3 @x 4 @ x 5 @x 6 76 7 7 6 76 6 7 6 7 5 4 @f 6 @ f6 @f6 @f6 @ f6 @f6 6 6 7 7 6 7 7 6 6 1 2 3 4 5 6 @x @x @x @x @x @x 6 Dx 5 7 6 ÿf5 …x o † 7 7 7 6 6 7 7 6 6 7 7 6 6 5 5 4 4 Dx 6

ÿf6 …x o †

Where Dxi is the di€erence between the guess and the true root (Ziÿx io ). The matrix containing all of the partials with respect to the variables is called the Jacobian (J) of the system of equations. Provided that the determinant of the Jacobian is non-zero it

Hydrodynamic simulation of semiconductor devices


may be inverted, and an iterative scheme similar to the one-dimensional case is possible. A shorthand version of this system of equations can be written as, JDx ˆ ÿF


where J is the Jacobian matrix, Dx is a column matrix containing the di€erence between the guess and the root (Ziÿxi), and F is a column matrix containing the functions (fi) evaluated at the previous solution. Next the system of di€erential equations for the general equation in one dimension has to be put in the proper form for the application of Newton's method. This is accomplished by the use of the general discretization equation introduced previously. A one-dimensional discretization equation with a simple source term will be used as an example of this method. The general di€erential equation and discretization equation for one of the unknown variables is dFlux…fi † ÿ S…fi † ˆ 0 dx


fi ˆ Fluxoe …fi † ÿ Fluxow …fi † ÿ DxS…fi † ˆ 0


Recall from the presentation of the general ¯ux Eqn. (18) that it depended upon the value of fi (the unknown) at two nodes. Two ¯ux expressions are present in the simple discretization equation above (Eqn. (35)). These have one node in common, the central node (fiP ), and two other nodes, the node to the east and to the west. For this example, assume that the source term (S) only depends upon the values of the unknown at the central node location. The values of the parameter fi at the other two nodes are also unknown (fiE and fiW ). Therefore, these values must also be included in the Jacobian. The expanded form of Eqn. (35) is   ÿ ÿ fi …fiP † ˆ Die fiP Z…ÿPie † ÿ fiE Z…Pie † ÿ Diw fiW Z…ÿPiw † ÿ fiP Z…Piw † ÿ DxS…fiP † ˆ 0


The above equation expresses the value of the ¯ux of the unknown at the node P in terms of its values at the nodes E, W and P. Forming the di€erential of Eqn. (36) with respect to the unknown at the three nodes results in the following equations,  ÿ  @fi …fiP † @ Die fiP Z…ÿPie † ÿ fiE Z…Pie † ˆ @ fiE @ fiE    ÿ ÿ @fi …fiP † @ Die fiP Z…ÿPie † ÿ fiE Z…Pie † ÿ Diw fiW Z…ÿPiw † ÿ fiP Z…Piw † ˆ @fiP @ fiP @S…fiP † @fiP  i ÿ i  ÿ@ DW fW Z…ÿPiw † ÿ fiP Z…Piw † ÿ Dx

@fi …fiP † ˆ @fiW

@ fiW


Similar equations can be readily found for the ¯ux of the unknown at nodes E and W, fi (fiE ) and fi (fiW ), which depend upon the central node (P) and only one surrounding node. Arranging these equations into a vector matrix form provides insight into how the values of the parameters at the nodes will interact. Substituting the partials into a matrix similar in form to Eqn. (32) obtains,


A. W. Smith and K. F. Brennan



@ fi fiE 6 @ fi 6 E

6 ÿ  6 @ f fi 6 i P 6 6 @ fiE 6 6 4 0

ÿ  @fi fiE @ fiP ÿ  @ fi fiP

@ fiP ÿ  @fi fiW @ fiP

3 2 3 2 i fi …fiE † 3 DfE 7 6 7 6 0 7 6 7 76 7 6 7 76 7 6 7 6 ÿ i7 7 6 7 6 7 7 6 7 @ fi fP 76 i 7 i 7 6 6 76 DfP 7 ˆ ÿ6 fi …fP † 7 i @ fW 7 7 6 7 76 7 6 7 ÿ i  76 7 6 7 6 @fi fW 56 7 6 7 5 4 5 4 @ fiW i i DfW fi …fW †


This example is for one equation at three nodes. In the approach taken here, only the nearest neighbor nodes are considered in the analysis. Therefore, the equations are solved at only three nodes in one dimension. In most cases multiple equations must be solved. If one considers again only the nearest neighbors, then in general, multiple equations must be solved at the same three nodes. Each block contains the partial derivatives of all n equations with respect to all n of the unknowns at the node being examined, i.e., 8 i 2 ‰1,nŠ > ÿ i > < @fi fk j 2 ‰1,nŠ …39† where k 2 fE,P,Wg > @ f jl > : l 2 fE,P,Wg As an example, consider the case for the hydrodynamic model in which six unknowns are present. The six equations describing these unknowns must be solved at each of the three nodes. In this case, each entry in the Jacobian matrix will contain a 6  6 block, making the size of the total Jacobian 18  18. It is important to realize that in the discretization of a device in one dimension many more than three nodes are required in general to cover the device domain. A system of similar equations is created for each node in the simulation domain. It is clear that as the number of nodes increase the size of the Jacobian matrix increases rapidly. The elements DfiE ; DfiP ; DfiW ; fi …fiE †; fi …fiP †; and fi …fiW † in Eqn. (38) are all vectors with a length equal to the number of variables in the problem. It should also be obvious that when the simulation domain is extended to two dimensions the size of the Jacobian matrix increases dramatically. In the two-dimensional simulation domain the values of the parameters will be dependent upon ®ve node locations if again only nearest neighbor interactions are considered. The extension to ®ve nodes follows readily from the fact that one needs now to consider solutions along two di€erent directions. Using the same example as above, that of a six equation hydrodynamic model, the size of each entry within the Jacobian expands from 18  18 to 30  30. A corresponding increase in the system Jacobian is also observed.

3.6. General Jacobian equation In the previous section it was shown that the partial derivatives of each of the model equations with respect to each variable needs to be solved at all the nodes in order to perform Newton's method. Therefore, a general and ecient method of calculating these partials is required. Here again, the power of selecting a general di€erential equation and, thereby, creating a general discretization equation will provide a concise numerical method.

Hydrodynamic simulation of semiconductor devices


As before, routines to determine only certain parameters associated with the speci®c equation will be required. The general method will be the same for all equations. The discussion for the form of the general partial di€erential will begin with the simple case of the one-dimensional zero source term discretization equation. Then the partial of the twodimensional source term augmented discretization equation will be explored. 3.6.1. General partial di€erential of the one-dimensional discretization equation without source terms. To create the Jacobian matrix, the partial derivatives of the general discretization equation (Eqn. (23)) are required. The creation of these partial derivatives is complex and will be illustrated in a number of steps, with each step adding complexity. This section demonstrates the procedure to obtain the partials for the one-dimensional discretization equation without source terms. The ®rst step is to begin with the general discretization equation for the ¯ux of the unknown between two node points without source terms, Eqn. (18), ÿ ÿ  ÿ  Flux…f†oi, j ˆ Di, j fi Z ÿ Pi, j ÿ fj Z Pi, j


Where Di, j Pi, j fi fj Z

is is is is is

the the the the the

di€usion coecient between nodes i and j Peclet number between nodes i and j unknown parameter at node i unknown parameter at node j Bernoulli function

Notice that the ¯ux has been de®ned as Fluxoi,j . This de®nition becomes important as terms in this form will be used in the next section to provide a shorthand notation for the source term discretization. To form the Jacobian matrix requires the partial derivative of Eqn. (40) at the nodes associated with the ¯ux. Therefore, de®ne Vk as one of the variables to be solved. The Jacobian for variable Vkn at node n (n $ [i, j], k $ [1, n] for one of the model equations) is given by @Flux…f†oi,j @Vkn


 ÿ ÿ ÿ   @ Di,j fi Z ÿ Pi,j ÿ fj Z Pi,j @ Vkn


which becomes @Flux…f†oi,j @ Vkn


ÿ ÿ   @Di,j ÿ fi Z ÿ Pi,j ÿ fj Z Pi,j @Vkn ‡ Di,j

ÿ ÿ !   @ fj ÿ  @Z ÿ Pi,j @ Z Pi,j @fi ÿ Z ÿ Pi,j ÿ k Z Pi,j ‡ fi ÿ fj @Vkn @ Vkn @Vkn @ Vn


This expansion is a result of applying the product rule of calculus. The next step is to evaluate the partials @Z…ÿPi,j † @Z…Pi,j † and k @Vn @Vkn


Substituting in the form for the Bernoulli function and carrying out the proper steps, the


A. W. Smith and K. F. Brennan

partial derivative for the positive argument function may be written as   @Z…Pi,j † @Z…Pi,j † 1 ÿ Z…ÿPi,j † ˆ Z…Pi,j † Pi,j @ Vkn @Vkn


A similar analysis for the negative argument results in   @ Z…ÿPi,j † @Z…Pi,j † 1 ÿ Z…Pi,j † ˆ Z…ÿP † i,j Pi,j @ Vkn @Vkn


Notice that only the partial derivative of the Peclet number between the two nodes will need to be calculated. To simplify writing of later equations the following de®nitions are made,     1 ÿ Z…Pi,j † 1 ÿ Z…Pi,j † and DZ…ÿPi,j † ˆ Z…ÿPi,j † …46† DZ…Pi,j † ˆ Z…Pi,j † Pi,j Pi,j The numerical implementation of the functions in Eqn. (46) will be explored further in a later section of this chapter. Making the substitutions called for by Eqns (44)±(46) into Eqn. (42) will produce the ®nal form of the partial di€erential equation @Flux…f†0i,j @ Vkn


 @Di,j ÿ fi Z…ÿPi,j † ÿ fj Z…Pi,j † k @ Vn   @fj @fi @ Pi,j @ Pi,j Z…ÿPi,j † ÿ k Z…Pi,j † ‡ fi DZ…ÿPi,j † ÿ fj DZ…Pi,j † …47† ‡ Di,j @ Vkn @Vn @ Vkn @ Vkn

Recall that Vkn is any one of the variables in the physical model, and in this case the nodes are in the set n $ [i, j]. The implementation of Newton's method required that the function ( fi) evaluate to zero. In the one dimensional case without source terms the Newton function ( fi) is fi ˆFlux…fi †e ÿ Flux…fi †w ˆ 0 Flux…fi †e ˆ Flux…fi †oP,E Flux…fi †w ˆ Flux…fi †oW,P


Applying the analysis of Eqn. (47) to ®nd the total Jacobian produces the following partial di€erential of fi. The partial for the node in the east (E) direction of fi is @fi @Die i ˆ …fP Z…ÿPie † ÿ fiE Z…Pie †† k @VE @VkE ‡ Die

! i @fiE @ P i i ÿ k Z…Pie † ‡ ke …fP DZ…ÿPie † ÿ fE DZ…Pie †† @ VE @VE


The values Pie , Die , Piw , and Diw are the values of the Peclet number and di€usivity for equation i at the node faces (the nodes are W, P, E). The value for the Jacobian at the


Hydrodynamic simulation of semiconductor devices

central node location (P) is given by  @Diw ÿ i  @fi @ Die ÿ i i i i ˆ Z…ÿP † ÿ f Z…P † ÿ f fW Z…ÿPiw † ÿ fiP Z…Piw † e e P E k k k @VP @ VP @ VP ! i  @Pie ÿ i i i @ fP i i i Z…ÿP e † ‡ k fP DZ…ÿPe † ÿ fE DZ…Pe † ‡ De @VkP @VP !  @ fiP @ Piw ÿ i i i i i i ÿ Dw ÿ k Z…Pw † ‡ k fW DZ…ÿPw † ÿ fP DZ…Pw † @VP @VP


The Jacobian element in the west direction (W) is given by  @ fi @ Diw ÿ i ˆ fW Z…ÿPiw † ÿ fiP Z…Piw † k k @VW @VW ‡


 @fiW @Piw ÿ i i Z…P † ‡ fW DZ…ÿPiw † ÿ fiP DZ…Piw † w k k @ VW @ VW



Recall that i $ [1, n] and k $ [1, n], therefore, Eqns (49)±(51) must be evaluated n2 times. In the two-dimensional model this analysis of the partial di€erentials is repeated in the other spatial direction and will also contribute two additional terms to the partial at the central node due to nodes above (N) and below (S) the central node. 3.6.2. General partial di€erential of the two-dimensional source term augmented discretization equation. Now that the Jacobian elements for the one-dimensional discretization equation with no source terms have been de®ned, this will form the basis for the derivation of the general partial di€erential equation for the two-dimensional source term augmented discretization equation (Eqn. (19)). The ®rst step is to re-write the source term augmented discretization formula in two dimensions, see Fig. 3 for the node cluster,      ÿÿÿÿ 4 de dw o o ~ ÿ Sw Lw Cÿ …Pw † ‡ f ˆ r  Flux ÿ S1DyY Fluxe ÿ Fluxw ‡ Se Le C‡ …Pe † ‡ Le Lw      dn ds ÿ Ss Ls Cÿ …Ps † ‡ ‡ DxX Fluxon ÿ Fluxos ‡ Sn Ln C‡ …Pn † ‡ Ln Ls ÿ SP DxDy…X ‡ Y ÿ XY† ˆ 0


Where Flux8 are the discretized ¯ux equations (de®ned in Eqn. (18)) in the four directions for the zero source term discretization. The reader can expand these terms, but as will be seen, the partial of these terms will produce results which contain terms similar to Eqn. (47). As in the previous section, de®ne Vkm as one of the variables to be evaluated. The Jacobian for the variable Vkm at node m is obtained by (m $ [E, W, N, S, and P]) taking the partial of the Newton function (fi) with respect to Vkm ,


A. W. Smith and K. F. Brennan

  @fi @ Fluxie @Fluxiw @ Sie de i ˆ Dy Y ÿ ‡ k Le C‡ …Pe † ‡ Le @ Vkm @Vkm @ Vkm @ Vm ‡

Sie Le

!   @C‡ …Pie † @Siw dw @Cÿ …Piw † i i ÿ k Lw Cÿ …Pw † ÿ ÿ Sw L w Lw @ Vm Vkm @Vkm

    ! @Y de dw i i i i i i Fluxe ÿ Fluxw ‡ Se Le C‡ …Pe † ‡ ÿ Sw Lw Cÿ …Pw † ÿ ‡ @V km Le Lw   @ Fluxin @ Fluxis @Sin dn i ÿ ‡ k Ln C‡ …Pn † ‡ ‡ Dx X Ln @Vkm @Vkm @ Vm

!!   i i i @ C …P † @S d @C …P † ‡ s ÿ n s ÿ ks Ls Cÿ …Pis † ÿ ÿ Sis Ls ‡ Sin Ln Ls @ V km @Vm @Vkm

    ! @X dn ds i i i i i i ÿ Ss Ls Cÿ …Ps † ÿ ‡ k Fluxn ÿ Fluxs ‡ Sn Ln C‡ …Pn † ‡ Ln Ls @Vm ÿ

   @SiP @X @Y @X @Y i DxDy…X ‡ Y ÿ XY† ÿ S DxDy ‡ ÿ Y ‡ X P @VK @Vkm @Vkm @Vkm @ Vkm m


where i $ [1, n] and k $ [1, n] as in the one-dimensional case where the source term was zero. Only three new partials have been introduced in this formulation, compared to the zero source term case, @ C‡ @Cÿ and @ Vkm @Vkm @ Si @SP and k @Vm @Vkm


The functional form of the partials for the parameter C produce the same equation regardless of whether C+ or Cÿ is being evaluated. The formula for these partials is given below ! @C @ Pim Z…ÿPim †Z…Pim † ÿ 1 ˆÿ k ÿ …55†  @ Vm Z…Pim † ÿ Z…ÿPim † 2 @Vkm Again, the only partial derivative that needs to be calculated is the partial of the Peclet number with respect to the unknown. The numerical implementation of the function in parentheses in Eqn. (55) is discussed further in a later section in this chapter. The partial derivatives of the parameters Y and X in Eqn. (53) are created by a sum of partials of the parameter Y (see the de®nition of Y and X in Eqn. (26)). The partial of the source term located at the central node (SP, the last partial in Eqn. (54)) is evaluated directly from the formula for the source term in the de®nition of the particular model equation. Notice that this term will be non-zero only for the partial with respect to the parameters at the central node. The other partials of the source terms required in the Jacobian involve the partial of the harmonic average of the source terms calculated at the control volume face (Si, the second partial in

Hydrodynamic simulation of semiconductor devices


Eqn. (54)). These partials are given by

  @Si @SI S2i di ˆ 1 ÿ Li @ Vki @ Vki S2I @ Si @ SP S2i di ˆ @ VkP @VkP S2P Li


for the surrounding elements and the central node, respectively. The de®nition of the harmonic average has been utilized for these partials. Notice that the partial of the harmonic average of a source term contains partials of the source terms at the surrounding nodes. Eqn. (56) is valid only for non-zero source terms. As de®ned in Pantankar(1) the harmonic source term was also only de®ned for non-zero source terms. From the preceding discussions the power of the general method presented will now be highlighted further. First, only the general source term ¯ux equation (Eqn. (52)) and the general partial derivative of this equation (Eqn. (53)) need to be programmed. The routines for the harmonic averages of certain parameters are included in the general code. Once programmed these equations become invariant, even when the problem is changed. Some of the parameters which are problem dependent are the di€usivity, Peclet number, and source terms. Additionally, the derivatives of these parameters with respect to the unknowns of the system must be supplied. Therefore, to add additional equations, such as those for satellite valleys, requires no modi®cation of the main code, only additions to the parameter equations. To change the program from the hydrodynamic model for semiconductor devices to some other system of equations only requires the creation of routines for D, P, S, d, (@D/@V), (@P/@V), (@S/@V), and (@f/@V) for each of the unknown variables in the new problem.

3.7. Numerical issues In all of the ¯ux equations Bernoulli functions arise in the discretization equations. The form of the Bernoulli function is Z…x† ˆ

x exp…x† ÿ 1


This function must be implemented with care on digital computers due to the exponential functions and the possibility of division by zero. In all of the preceeding sections of this chapter we have followed the approach given in Pantankar(1). At this point a deviation is required which is special to semiconductor device simulation. In numerical heat and mass transfer problems, the numerical values of the numbers do not vary over a large range and therefore the hybrid scheme of Pantankar(1) is sucient. In semiconductor device simulation the carrier concentrations can vary by more than 20 orders of magnitude making it necessary to have very accurate values for the Bernoulli functions. If the hybrid scheme is used, the numerical error for the minority carrier concentration becomes excessive. Alternatively, we have employed a di€erent approach in which the Bernoulli function is approximated in several di€erent ranges by di€erent functions. These approximations have been utilized within the calculations presented in this chapter and are

318 given as,

A. W. Smith and K. F. Brennan

8 ÿx x R ÿ 36:0 > > > > > > x > > ÿ36:0 < x R ÿ 8:0  10ÿ7 > > > exp…x† ÿ 1 > > > > x > > ÿ8:0  10ÿ7 < x R 2:0  10ÿ7 < 1:0 ÿ 2 Z…x† ˆ > > x exp…ÿx† > 2:0  10ÿ7 < x R 36:0 > > > 1 ÿ exp…ÿx† > > > > > x exp…x† 36:0 < x R 746:0 > > > > > > : 0:0 746:0 < x


The usage of the above formulations ensures that no singularities are obtained during the evaluation of Z(x). The limits on the values of the argument (x) given above were calculated for a personal computer using the algorithms suggested in Selberherr(3). In addition, the ¯ux equations and the Jacobian matrix also contain the Bernoulli function with the negative value of the Peclet number as the argument. Applying the above analysis would be computationally expensive, therefore, the following formula was used Z…ÿx† ˆ x ‡ Z…x†


Eqn. (59) is an identity of the Bernoulli function. Using this formulation within a computer code requires only one extra addition statement, instead of ®ve conditional IF statements and the exponential evaluations contained in Eqn. (58). The Jacobian matrix contains functions which were created when the derivative of the Bernoulli function was taken with respect to the unknown variable which may also pose mathematical diculty in their evaluation. These functions were denoted by DZ(P) and the formula for this function is (recall that the partial of the Peclet number (@P)/(@V) was included in the general Jacobian equation) DZ…P † ˆ Z…P †

1 ÿ Z…ÿP † P


From the equation it is clear that as the Peclet number approaches zero, the computer representation will experience a singularity. To avoid this the following formulae were utilized 8 1 P > > jPjR1:0  10ÿ3 ÿ ‡ < 2 6:0 …61† DZ…P† ˆ > > : Z…P † 1 ÿ Z…ÿP † elsewhere P As in the case with the Bernoulli function, the value of DZ for the negative argument of the Peclet number is also required. As an alternative to using the above evaluation for the negative argument the following formula is used DZ…ÿP † ˆ 1 ‡ DZ…P †


This function is again a result of using an identity of the Bernoulli function. The use of this formula is suggested so as to save computation time. When the source term augmented ¯ux equations and general Jacobian are used, two additional functions were introduced which could lead to numerical instabilities. These two

Hydrodynamic simulation of semiconductor devices


functions are given by C‡ …P † ˆ

Z…P † ÿ 1 Z…ÿP † ÿ 1 and Cÿ …ÿP † ˆ P P


A quick observation shows that these are similar to Eqn. (60) divided by a Bernoulli function. However, the use of these equations could be computationally intensive and can lead to singularities as the Peclet number becomes large. A more ecient method would be to use 8 Z…P † ÿ 1 > > jPjr36:0 < P …64† C‡ …P † ˆ > DZ…ÿP † > :ÿ elsewhere Z…ÿP † Using the identity of the Bernoulli function the negative argument parameter is Cÿ …ÿP † ˆ 1 ‡ C‡ …P †


Finally, the derivative of the functions listed above (C+ and Cÿ) are required for the generalized Jacobian matrix. Fortunately, the formula for the partial derivative is the same for both C+ and Cÿ, DC…P † ˆ

Z…P †Z…ÿP † ÿ 1 P2


The formulation which avoids the singularity for the computer implementation of this equation is 8 ÿ…1 ÿ 0:052P 2 † > > jPjR1:0  10ÿ3 < 12 …67† DC…P † ˆ > Z…P †Z…ÿP † ÿ 1 > : elsewhere …Z…P † ÿ Z…ÿP †† Although the use of the formulations presented above require longer code, the e€ort is required to ensure the computational integrity of the program. 3.8. Node position and general boundary conditions The placement of nodes is an important consideration for the construction of the simulation domain. In general, the nodes are not equally spaced throughout the simulation domain. This is because, in some regions, node placement must be very ®ne in order to capture sudden changes in the potential, carrier concentration, etc. However, in other regions of the device, the important simulation quantities do not change suddenly and hence do not require extensive grid re®nement. Fine grid re®nement in regions which do not require it wastes computational resources. Therefore, variable grid spacing is necessary in the vast majority of device simulations. One method is to de®ne the placement of the nodes ®rst, then de®ne the control volume faces midway between the nodes. In this case, the node is not necessarily in the center of the control volume. If the node spacing is non-uniform, the node is clearly not in the center of the volume. This method is a little more accurate in terms of the interpolated pro®le. However, the assumption of a source term uniform over the whole control volume becomes questionable since the node is not at the center. A second method of grid generation is to place the node at the center of the control volume. In this case the node volume faces do not lie half way between the nodes. This method o€ers the advantage that the control volume faces can be placed at material interfaces, then the node locations follow later. A


A. W. Smith and K. F. Brennan

®nal advantage of the latter method of mesh generation is related to the de®nition of the boundaries. The ®rst method leads to half sized control volumes which require special discretization equations. In the latter method the volumes completely ®ll the simulation volume. A control volume of zero width is then de®ned for the speci®cation of the boundary condition. As previously stated, the solution of the general di€erential equation requires the speci®cation of the boundary conditions. As is always the case for partial di€erential equations, a unique and well-behaved solution requires the speci®cation of the boundary conditions along all of the surfaces of the device. These boundary conditions may be Dirichlet (®xed values) or Neumann in nature (¯ux or derivative) or a mixture of both. The boundary conditions must also be discretized and the partials created for the Jacobian matrix. The matix of equations created must be iteratively solved to provide the solution. Iteration is required due to the non-linear nature of the equations. This completes our discussion of the solution of the general di€erential ¯ux equation. Using the techniques outlined above, a speci®c hydrodynamic model can now be formulated. We begin with a brief history of hydrodynamic simulation. Next we derive the details of the hydrodynamic simulation and present some calculated results. 4. H Y D R O D Y N A M I C M O D E L

4.1. Background and history A brief history and introduction to the hydrodynamic model will now be presented. Often major distinctions have been highlighted between hydrodynamic models and energy balance models. However, these models are really very similar and only di€er in terms of where within the derivation of the balance equations assumptions are made on the energy ¯ux and distribution function. Due to this similarity the hydrodynamic and energy balance models will be generically called hydrodynamic models. The ®rst reference to a hydrodynamic model is found in the work of Stratton(20) in which the basic energy balance equations were derived from the Boltzmann transport equation. Assuming a speci®c form for the relaxation parameter and expanding the distribution function in spherical harmonics, Stratton produced a set of coupled di€erential equations for the carrier momentum and energy. The equations were applied to simulate the performance of a Schottky diode. In 1970, Blotekjaer(21) introduced an energy balance model in which no assumption for the form of the relaxation parameter was required. To produce a simple set of equations the distribution function was again expanded in spherical harmonics. The di€erential equations were created by using the average value theorem of di€erential calculus(22). Although the equations of Stratton and Blotekjaer are valid under the assumptions used in their derivation, Stratton's equations hold more physical meaning since the average value theorem was not employed(22). The set of di€erential equations for the energy balance model were neglected in device simulation for some time after the work of these two authors. During the decade of the 70's attempts were made to simulate devices using extensions of the drift±di€usion equations, without resorting to the full energy balance models. These models used a Joule heating term (JE) to describe the energy transfer from the carriers to the lattice. In the work by Gaur and Navon(23), the standard drift±di€usion equations were augmented by Fourier's heat ¯ow equation for the lattice energy. The Joule heat term was added to Fourier's equation as a source term(23). In 1977, this work was extended by(24) adding the e€ect of the thermal di€usion of the carriers to the model while retaining the Joule heating term. These models predicted device heating in regions of large potential gradients and high current ¯ow. In 1978, the e€ect of the energy transfer due to

Hydrodynamic simulation of semiconductor devices


recombination events was added to the Joule heating term by Adler(25). Therefore, local device heating was predicted in regions of the device where the recombination was high. The use of the Joule heating term assumes that the carriers are in thermal equilibrium with the lattice at all times. In modern devices this assumption is no longer valid, therefore, the full set of energy balance equations must be solved. In 1981, Cook and Frey(26) used a simple energy balance analysis to describe ballistic transport in GaAs devices. This early work showed the importance of including the di€usion e€ect for the charge carriers even in regions of the device where carrier ¯ux is dominated by energy e€ects. A more re®ned model was presented by these authors in 1982 for Si and GaAs MESFET's(27). These models required detailed knowledge of the device structure to provide simpli®cations to the numerical methods used for the simulation. Cook later presented a more general model for Si BJT's which provided a new phenomenological energy relaxation time model to describe the energy interaction between the electrons, holes, and lattice(28). Mains, Haddad, and Blakey used an energy balance model to simulate GaAs IMPATT diodes in 1983(12). The author's simulations su€ered from numerical instability due to the non-linear nature of the di€erential equations. The phenomenological energy relaxation model of Cook was also replaced by relaxation parameters determined from steady state Monte Carlo calculations. A numerical technique to stabilize the non-linear nature of the energy balance equations was proposed by Tang in 1984(14). The approach is similar to that used in the Scharfetter± Gummel technique for the isothermal current equations(13). In addition, the Scharfetter± Gummel method was extended to account for non-isothermal conditions in the current equations. A similar derivation was proposed by McAndrew et al.(29) in 1985 and, yet, again in 1987 by Szeto and Reif(30). These extensions, in a form similar to the standard Scharfetter±Gummel method, allowed the incorporation of the energy balance equation into simulation models with a smaller amount of e€ort than previously possible. The complete set of six partial di€erential equations constituting the hydrodynamic model was proposed in 1985 by Wang(31). Prior to this time, energy balance models assumed that the lattice was an in®nite sink for carrier energy (Tl constant) or as in the case of the Joule heating models the carriers were assumed to be in equilibrium with the lattice. Wang reintroduced the lattice temperature equation and coupled this to the carrier energy equations through the source terms. Therefore, a loss in energy of the electron system was re¯ected by a gain in energy of either the hole system or the lattice. This was not a local e€ect (such as the Joule term) but was related to the energy relaxation time of the charge carriers. Additionally, Wang introduced an impact ionization formulation which was dependent upon the carrier energy, not the current as in most drift±di€usion models. In 1986 a multi-dimensional discretization scheme for the hydrodynamic model was introduced by Rudan and Odeh(32). The device domain was discretized into triangular elements and discretization equations in the form of the Scharfetter±Gummel formulation were presented. A very similar model, in terms of the discretization equation and set of di€erential equations, was introduced in 1988 by Forghieri(33). A separate two-dimensional model for GaAs MESFET's was introduced by Snowden and Loret(34) in the previous year. Also in 1987, McAndrew et al.(35) and Azo€ (36) re-derived the energy balance equations based on the perturbation approach to the distribution function. By making assumptions and simpli®cations both models were shown to produce forms similar to those arrived at by Stratton and Blotekjaer. In a second paper Azo€ (37) proposed a closed form method for solving the energy and momentum conservation equations. Azo€ also extended the energy balance equations to non-parabolic heterostructure and degenerate semiconductors(38). The interest in hydrodynamic modelling and its application grew even more rapidly after 1987, as characterized by the large amount of literature published on this method. For


A. W. Smith and K. F. Brennan

instance, Ou and Tang(39) modeled sub-micron Si BJT's in 1987 using the approach previously proposed by Tang(14). Three-dimensional simulations of VLSI devices using an energy balance model were undertaken by McAndrew et al. in 1988(40). The importance of the thermal di€usion e€ect of the carriers was veri®ed in that work. Sub-micron GaAs MESFET's were again simulated in 1988 by Feng and Hintz(41) using an extension similar to that proposed by Tang to stabilize the energy and current ¯ow equations. New numerical treatments are being explored such as shock capturing algorithms to stabilize the non-linearities(42). For more recent applications of the hydrodynamic model and topics related to the hydrodynamic model the reader should refer to references(2, 38, 42±47). Recently, non-parabolic formulations of the hydrodynamic model have been formulated. Due to the application of the model to III±V narrow band gap materials the inclusion of the non-parabolicity of the bandstructure is necessary to accurately model the device physics(48). One of the non-parabolic formulations is presented in detail within this chapter. The reader can reduce the form to parabolic bands by setting a to zero. 4.2. Derivation of the transport model In Section 3, the Boltzmann equation was derived and the basics of its solution were presented. In this section, we illustrate the steps taken in the solution of the Boltzmann equation using the method of moments. When possible, analogies to the drift±di€usion model and other hydrodynamic models in the literature will be drawn. The ®nal result of applying the method of moments to the Boltzmann transport equation produces a system of non-linear partial di€erential equations which describe the number, momentum, and energy (temperature) of the electron and hole sub-systems, and provide the energy (temperature) of the lattice. A general moment equation can be formulated from the Boltzmann equation after applying several vector and tensor identities(49, 50). We start by multiplying the Boltzmann equation by a general function, Y, as, @ ~ r ~ r  … f Y~ug † ÿ f …~ug  r ~ k †Y ˆ Y ÿ… f ÿ fo † ~ r †Y ÿ f …F …Yf † ‡ r @t h t


The quantity Y is the parameter for which the expected value is to be computed. Following the example of Woolard et al.(51) the following values of Y are used: 1 for the particle continuity equation, ug, the carrier drift velocity, for the ¯ux equation, W, the carrier energy, for the energy equation, and Wug for the energy ¯ux equation. Setting Y = 1 and Y = W in Eqn. (68) gives the following conservation equations for the carrier concentration and carrier energy @f ~ r  … f ~ug † ˆ ÿ… f ÿ fo † ‡r @t t


@ ~ r  … fW~ug † ‡ f …r ~ r ec  ~ug † ˆ W ÿ… f ÿ fo † … fW † ‡ r @t t



In the derivation of Eqns (69) and (70) we have made use of the fact that r and k, the position and wavevector, are assumed to be independent, therefore the order of the di€erentials can be switched. Also the de®nition of the group velocity has been utilized 1~ ~ug ˆ r kW h


Hydrodynamic simulation of semiconductor devices


Eqns (68)±(70) are valid for all choices of the dispersion relations governing the energy vs. wavevector, E(k), relationship. In this chapter, we will restrict discussion to only parabolic and non-parabolic formulations of the E(k) relation. 4.2.1. Particle ¯ux equations. The ¯ux equation for the hydrodynamic system can be developed from the general conservation equation substituting in for the moment (Y = ug). The following general conservation equation is found,   2 ~~ @ 1 ~ I ‡ ~ ~ r  … f v~g v~g † ÿ f h …kk†  r …~ vg f † ‡ r ÿf F @t m…k† m…W † m…k† ~~ ~  3…krk m…k†† ˆ ÿ… f ÿ fo † v~g ‡f F m2 …k† t


Notice that the third, fourth and ®fth terms on the left hand side are tensor products. The identity matrix (I) is the tensor in the fourth term. The factor of 3 in the ®fth term is due to the order of parenthesis in the original moment equation. This moment equation can not be processed further until some functional form of the e€ective mass is assumed, which depends on the choice of the dispersion relation. In this chapter, we consider the Kane formulation of W(k) relations given by the following equation W ‡ aW 2 ˆ

…hk†  2 2m


The corresponding k or W dependent e€ective masses are given by v ! u 2 2 u h k m…W † ˆ me …1 ‡ 2aW † m…k† ˆ me t1 ‡ 4a 2me


where we have used the momentum as de®ned by hk/m.  Eqn. (74) applies to the Kane or a formulation of the dispersion relation. The mass formulation was derived earlier(51) and is a result of the momentum±velocity relation given above. For this reason it is not the usual de®nition of the e€ective mass. Substitution of Eqn. (74) into Eqn. (72) will lead to a nonparabolic particle ¯ux equation. The further assumption of equi-partition of energy is also required. The resulting ¯ux equation is (the positional subscript on all of the gradient operators has been removed)      ~  @ ÿ 2 ~ W…1 ‡ aW † 4 W…1 ‡ aW † rm e ~ug f ‡ r f ˆ f 1 ‡ aW…1 ‡ aW † @t 3me 3 …1 ‡ 2aW †2 …1 ‡ 2aW †4 m2e


  4   W 2 ÿ 1 ‡ …1 ‡ aW †2 ! f 3 ~ ~ e ˆ ÿ… f ÿ fo † ~ug ra ‡ re f 4 3 t me …1 ‡ 2aW † me …1 ‡ 2aW †



A. W. Smith and K. F. Brennan

This equation must be integrated over all k space or equivalently over energy using the density of states(52). Before the integration is performed one more assumption must be made, that the relaxation time is independent of k or W. Typically, both sides of Eqn. (75) are multiplied by the relaxation time before the k space integration is performed. If the constant relaxation time assumption is not made then the energy dependence of the relaxation time must be moved through the gradient operator on the second term of Eqn. (75) and a term accounting for the gradient of the relaxation time must be re-created. For instance, substituting in an energy dependent relaxation time, t(W ), into the second term above obtains,     2 ~ W…1 ‡ aW † 2 ~ W…1 ‡ aW † r f r f t…W † ˆ t…W † 3me 3me …1 ‡ 2aW †2 …1 ‡ 2aW †2 …76† 2 W…1 ‡ aW † ~ rt…W † f ÿ 3me …1 ‡ 2aW †2 Notice that an additional term involving the gradient of the relaxation time now appears. Since in most cases the energy dependence of the relaxation time is not a constant power of the energy this adds more complexity to the derivation and is commonly ignored. However, Monte Carlo results indicate that the relaxation time is nearly a constant over a certain energy range so there is some justi®cation in neglecting the term involving the gradient of the relaxation time. Changing the integration from k space to energy space, substituting for the mobility, (m = t/m), and recognizing that ug fo is equal to zero, gives the following integral equation   4 3=2 # " … 1 ‡ aWY ~ e … …WY† rm @ …n~vg † 2m ~ …WY†3=2 3 dW ÿ m ‡ r g f g f dW t 3 C me @t 3 C   4 … … W 2 ÿ 1 ‡ Y 2 …WY†1=2 …WY†1=2 3 ~ ~ dW ‡ m re g f dW ˆ n~ vg ‡ mrag f c C3 C2


Where the short hand notations gˆ

1 2p2

2me h2


and Y ˆ …1 ‡ aW †; C ˆ …1 ‡ 2aW †


have been used. To produce closed form solutions for the integral in Eqn. (76) the binomial expansion must be used repeatedly and all terms of order a2 or higher are set to zero. The ¯ux equation for the a formulation then becomes " … #     ~ e … rm @ …n~vg † 2m ~ aW 19aW 3=2 ‡ r g f W g f W 3=2 1 ÿ 1ÿ dW ÿ m dW t me @t 3 2 6 …

~ ‡ mrag f W 5=2

   … 1 5aW 7aW 1=2 ~ ‡ dW ‡ mrec g f W dW ˆ n~ vg 1ÿ 3 6 2


Hydrodynamic simulation of semiconductor devices


4.2.2. Energy ¯ux equations. To obtain the energy ¯ux equation, the substitution, Y = Wug, is made within the general conservation equation given by Eqn. (68). This becomes,   ~ @ h2 …k~k† 1 ~  ~ug ~ug ~ ~ ~ …W~ug f † ‡ r  … fW~ug ~ug † ÿ f ~ug ~ug  rW ÿ fW r ÿ fF @t m…k† m…W † ~ ~  3W…krk m…k†† ~  W I ‡ fF ÿ fF m…k† m2 …k† ÿ… f ÿ fo † W~ug ˆ t


Notice that the fourth through seventh terms on the left hand side are tensor products, I again being the identity matrix. As in the case of the particle ¯ux moment the energy ¯ux moment equation can not be processed further until some functional form of the e€ective mass is assumed. Using the same assumptions as for Eqn. (79) an energy ¯ux equation can be formulated. At this point in the derivation a recursion relation must be formulated, a distribution function assumed, or some other mathematical method {Minimum±Maximum theorem} must be used to provide mathematical closure for Eqn. (79) and the energy ¯ux equation. Since, a goal of this paper is the formulation of models suitable for the numerical simulation of devices the ®rst two options are explored(53). A recursion relation would allow moments of higher order to be approximated by lower order moments, the lower orders are calculated from the conservation equations. This option does not require that a speci®c form of the distribution function be used unless the recursion relations are based on a speci®c distribution. However, in the case of nonparabolic bands the standard recursion relations may no longer be applicable(51). Therefore, this option was not pursued. The other option, and the one chosen for this work, is to assume a speci®c form for the distribution function; higher moments can then be calculated based on the known distribution function. Some of the choices for the distribution function include heated Maxwellian, shifted and heated Maxwellian, heated Fermi±Dirac, or shifted and heated Fermi±Dirac. Since the Maxwellian distributions can be recovered by relaxing the degeneracy, the Fermi±Dirac distributions were the only ones considered for this work. In the parabolic band formulation the energy of the shifted Fermi±Dirac distribution is 3 F3=2 …Z† me 2 ‡ v ˆ KTe 2 d 2 F1=2 …Z†


In the non-parabolic formulation(48) the above relationship for the energy will no longer be true due to the change in the density of states. Also higher order powers of the energy are required to close the relationships in the above formulation. This will require cross product terms involving the temperature and the velocity. Due to these limitations and the fact that all the formulations break down as the energy rises, the heated Fermi±Dirac distribution was used to close the relationships. The ¯ux equation in the binomial a (non-parabolicity)


A. W. Smith and K. F. Brennan

formulation is


3 10 6 7 F21 …Zc † ‡ aKTe F3 …Zc †F1 …Zc † 6 7 4 2 2 2 ~ 6 ÿn~ vg ˆ mKTe 6 !7 7rn F 1 …Zc †F …Zc † 3 ÿ2 4 5 15 2 F1 …Zc †Fÿ1=2 …Zc † ‡ aKTe F21 …Zc † ‡ 4 2 2 2 2

3 21 aKTe F3 …Zc † mnKTe 6 7~ 4 2 ! ‡ mn4 2 5rec ÿ 15 15 F1 …Zc † ‡ aKTe F3 …Zc † Fÿ12 …Zc † ‡ aKTe F12 …Zc † 4 2 2 4 F1 …Zc † ÿ



6 F21 …Zc † ‡ 20aKTe F5 …Zc †Fÿ1=2 …Zc †7 ~ 6 7 rm 2 6 ! 7  6F12 …Zc † ‡ 2 7 m 15 4 5 2 F1 …Zc † ‡ aKTe F3 …Zc † 4 2 2


mnKTe !

15 F1 …Zc † ‡ aKTe F3 …Zc † 4 2 2

15 Fÿ12 …Zc † ‡ aKTe F1 …Zc † 4 2




6 65 2 6 62 Fÿ12 …Zc †F32 …Zc † ÿ F12 …Zc † ÿ 4

F31 …Zc † 2

2 F1 …Zc † ‡ 2

15 aKTe F3 …Zc † 4 2

7 ~ e 7 rT !7 7 Te 5

2 6 10F1 …Zc †F3 …Zc † 6 175 2 2 !6 !ÿ aKTe F7 …Zc † ÿ 6 4 2 15 15 4 F1 …Zc † ‡ aKTe F3 …Zc † Fÿ12 …Zc † ‡ aKTe F1 …Zc † 4 4 2 2 2 4mn…KTe †2

3 F21 …Zc †F3 …Zc † 2 ! 2

‡ F1 …Zc † ‡ 2

15 aKTe F3 …Zc † 4 2

Fÿ12 …Zc † ‡

15 aKTe F1 …Zc † 4 2

7 7 ~ !7 7ra 5


By making similar assumptions, substitutions, and approximations {binomial expansion, Fermi±Dirac statistics, equipartition of energy...} the energy ¯ux equations for the Kane or non-parabolicity (a) dispersion relationships become

! 3 15 7 6F32 …Zc † F12 …Zc † ‡ 4 aKTe F32 …Zc † ÿ 4 aKTe F12 …Zc †F52 …Zc †7 7 5m…KTe †2 6 ~ ~ 6 ! !7 ÿS ˆ 6 7rn 2 15 15 4 5 F1 …Zc † ‡ aKTe F3 …Zc † Fÿ1 …Zc † ‡ aKTe F1 …Zc † 4 4 2 2 2 2 2

2 6 5m…KTe † 6 6 ‡ 6 2 4

3 23 2 F3 …Zc † ÿ aKTe F5 …Zc † 7 7 4 2 2 ~ c ÿ 5mn…KTe † !7 rE 7 2 15 5 F1 …Zc † ‡ aKTe F3 …Zc † 4 2 2

2 6 6 6 6 4

3F3 …Zc †



15 2 Fÿ12 …Zc † ‡ aKTe F1 …Zc † 4 2

14aKTe F7 …Zc † 2

15 F1 …Zc † ‡ aKTe F3 …Zc † 4 2 2


3 21aKTe F5 …Zc †F1 …Zc † 2 !2

ÿ 8 F1 …Zc † ‡ 2


15 aKTe F3 …Zc † 4 2

5mn…KTe †2 15 2 F1 …Zc † ‡ aKTe F3 …Zc † 4 2 2

Fÿ12 …Zc † ‡

15 aKTe F1 …Zc † 4 2

7 ~ 7 rm !7 7 m 5


2 3 6 F1 …Zc †F5 …Zc † ÿ aKTe F1 …Zc †F3 …Zc † 65 21 20 2 2 2 6 !2  6 F5 …Zc † ÿ aKTe F7 …Zc † ‡ 8 2 15 42 2 F1 …Zc † ‡ aKTe F3 …Zc † 4 2 2 ! 3 25 101 aKTe F1 …Zc †F5 …Zc †7 3F3 …Zc † F1 …Zc † ‡ aKTe F3 …Zc † ÿ 4 8 2 2 2 2 2 ~ e 7 rT 7 ! ÿ 7 Te 15 5 2 Fÿ12 …Zc † ‡ aKTe F1 …Zc † 4 2 2 5mn…KTe †3


2 F1 …Zc † ‡ 2

15 aKTe F3 …Zc † 4 2

6 6 !6 6 4

3F23 …Zc † 2

2 Fÿ12 …Zc † ‡


15 aKTe F1 …Zc † 4 2

3 7 105aKTe F1 …Zc †F3 …Zc †F5 …Zc † 7 2 2 2 ~ 7ra ! ! ‡ 7 15 15 5 16 F1 …Zc † ‡ aKTe F3 …Zc † Fÿ12 …Zc † ‡ aKTe F1 …Zc † 4 4 2 2 2

315 aKTe F9 …Zc † 16 2 …83†


A. W. Smith and K. F. Brennan

Again, the reader can verify that the non-parabolic formulations reduce to the parabolic case as the non-parabolicity factors are decreased. As in the case of the particle ¯ux equations the energy ¯ux equations can be discretized using normal techniques(54).

4.3. System of equations for the hydrodynamic model The ®nal form of the system of six balance equations for the hydrodynamic model can now be given as: Electron continuity: 

@n @t

~  …n~vn † ˆ ÿR ‡ G ‡r


~  … p~ ‡r vp † ˆ ÿR ‡ G


Hole continuity: 

@n @t

Electron energy balance: @…nWn † ~ ~ c ÿ n~ ~ ˆ ÿ…nWn ÿ nWo † ‡ r  …n~ qn † ‡ n~ vn  re vn  rw @t t3n


Hole energy balance: @ … pWp † ~ ~ v ÿ w ÿ Eg † ˆ ÿ… pWp ÿ pWo † ‡ r  … p~ qp † ‡ p~ vp  r…e @t t3p


Poisson's equation: ~ ˆ q … p ÿ n ‡ N ‡ ÿ N ÿ† ~  …xr rc† r a d xo


Lattice heat ¯ow:  rcp

@ Tl @t

~  …krT ~ 1 † ˆ …nWn ÿ nWo † ‡ … pWp ÿ pWo † ‡r t3n t3n

The auxiliary equations which must accompany these are: Electron ¯ux equation: Given above as Eqn. (82)



Hydrodynamic simulation of semiconductor devices

Hole ¯ux equation: 2

3 10 6 7 F21 …Zv † ‡ aKTh F3 …Zv †F1 …Zv † 6 7 4 2 2 2 ~ 7rp ÿp~ vg ˆ mKTh 6 ! 6 Fÿ12 …Zv †F3 …Zv † 7 4 5 15 2 2 F1 …Zv †Fÿ12 …Zv † ‡ aKTh F1 …Zv † ‡ 4 2 2 2 2

3 21 aKTh F3 …Zv † mpKTh 6 7~ 4 2 ! ÿ mp4 2 5rev ÿ 15 15 F1 …Zv † ‡ aKTh F3 …Zv † Fÿ12 …Zv † ‡ aKTh F1 …Zv † 4 2 2 4 2 F1 …Zv † ÿ



6 F21 …Zv † ‡ 20aKTh F5 …Zv †Fÿ12 …Zv †7 ~ h 6 7 rm 2 2 ! 7 6 6F12 …Zv † ‡ 7 mh 15 4 5 2 F1 …Zv † ‡ aKTh F3 …Zv † 4 2 2


mpKTh !

15 F1 …Zv † ‡ aKTh F3 …Zv † 4 2 2

15 Fÿ12 …Zv † ‡ aKTh F1 …Zv † 4 2




6 65 2 6 62 Fÿ12 …Zv †F32 …Zv † ÿ F12 …Zv † ÿ 4

F31 …Zv † 2

2 F1 …Zv † ‡ 2

15 aKTh F3 …Zv † 4 2

7 ~ h 7 rT !7 7 Th 5

2 6 10F1 …Zv †F3 …Zv † 6 175 2 2 !6 !ÿ aKTh F7 …Zv † ÿ 6 4 2 15 15 4 F1 …Zv † ‡ aKTh F3 …Zv † Fÿ12 …Zv † ‡ aKTh F1 …Zv † 4 4 2 2 2 4mp…KTh †2

3 F21 …Zv †F3 …Zv † 2 ! 2

‡ F1 …Zv † ‡ 2

15 aKTh F3 …Zv † 4 2

Fÿ12 …Zv † ‡

15 aKTh F1 …Zv † 4 2

Electron energy ¯ux equation: Given on p. 327 as Eqn. (83).

7 7 ~ !7 7ra 5



A. W. Smith and K. F. Brennan

Hole energy ¯ux equation: 2

! 3 15 7 6F32 …Zv † F12 …Zv † ‡ 4 aKTh F32 …Zv † ÿ 4 aKTh F12 …Zv †F52 …Zv †7 7 5m…KTh †2 6 ~ ~ 6 ! !7 ÿS ˆ 6 7rp 2 15 15 4 5 F1 …Zv † ‡ aKTh F3 …Zv † Fÿ12 …Zv † ÿ aKTh F1 …Zv † 4 4 2 2 2 2 3 6 F …Z † ÿ 23 aKTh F …Z † 7 v v 5 3 7 5mpKTh 6 4 2 2 ~ 6 !7 ÿ 6 7rev 2 15 4 5 F1 …Zv † ‡ aKTh F3 …Zv † 4 2 2 2 6 5mp…KTh †2 6 6 ÿ 6 2 4


8 F1 …Zv † ‡ 2

14aKTh F7 …Zv †



15 2 Fÿ12 …Zv † ‡ aKTh F12 …Zv † 4 21aKTh F5 …Zv †F1 …Zv † 2 !2



3F3 …Zv †

15 aKTh F3 …Zv † 4 2

5mp…KTh †2 15 2 F1 …Zv † ‡ aKTh F3 …Zv † 4 2 2 2

Fÿ12 …Zv † ‡

15 F1 …Zv † ‡ aKTh F3 …Zv † 4 2 2 3

15 aKTh F1 …Zv † 4 2


7 ~ h 7 rm !7 7 mh 5


3 6 F1 …Zv †F5 …Zv † ÿ aKTh F1 …Zv †F3 …Zv † 65 21 20 2 2 2 !2 6 62 F52 …Zv † ÿ 8 aKTh F72 …Zv † ‡ 15 4 F1 …Zv † ‡ aKTh F3 …Zv † 4 2 2 ! 3 25 101 aKTh F1 …Zv †F5 …Zv †7 3F3 …Zv † F1 …Zv † ‡ aKTh F3 …Zv † ÿ 4 8 2 2 2 2 2 ~ h 7 rT 7 ! ÿ 7 Th 15 5 2 Fÿ12 …Zv † ‡ aKTh F1 …Zv † 4 2 2 5mp…KTh †3


2 F1 …Zv † ‡ 2

15 aKTh F3 …Zv † 4 2

6 6 !6 6 4

3F23 …Zv † 2

2 Fÿ12 …Zv † ‡


15 aKTh F1 …Zv † 4 2 3

7 105aKTh F1 …Zv †F3 …Zv †F5 …Zv † 7 2 2 2 ~ ! !7 ‡ 7ra 15 15 5 16 F1 …Zv † ‡ aKTh F3 …Zv † Fÿ12 …Zv † ‡ aKTh F1 …Zv † 4 4 2 2 2

315 aKTh F9 …Zv † 16 2


Hydrodynamic simulation of semiconductor devices


By substituting the ¯ux equations into the conservation relations, equations for the six unknown parameters (C, n, p, Te, Tl, and Th) are obtained. Eqns (88) and (89) are Poisson's equation and the lattice heat ¯ow equation respectively. The parameters xr and xo in Poisson's equation are the relative permitivity and the permitivity of free space, respectively. The terms introduced in Eqn. (89), r and cp, are the material density and heat capacity. Notice that losses in carrier energy in Eqns (85) and (86) are generation terms in Eqn. (89). The parameters in these equations are not constants, but are actually functions of space, temperature, and time. Also, note that this derivation used the Fermi±Dirac distribution. The rationale for retaining the Fermi functions of di€erent order in all of the hydrodynamic equations is quite simple. First Fermi statistics are the proper form. Boltzmann statistics are an approximation. Second, as degeneracy is approached (the reduced Fermi energy Z equal to zero) the values of the Fermi integrals of various orders begin to separate from each other. Since the ratio of these functions appear in all of the auxiliary equations, setting these ratios to one is not acceptable. Therefore, any simulation program which fails to take these terms into account will not mimic the true physical system since the e€ect of certain parameters, which are multiplied by these ratios, are not accurately represented. Notice that once the auxiliary equations are substituted into Eqns (84)±(88) through Eqn. (89) conservation equations are formed. These equations may now be manipulated into the form of the general di€erential equation, Eqn. (11). Once this task is performed the remainder of the general fomulation may be applied such as the Peclet number, di€usion coecient, and source terms. In the computer simulation of the ¯ux equations, the code must accurately represent the asymptotic limits of the equations without introducing the possibility of division by zero. For instance, in both the carrier and carrier energy ¯ux equations the di€usivity and Peclet number contain terms of the form   m…L† …92† ln m…0† As m(L) approaches m(0) this term will evaluate to zero, but no numerical singularities will arise in the computer implementation because this term never appears in a denominator. However, the equations also contain terms with   Te …L† ln Te …0† …93† A rTe in the Peclet number and the inverse in di€usivity which may lead to numerical singularities. A is the drift component of the Peclet number. In the limit as Te(L) approaches Te(0) this term should reduce to (L/Te(0)), however, in the computer formulation both the numerator and denominator may evaluate to zero, producing division by zero. A suitable approximation to avoid computer division by zero as Te(L) approaches Te(0) is 8   Te …L†   > > ln > Te …L† > > Te …0† ln < jTe …L† ÿ Te …0†jR0:01 Te …0† rTe …94† ˆ > rTe > > L …Te …L† ÿ Te …0††L > > : ÿ elsewhere Te …0† Te …0† This approximation, normalized by (Te(0)/L), is shown in Fig. 4 for two values of


A. W. Smith and K. F. Brennan

FIG. 4. Normalized approximation used to avoid numerical singularities due to the temperature terms in the exponent of the discretized ¯ux equations.

temperature (300 and 600 K plotted as diamonds and triangles, respectively) at the origin. The original normalized function is also plotted as a solid line in Fig. 4 in both cases. It is clearly seen that the approximation closely matches the function over the range speci®ed. Notice that the range of plotted temperature di€erentials in Fig. 4 is very small, however, the function can deviate from unity for larger temperature di€erentials such as those seen in device simulations. A rather simplistic approach would be to replace Eqn. (93) by a term representing the average temperature, A

2L Te …0† ‡ Te …L†


However, the inclusion of this term in the exponential term may lead to unacceptable errors as the temperature gradient increases. The error generated by the use of Eqn. (95) in the exponential term of the ¯ux equation is depicted in Fig. 5 for (A = 500/L where A is given as (c(L) ÿ c(0))/(L*K)), a reasonable number which represents a potential variation of 050 meV between two nodes, as a function of temperature. The temperature gradient term was not explicitly taken into account in this error calculation because the electric ®eld component is approximately ®ve times larger than the temperature component. However, depending upon the sign of the temperature gradient the error may be slightly larger or smaller as the temperature e€ect is included in the parameter (A). 4.4. Boundary conditions The discretization and subsequent formulation of the algebraic equations for the di€erential equations has not been formally completed. As observed from Fig. 1B, most of the nodes are contained within the interior of the device. However, some nodes are positioned on the boundary of the device. These boundary nodes will not have nodes

Hydrodynamic simulation of semiconductor devices


FIG. 5. Absolute value of the error if the average nodal temperature is used in the exponential as compared to using the exact formulation.

surrounding them from which algebraic equations for the parameters can be created. At these nodes additional information is required to de®ne the values of the parameters. In some cases, the value of the parameter is known explicitly at the boundary in question. In other cases, the ¯ux of the unknown is speci®ed. Mathematically these types of boundary conditions are classi®ed as Dirichlet and Neumann, respectively. The device designer is not speci®cally interested in the mathematical terminology of boundary conditions. Therefore, the boundary conditions are de®ned by physical considerations, not mathematical con®guration. Boundary conditions are the method by which the device designer is able to connect the domain outside the simulation area to the parameters being solved for within the discretized domain. The choice of boundary conditions will also a€ect the simulation time of the device. If the parameter is speci®ed by a ®xed value at the boundary then the interior nodes will have no e€ect on the parameter. Therefore, the rate of convergence is increased. On the other hand, if the ¯ux of the parameter is speci®ed many iterations are required to ensure that self-consistency has been achieved. In addition to specifying the boundary condition type, the choice of node structure must also be selected at each boundary. There are two types of boundary nodes permissable; a half a control volume or a zero width control volume. A half control volume has only one distance speci®ed. In this type of boundary condition the program must include the regular source term for the volume as well as any terms from the boundary condition. The second type, a zero volume node, does not have any of the associated source terms from the partial di€erential equation because the volume is identically zero. The only terms which need to be included are the boundary parameters. In addition, the zero width control volume may be easily used as interior nodes within the simulation domain to simulate abrupt heterostructures.


A. W. Smith and K. F. Brennan

4.4.1. Thermal boundary conditions. The thermal boundary conditions for the carrier energy equations will be assumed to be the same as those for the lattice temperature equation, as suggested by McAndrew(12). This implies that at the contacts the carriers instantly relax to the lattice temperature. The three boundary conditions available for the lattice temperature are: 1. A given boundary temperature. 2. A given heat ¯ux through the boundary (a heat ¯ux of zero implies symmetry). 3. The boundary heat ¯ux is speci®ed using a heat transfer coecient and ambient external temperature. These boundary conditions should be capable of modelling most cases encountered. 4.5. Electrical boundary conditions Compared to the thermal boundary conditions, the electrical boundary conditions are quite complex. For the carrier concentrations the boundary conditions are treated di€erently if the carrier is a majority type as compared to a minority type carrier. Then if the current boundary conditions are given in terms of an external resistor, the carrier concentration boundary conditions are coupled to the boundary condition for the potential (Poisson's equation) as well as to the other charge carrier through surface recombination velocities. i. Ohmic Contacts. Usually one assumes that an ohmic contact is ideally conducting. Therefore, no voltage is dropped across the contact. The boundary condition for the potential then becomes cbc ÿ cbi ˆ ca


Where cbc is the applied voltage at the contact and cbi is the built in potential due to the doping concentration. For the carrier concentrations one choice of the boundary conditions is given below p C 2 ‡ 4n2i ‡ C nˆ 2 p 2 C ‡ 4n2i ÿ C pˆ 2 where C ˆ Nd ÿ Na


This corresponds to an in®nite surface recombination velocity at the contact surface. ii. Free Surface. For a free surface, or the surface between a semiconductor and an insulator the boundary condition for the electrostatic potential is esemi

@c ˆ Qint @Xi


The parameter Qint is the charge associated with the interface. If there is no interface charge then the boundary condition becomes @c ˆ0 @ Xi


Hydrodynamic simulation of semiconductor devices


The boundary conditions for the carrier concentrations will depend upon satisfying the continuity equations at the surface with respect to the surface recombination velocity, Jn ˆ ÿ q Rsurf Jp ˆ q Rsurf


Where the user determines the value of the surface recombination velocity (Rsurf). iii. Symmetry Boundaries. By the de®nition of symmetry the following boundary conditions for the electrostatic potential, electron concentration, and hole concentration must exist @c ˆ0 @ Xi @n ˆ0 @ Xi @p ˆ0 @ Xi


It is important for this type of boundary condition to be applied with some care. Placing a symmetry boundary condition too close to a region in which it is not true physically could lead to the simulation producing incorrect results. 4.6. Advanced physics Derivation of the conservation equations under the parabolic band approximation results in rather simple coecients on the forcing terms in the ¯ux equations. By relying on the parabolic band approximation higher order energy transport e€ects due to variation in the band structure are neglected. Interest in accounting for band structure e€ects in hydrodynamic device simulation has begun to grow because parabolic models can not adequately account for high energy e€ects in semiconductors with non-parabolic band structures. Non-parabolic band formulations have a history dating back to the 1950s(48, 55±57). However drift±di€usion models and more speci®cally hydrodynamic simulators with non-parabolic band formulations are a very recent topic of research. Therefore they are presented in this paper. The Kane formulation presented here is not the only non-parabolic formulation, and possibly not even the optimal choice under certain conditions(58, 59). However, it does provide a common ground to other models presented in the literature and is currently the most widely used non-parabolic formulation. The reader should recall that the parabolic formulation can be recovered by setting all of the non-parabolicity factors to zero. The basic hydrodynamic model described in the previous section is the minimum required for a successful simulation of advanced semiconductor devices. Additional physical models can be added to the simulation code to enchance the applicability to other advanced devices. Recall that a simulation program may not capture the workings of a physical phenomena unless the mathematical equations which govern it are included. Some of the enhancements which can be added to the standard hydrodynamic formulation include e€ects such as band gap narrowing, charge injection, acoustic charge transport, anisotropic materials, external circuits, photogeneration, photon recycling, thermionic emission, impact ionization, and tunneling, etc.. Sub-models may be required to complete the description of some of these models. For instance, photogeneration would require re¯ection, refraction,


A. W. Smith and K. F. Brennan

and absorption models. Boundary conditions can be added to simulate Schottky barriers, vacuum emission, pinned fermi levels, and charged oxide layers. The reader can consult the literature for the mathematical description for these and other models. Recall that due to the general numerical nature of the model presented to this point, the addition of new physical models is fairly straight forward. Source terms in addition to radiative(60), thermal(61), and Auger(62) and inverseÐAuger(63) generation and recombination mechanisms(61) such as photon recycling, photogeneration and tunneling have their mathematical formula added to the existing code for source terms and the partials with respect to the hydrodynamic variables are incorporated into the partial of the source term routine. The same is true for forcing terms such as the band gap narrowing term which acts as a convective term.

4.7. Material parameters The material parameters which are important for the successful hydrodynamic simulation of semiconductor devices are essentially the same as in the drift±di€usion equations with the addition of energy relaxation times and coupling terms to the lattice and hole systems. The energy relaxation times can be taken from constant ®eld Monte Carlo simulations(41) or ®xed values based on material parameters(14). Unfortunately, at present, these parameters are only known for a select group of materials with any accuracy (Si, GaAs). As the importance of other materials increases (GaN, InGaAs, . . . ) and the need for more accurate simulations arises more of these parameters will become available. In certain cases material parameters may be estimated by interpolating between values for the constituent binary materials for ternary compounds, or calculated based on other material parameters which are better characterized. For instance, the a or non-parabolicity factor, used in the ¯ux equations can be written in terms of known parameters of the semiconductor as(64)  2 me 1ÿ mo a1 …102† Eg Notice that the non-parabolicity factor will increase as the band gap or the e€ective mass of the band decreases. Narrow band gap materials such as InAs would have a very high a value (06 to 10) and therefore will show more of an e€ect in the simulation. One other set of parameters which is important to consider is the density of states of the electrons and holes. The density of states is a function of the respective carrier energy, or equivalently, its temperature. Therefore, the density of states must be re-calculated at each step of the simulation where the temperature changes. This dependency is a result of the transformation used to make the distribution function integrable with the density of states. This dependency is typically neglected in most hydrodynamic models. As just stated the material parameters which are important for the successful simulation of the devices are the non-parabolicity factors and the energy relaxation times. For silicon using a bandgap of 1.124 eV and an e€ective mass of 0.326, the non-parabolicity factor can be estimated, from Eqn. (102), to be 0.4039 eVÿ1. The low ®eld mobility values for Si are 1,332.2 and 380.6 (cm2/V/sec) for doping densities of 2.0  1015 cmÿ3 and 5.0  1017 cmÿ3, respectively. The energy relaxation time in silicon is set to 0.2 ps for the electrons and 20 fs for the holes. Using GaAs with a bandgap of 1.424 eV and an e€ective mass of 0.070 the non-parabolicity factor is estimated to be 0.60736 eVÿ1. The low ®eld electron and hole mobility values used for GaAs are 7,940.9 and 2,972 (cm2/V/sec) respectively for the doping densities similar to Si. The energy relaxation time in GaAs is approximately 0.1 ps for both

Hydrodynamic simulation of semiconductor devices


electrons and holes at these doping concentrations. These material parameters are used in the simulations presented in the next sections. 5. A P P L I C A T I O N O F T H E M O D E L

The model described in the previous sections can be applied to various devices. In this chapter it will be applied to two di€erent device structures as a means of illustrating the method and to contrast the di€erence between the hydrodynamic and drift±di€usion methods. The ®rst device examined is the standard device for testing hydrodynamic simulations, the ballistic diode. The ballistic diode is a single carrier, one-dimensional device. Being that the material parameters and properties of GaAs and Si are relatively well known, here we examine the behavior of a ballistic diode made from either GaAs or Si. The one dimensionality and unipolar nature of the device makes the model simple in that only one energy equation needs to be solved. For this reason, the ballistic diode was the ®rst device simulated by a hydrodynamic simulation program(26). In order to provide a complete solution, the simulation of the hole and lattice temperature equations will be included in the simulation of the ballistic diode presented here. The second device to be presented is a heterojunction bipolar transistor (HBT). The bipolar nature of the HBT necessitates the use of both carrier energy equations. In addition, the structure of the device is such that its simulation must be two-dimensional in order to capture physical e€ects such as current spreading at the emitter/base junction or current crowding at the base contact. The other e€ect of interest in the HBT simulation is the energy overshoot in the base region due to the heterojunction discontinuity. A description of each device is followed by the simulation results and a fairly concise explanation of the physics involved. 5.1. Silicon ballistic diode Figure 6 displays the structure and doping concentrations of the Si ballistic diode considered here. The device consists of a lightly doped n-type region sandwiched between two heavily doped n-type regions. The length of the lightly doped region is relatively short to show the e€ect of carrier heating. Electrons are injected from the highly doped region to the left and accelerated by the electric ®eld, Fig. 7, to the anode. The overall bias assumed to be applied to the ballistic diode for these calculations is 1 volt. In the lightly doped region the electric ®eld is substantial and acts to accelerate the carriers. Due to the

FIG. 6. Structure and doping concetrations used in the ballistic diode simulations.


A. W. Smith and K. F. Brennan

FIG. 7. Electric ®eld in a silicon ballistic diode baised to 1 V, this is very similar for all of the simulations (hydrodynamic, drift-di€usion, with and without non-parabolicity e€ects).

FIG. 8. Electron energy for the silicon ballistic diode biased to 1 V for the drift-di€usion with parabolic bands, hydrodynamic with parabolic bands, drift-di€usion with non-parabolic bands, and hydrodynamic with non-parabolic bands.

Hydrodynamic simulation of semiconductor devices


combination of the high electric ®eld and relatively short length, the carriers transit the lightly doped region quickly. Consequently, they su€er few scattering events and therefore gain more energy from the ®eld than they lose via lattice scatterings. As the carriers enter the highly doped region the electric ®eld suddenly drops. As a result, the acceleration dramatically decreases. Coupled with the enhanced scattering rate from the large concentration of ionized impurities, the carriers thermalize. Figure 8 shows the calculated electron energy pro®le of the device using both a drift± di€usion simulation and the hydrodynamic model. In addition, the result including the e€ect of the non-parabolicity factor within both simulation types is included in the ®gure. The energy is taken as the sum of the kinetic energy, 12 mv2, and the thermal energy, 3/2 kT. In the drift±di€usion simulation, only the kinetic energy changes; the carrier temperature is always kept equal to the lattice temperature. In contrast, within the hydrodynamic simulation the carrier temperature is determined during the course of the simulation. In the case of the drift±di€usion simulation the energy is essentially constant throughout the simulation domain for both the parabolic and non-parabolic formulations since the kinetic energy term changes very little. The most substantial di€erence is in the thermal energy of the carriers as re¯ected by the hydrodynamic simulation results. The energy pro®le calculated using the hydrodynamic simulation displays a pronounced peak near the anode section of the device. Notice also that the energy peak of the parabolic hydrodynamic simulation is larger than the non-parabolic case. From the form of the ¯ux equations (carrier and energy) it can be seen that the non-parabolicity term will tend to decrease the current, velocity, and e€ective mobility. Figure 9 displays the velocity pro®les calculated for the same situations, i.e., drift±di€usion and hydrodynamic both with and without non-parabolicity e€ects. These ®gures display quantities that re¯ect the internal characteristics of the device but which cannot be readily

FIG. 9. Electron velocity for the silicon ballistic diode biased to 1 V for the four di€erent simulation conditions.


A. W. Smith and K. F. Brennan

FIG. 10. Current voltage characteristics for the silicon ballistic diode utilizing the four di€erent simulation conditions: drift-di€usion with parabolic bands, hydrodynamic with parabolic bands, drift-di€usion with non-parabolic bands, and hydrodynamic with non-parabolic bands.

experimentally measured. The important information which can be measured are the output characteristics. Figure 10 shows the calculated current voltage characteristics of the Si ballistic diode. The current is lower for the hydrodynamic simulations as compared to the drift±di€usion model. This is due to enhanced scattering as the carrier energy increases which is not accounted for in the drift±di€usion simulation. The same trend is observed for the simulation which includes the non-parabolic terms. The non-parabolicity restricts the number of allowed energy levels which are available for transport. The carriers are therefore scattered more lowering the observable current. The variables in the hydrodynamic model are the carrier and lattice temperatures. As discussed, the total energy shown in the ®gures duscussed above is based on these temperatures as well as the kinetic energy. In order to more clearly assess the behavior of the device the electron, hole, and lattice temperatures for the Si ballistic diode are depicted in Fig. 11 for the case where the hole temperature has been assumed to be equal to the lattice temperature at the injecting contact. In these calculations, the energy bands are assumed to be parabolic and the applied bias is taken as 0.4 V. In addition, the boundary condition on the energy is ®xed at the anode; the carrier energies are forced to their equilibrium values at the contact. Comparison of Figs 11 and 7 shows that the peak in the electron temperature is not coincident in position with the maximum value of the electric ®eld; i.e. the carriers attain their maximum temperature prior to experiencing the peak in the electric ®eld. A similar e€ect is observed for the electron energy as shown by Fig. 8. This di€erence in the maximum in the ®eld position and carrier energy is due to the choice of boundary condition as stated above as well as the fact that the electrons enter a highly doped region after traversing the device.

Hydrodynamic simulation of semiconductor devices


FIG. 11. Calculated temperatures for the silicon ballistic diode with parabolic bands baised to 0.4 V, using a ®xed temperature boundary condition for the holes (minority carrier) at the anode end.

Rudan et al.(32) have observed a similar e€ect. They attribute the di€erence in the positions of the maximum temperature and ®eld to the fact that there is a back di€usion of the electrons from the highly doped region which acts to counterbalance the drift component of the current to some extent. As a result, the peak electron energy is o€set from the peak ®eld. The depression in the electron temperature prior to the reverse electric ®eld arises from the fact that the electrons lose energy going up the potential barrier. The change in lattice temperature in Fig. 11 with position is very small (<0.1 K) and is not noticeable on this scale. The holes lose energy near the injecting contact due to cooling by the electric ®eld, similar to the behavior of the electrons at the other end of the device. There is some degree of uncertainty for the boundary condition on the hole temperature. We have identi®ed two possibilities. These are that the hole temperature returns to the lattice temperature at the contact or that no minority carrier energy ¯ows into the device from the contact. Currently, little information is available to decide which of these two possibilities presents the most consistent picture of hole transport. The ®rst choice is that the hole temperature rises to the lattice temperature at the contact location due to the boundary condition imposed on the hole energy ¯ux equation. The corresponding calculated hole temperature is shown in Fig. 11. As can be seen from the ®gure, the hole temperature is ®xed at the lattice temperature at the contact located at 0.6 mm. The hole temperature then dips below the lattice temperature due to the retarding potential barrier. Finally, the hole temperature recovers to the lattice temperature after the peak in the potential barrier. If the boundary condition for the hole energy ¯ux equation is changed to the case where no minority carrier energy ¯ows into the device from the contact then the temperature pro®les shown in Fig. 12 are the result. In the carrier temperature calculations shown in Fig. 12, the boundary condition on the electron temperature at the contact remains the same. These may be compared to the temperatures for the ®xed boundary condition shown in Fig. 11.


A. W. Smith and K. F. Brennan

FIG. 12. Calculated temperatures for the silicon ballistic diode with parabolic bands biased to 0.4 V, using a zero ¯ux energy boundary condition for the holes (minority carrier) at the anode end.

Switching to the ¯ux boundary condition for the holes makes no di€erence in the output characteristics of the device since the output characteristics are dominated by the electrons. However, it does produce a di€erent temperature pro®le for the holes. The holes remain cold in the region of the contact and do not return to the lattice temperature until after the maximum in the ®eld is reached. This last ®gure shows some of the ambiguity of the hydrodynamic and other models of the same type, i.e. a mis-applied thermal boundary condition can lead to di€erent results just as a mis-applied symmetry boundary in the electrical model can lead to erroneous results. Unfortunately, little work has been done on the derivation of appropriate thermal boundary conditions for minority carrier current ¯ow and the correct choice of boundary condition cannot be ascertained from the output characteristics alone. Therefore, at this time, the correct choice for the hole carrier energy boundary condition is unknown. Switching the doping type of the ballistic diode from n-type to p-type produces similar temperature plots with the electron and hole temperatures switched. However, the overall magnitude of the hole temperature in a p-type device is lower than the corresponding electron temperature in the n-type device. The reason for this is the hole e€ective mass is higher and more tightly coupled to the lattice. Therefore, it is harder to raise the hole temperature to levels similar to the electron temperature. 5.2. Gallium arsenide ballistic diode The structure and doping concentration of the GaAs diode are the same as for Si, except for the doping concentration. The current voltage output characteristics of the diode are shown in Fig. 13. In the case of GaAs, the calculated device characteristics comparing the drift±di€usion to the hydrodynamic models are signi®cantly di€erent, as seen by Fig. 13. The non-parabolic formulation predicts a signi®cantly lower current than the parabolic

Hydrodynamic simulation of semiconductor devices

FIG. 13. Current voltage characteristics for the GaAs ballstic diode utilizing the four di€erent simulation conditions: drift-di€usion with parabolic bands, hydrodynamic with parabolic bands, drift-di€usion with non-parabolic bands, and hydrodynamic with non-parabolic bands.

FIG. 14. Electron energy for the GaAs ballistic diode biased to 1 V for the drift-di€usion with parabolic bands, hydrodynamic with parabolic bands, drift-di€usion with non-parabolic bands, and hydrodynamic with non-parabolic bands.



A. W. Smith and K. F. Brennan

FIG. 15. Electron velocity for the GaAs ballistic diode biased to 1 V for the four di€erent simulation conditions.

FIG. 16. Electron energy for the GaAs ballistic diode biased to 2 V for the drift-di€usion with parabolic bands, hydrodynamic with parabolic bands, drift-di€usion with non-parabolic bands, and hydrodynamic with non-parabolic bands showing the limitation of the Kane formulation.

Hydrodynamic simulation of semiconductor devices


formulation for the hydrodynamic model. In fact, the curvature of the current-voltage characteristic is di€erent between these two cases at high applied bias. The energy and velocity plots for the GaAs diode biased at 1 volt are shown in Figs 14 and 15. Unlike Si, the GaAs diode shows a large di€erence in the velocity throughout the low doped region between the hydrodynamic and drift±di€usion simulations, Fig. 15. There is also a trend towards saturation of the velocity in the middle of the undoped region. At the anode end of the device using the parabolic formulation a very strong velocity overshoot is observed for the hydrodynamic simulation. Though velocity overshoot is again predicted by the hydrodynamic model with the non-parabolic band formulation, it is not as pronounced as for the parabolic case. The energy plots in Fig. 14 qualitatively resemble the results calculated for Si, but have some important di€erences. First, near the source end of the device, the non-parabolic forms predict a greater energy than the parabolic case. At the end of the device, the trend is the same as in silicon; the parabolic case exhibits a higher energy than the non-parabolic case. However, the di€erence is greater between the two cases in GaAs than in silicon. The electrons are heated as they reach the anode end of the low doped region. The energy then falls to the lattice energy just after the high doped region is reached. A small distance is required to allow the carriers to scatter and lose their energy, relaxing back to the lattice energy. The current voltage characteristic of Fig. 13 shows the drawback of using the Kane formulation for band non-parabolicity in the hydrodynamic model. The current actually decreases as the voltage increases for this formulation. The reason for this trend is that the model has been applied to an energy range beyond the limit of the assumptions made during the derivation of the model at high biases. The binomial approximation was used in the derivation of the closed form coecients of the model. Azo€ (38) predicted that the binomial expansion of the di€usion term becomes inaccurate when the condition 2aW < 1 no longer holds. However, when non-parabolicity e€ects are included in the ®eld term the constraint is magni®ed by a factor of two, 4aW < 1. For the material parameters discussed here this limits the energy range to 0.4115 eV. Under a 1 volt bias, the peak average energy for either the parabolic or non-parabolic formulations is predicted to be signi®cantly less than 0.4115 eV. Therefore, one would expect that the non-parabolic model would apply. However, at a higher bias, the total energy in the GaAs device increases. Figure 16 shows the total energy of the GaAs ballistic diode biased to 2 volts. The energy peak of the Kane formulation is less than 0.4115 eV. However, the parabolic formulation reaches a peak value of 00.7 eV showing a large discrepancy between the two models. Previous work, using another band non-parabolicity formulation(59, 65), tracks the parabolic model rather than the Kane model. As a result, the predicted value of the Kane, non-parabolic model is somewhat suspect under these conditions since it fails to follow the parabolic case and alternative non-parabolic band formulation based models. 5.3. Heterojunction bipolar transistor, HBT, simulation As previously stated the HBT utilizes both carriers to perform its function. As such, the energy balance equations must be solved for both electrons and holes. However, the same trends hold as in the ballistic diodes already presented. That is, non-parabolicity and energy e€ects tend to reduce the output characterisics of the simulated device. This is true because the energy and non-parabolicity e€ects tend to reduce the e€ective mobility and di€usivity of the material, thereby lowering the calculated current. In the calculations discussed below, only the parabolic model results will be presented. This is because the energies to which the carriers are heated are, for the most part greater than the limit to which the Kane non-parabolicity model is valid. The ®rst HBT to be simulated is an AlGaAs/GaAs device with a fairly thick base region, 0.1 mm thick, of heavily p-type doping concentration 1.0  1019 cmÿ3, as shown in Fig. 17.


A. W. Smith and K. F. Brennan

FIG. 17. Device structure for the thick base heterojunction bipolar transistor (HBT).

Inspection of Fig. 17 clearly indicates the two-dimensional aspects of the device. Notice that the lateral dimensions of this structure are relatively small; the emitter, base and emitter±base separation regions are all 0.3 mm in length. The calculated results using both the drift±di€usion and hydrodynamic simulations will be presented for two di€erent bias conditions. In both cases the base±emitter voltage is varied at ®xed collector±base voltage.

FIG. 18. Gummel plot for the thick base HBT at two di€erent collector biases (2 and 5 V).

Hydrodynamic simulation of semiconductor devices

FIG. 19. Electron concentration for the thick base HBT at 2 V collector bias and 1.3 V base emitter bias. Notice due to the small lateral dimension the electron concentration is not dictated by the doping.

FIG. 20. Electron temperature for the conditions in Fig. 19. Notice that the electron temperature is not high in the base due to the heavy doping and the large dimension.



A. W. Smith and K. F. Brennan

FIG. 21. Hole temperature for the conditions in Fig. 19. The hole temperature does rise in the base-collector junction highlighting the need for inclusion of this e€ect in the simulation.

FIG. 22. Electron current pro®le for the conditions in Fig. 19 showing the two-dimensional e€ects of the device geometry.

Hydrodynamic simulation of semiconductor devices


The ®rst case is with the collector±base voltage ®xed at 2 volts and the second is a 5 volt collector±base bias. Figure 18 shows the calculated Gummel plot, showing drift±di€usion and hydrodynamic calculated base and collector currents as a function of the base±emitter voltage, for the two collector±base biasing conditions. The top four curves plotted in Fig. 18 correspond to the case with a 5 volt collector±base bias, while the lower four curves apply to the 2 volt base±collector bias. There is little di€erence in the calculated base and collector currents using the drift±di€usion and the hydrodynamic simulations. Under either base±collector biasing conditions, the drift±di€usion and hydrodynamic based calculations track one another quite closely. As we will see below, this is due principally to the relatively wide base width of the present structure wherein the e€ects of non-stationary transport are not substantial. Additionally, since the current mainly re¯ects the behavior of the mean of the distribution rather than its tails, the drift±di€usion and hydrodynamic methods produce comparable results. Figures 19±22 show the hydrodynamic calculated electron concentration, temperature, total energy, hole temperature and the electron current ¯ow pattern at a 1.3 V base±emitter and a 2 volt collector±base bias. The electron concentration behaves as expected; the electron concentration within the emitter and collector regions is high while it is very much lower in the base region and in the depletion regions. Notice that in this case, the electron concentration within the base is signi®cantly larger than its equilibrium concentration. This is due in part to the relatively small lateral dimensions of the device as well as the applied biases chosen. Since the lateral dimension of the base is relatively small, the injected

FIG. 23. Electron concentration for the thick base HBT at 2 V collector bias and 1.9 V base emitter bias.


A. W. Smith and K. F. Brennan

electron concentration from the emitter spreads throughout the base. The increase in the electron temperature near the heterojunction interface is clearly shown in Fig. 20. Inspection of Fig. 20 shows strong electron heating within the base±collector depletion region due to the high electric ®eld present. The electron temperature exceeds 5000 K in this region. Notice that the electrons and holes both are heated in the base collector region as shown by Figs 20 and 21. Most hydrodynamic simulations do not include the hole temperature equation, instead opting to apply the drift±di€usion equation for the holes. As noted, the hole temperature does change so ignoring this e€ect may lead to errors in the simulation. The electron current ¯ow pattern for the HBT structure of Fig. 17 under 1.3 V base± emitter and 2.0 V collector±base bias is shown diagrammatically in Fig. 22. The importance of the two-dimensional simulation is apparent from this ®gure. Notice the degree to which the electron current components are two-dimensional throughout the entire structure. The electron concentration and temperature calculated using the hydrodynamic simulation are illustrated in Figs 23 and 24 for the HBT structure shown in Fig. 17 at a 1.9 V base± emitter bias and a 2.0 V collector±base bias. The increased base±emitter bias acts to reduce the potential barrier between the emitter and base region resulting in a substantial increase in electron injection into the base. This can be clearly seen by comparing the electron concentrations under 1.3 and 1.9 V base±emitter bias as shown by Figs 19 and 23. The greater electron concentration within the base under higher base±emitter bias is to be expected since the larger base±emitter bias reduces the potential barrier at the emitter±base junction which in turn enables more electrons to overcome the potential barrier. By

FIG. 24. Electron temperature for the conditions in Fig. 23. Notice in this case the temperature is depressed as compared to Fig. 20.

Hydrodynamic simulation of semiconductor devices


reducing the potential barrier to injection, much colder electrons can now enter the base than before. As a result, the average electron temperature in the depletion layer is lower at 1.9 V than 1.3 V base±emitter bias as re¯ected by comparison of Figs 20 and 24. It is interesting to notice that in this device structure, at either base±emitter bias, the electron temperature within the base is close to that in the emitter, little e€ect of the injection condition on the energy is observed. Inspection of both Figs 20 and 24 shows that the energy within the base region is relatively low. The electron temperature is changed substantially mostly by the action of the collector±base bias. Since the electron temperature within the base is close to its equilibrium value, little di€erence between the predicted currents from the drift±di€usion and hydrodynamic simulations should occur, which is precisely what is observed. Next we examine the e€ect of increasing the collector±base bias to 5 V. The electron concentration for the device operated at a collector±base bias of 5 V and a base±emitter bias of 1.3 V is plotted in Fig. 25. As the collector bias is increased, the electrons are heated more at the base-collector junction. As can be seen from Fig. 26, the electron temperature is substantially higher in the base±collector region at higher collector bias as expected. At this bias condition there is not much di€erence in the output characteristics comparing the drift±di€usion to the hydrodynamic solutions. This is due to the ®eld dependent mobility in the drift±di€usion closely matching the e€ect of the energy scattering of the hydrodynamic model. The calculations were repeated for the case where the base±emitter is biased to 1.9 V and the base±collector is biased to 5.0 V. As before, the higher base±emitter bias

FIG. 25. Electron concentration for the thick base HBT at 5 V collector bias and 1.3 V base± emitter bias.


A. W. Smith and K. F. Brennan

FIG. 26. Electron temperature for the conditions in Fig. 25. Notice in this case the temperature is quite high but does not lead to more current ¯ow in the device.

FIG. 27. Device structure for the thin base heterojunction bipolar transistor (HBT).

Hydrodynamic simulation of semiconductor devices


leads to increased electron injection into the base region of the device. Again, the resulting electron temperature is lower with higher base±emitter bias. The second HBT examined has a thin base, 0.05 mm, and is lightly doped, 1.0  1018 cmÿ3 to show base pushout. A sketch of the device structure is shown in Fig. 27. Notice that this structure is signi®cantly larger in the lateral dimension than that shown in Fig. 17. The di€erent structure was chosen to enable comparison with the work of Yokoyama(66). Our calculations for the current agree fairly well with those of Yokoyama. Exact agreement is not possible since we do not know the full details of the parameters used in his model. Nevertheless, the reasonably close agreement provides some level of validation for our model. The biasing scheme is slightly di€erent than the ®rst HBT simulated in that the emitter collector bias, Vec, is ®xed at 2 V and the base±emitter voltage is varied. The Gummel plot for the second HBT structure is shown in Fig. 28. In this case, all three currents, base, emitter and collector, are plotted as a function of the base±emitter voltage calculated using both the drift±di€usion and hydrodynamic models. The collector voltage is held ®xed in all of these calculations at 2.0 V. The dashed lines shown in the ®gure correspond to calculations made with the drift±di€usion model while the solid lines represent the calculations using the hydrodynamic model. The unusual peak in the base current calculated using the drift±di€usion model which appears at low base±emitter voltage is apparently due to numerical imprecision. At these low biases, the value of the calculated current is less than the tolerance allowed in the solution. Consequently, it is likely that there exists substantial numerical error in this range. At higher biases, the behavior of the currents is better. Notice that at higher biases, that the sum of the currents obeys the node law in both models as it should. It is interesting to note that there exists a much more noticeable di€erence between the currents calculated using the drift±di€usion and hydrodynamic models for the second HBT device. The primary reason why there is

FIG. 28. Gummel plot for the thin base HBT, collector, emitter, and base currents are included for both the drift-di€usion and hydrodynamic models.


A. W. Smith and K. F. Brennan

FIG. 29. Electron concentration for the thin base design with a 2 V collector bias and 0.6 V baseemitter voltage. In this case the electron concentration is very close to the equilibrium concentration in the base region, in contrast to the thick base design.

FIG. 30. Electron temperature for the thin base design biased as in Fig. 29. Notice that the electron temperature is high throughout the base region of the device.

Hydrodynamic simulation of semiconductor devices


now a di€erence in the calculated currents is due to the fact that the second HBT device is much more susceptible to non-stationary transport e€ects. Figures 29 and 30 again show the electron concentration, and temperature within the device structure. Inspection of Fig. 29 shows that the electron concentration is very low, near its equilibrium concentration in the base near the base contact. This is di€erent from the ®rst HBT device shown in Fig. 17. In the ®rst HBT device, the lateral dimensions are much smaller so that the injected electrons from the emitter can spread throughout the base region even at relatively low base±emitter bias. In the second HBT device, since the lateral dimension of the base is 4.0 mm, the injected electrons from the emitter do not fully spread throughout the base, resulting in a very steep gradient for the electron concentration in the base. As in the ®rst HBT simulation the electrons and holes are heated at the base-collector junction. Perhaps the most interesting feature shown by the HBT analysis is that the electron temperature remains quite hot in the thin base device while in the thick base structure the electrons thermalize in the base. This can be readily seen by comparing the electron temperature plots for the two cases. As a result of the much higher electron temperature in the base, the overall electron temperature throughout the collector depletion region is higher. Comparing the calculated currents in both devices to that calculated using the drift±di€usion method, indicates a much greater di€erence between the two models in the narrow base device than the wide base device. Consequently, usage of the hydrodynamic simulation is of greatest importance in the narrow base device where the e€ects of carrier heating are most pronounced.

FIG. 31. Electron temperature for the thin base design biased to 2 V on the collector and 1.6 V on the base-emitter. Notice that the temperature is depressed under the base contact and is much lower overall as compared to the lower bias condition.


A. W. Smith and K. F. Brennan

It is interesting to examine the e€ect of increasing the base±emitter bias in the second HBT structure at a constant collector bias of 2.0 V. The resulting electron temperature within the device under these conditions is shown in Fig. 31. Comparing the electron temperature at the higher base±emitter bias to that shown in Fig. 30 for the lower base±emitter bias reveals that the electron temperature is substantially lower throughout the structure. This can be explained by consideration of the potential plots for the device. Figures 32 and 33 show the electrostatic potential within the device for both 0.6 and 1.6 base±emitter biases at ®xed collector bias. Comparison of the two plots shows that the higher base±emitter bias results in a smaller potential drop across the entire device. As such the device is closer to ¯at band leading to reduced carrier heating within the base±collector depletion region. Finally, the relationship between the electric ®eld pro®le and the carrier temperature within the HBT is examined. The calculated carrier temperatures and the corresponding electric ®eld pro®le for the thin base HBT structure biased at a base±emitter voltage of 0.6 V and a collector bias of 2.0 V are shown in Fig. 34. As can be seen from the ®gure, in this case, the electron carrier temperature peaks well after the peak in the electric ®eld within the collector region in contrast to what was observed above for the ballistic diode. The principal di€erence in this case is that the collector bulk region is very much longer and has a much lower doping concentration than the anode of the ballistic diode. As such, the boundary condition and e€ects of electron di€usion do not strongly in¯uence the carrier dynamics; the collector current is drift dominated within the depletion region. It should be noted that the present hydrodynamic model does not fully recover the physics of the HBT structure. Since the present approach uses a single valley to model the band structure, it does not properly account for the full e€ects of ballistic transport and velocity

FIG. 32. Potential pro®le for the bias conditions in Fig. 29. There is a potential drop from the emitter to the base.

Hydrodynamic simulation of semiconductor devices

FIG. 33. Potential pro®le for the higher base-emitter voltage. There is little potential di€erence between the base and the emitter, the device is approaching a ¯at band condition.

FIG. 34. Electron temperature, hole temperature, and electric ®eld pro®le through the device from the emitter though the base into the collector. The pro®les show the temperature maximums are not aligned with the ®eld maximum.



A. W. Smith and K. F. Brennan

overshoot. How the single valley formulation limits the accuracy of the method, when applied to compound semiconductors, can be understood as follows. In a material such as GaAs, in which intervalley transfer greatly alters the velocity of the electrons, proper accounting for the onset of this transfer is essential in order to correctly predict the extent and magnitude of the velocity overshoot. In a single valley model, the carrier velocity is determined through an energy dependent mobility. Though this is far more capable of capturing non-stationary transport e€ects than simply using a ®eld dependent mobility, as is done in drift±di€usion, it is still somewhat inadequate. In this approach, there is no time dependence to the transfer; once a carrier attains a certain energy, its mobility and hence velocity is immediately ®xed to the mobility which corresponds to that energy. Physically, what happens is somewhat di€erent. A carrier does not immediately transfer into the satellite valleys once it has attained sucient energy. Instead, there is some time lapse until the electron su€ers a phonon scattering before it transfers into the secondary minima. In a single valley hydrodynamic simulation this time lapse is not included. In order to overcome this limitation, a two or more valley model would need to be implemented. In a two valley simulation, an energy dependent relaxation time governs the time in which the carrier concentration changes within the secondary minima. Using this approach, the extent of velocity overshoot and ballistic transport can be more accurately modeled.

6. C O N C L U S I O N S

In this paper, an introduction to hydrodynamic modeling of semiconductor devices has been presented. The method is ®rst generalized beyond the hydrodynamic model in order to provide a fundamental framework for developing a numerical model at any level of complexity. The basic issues involved in constructing the mathematical model are presented in generality. These are: (1) the determination of the basic ¯ux equations which govern the system, (2) discretization of the basic ¯ux equations with and without source terms, and (3) application of Newton's method to solve these equations within the simulation domain. This approach is then applied to the speci®c case of hydrodynamic semiconductor device simulation. Again, the basic ¯ux equations are derived and discretized. A system of algebraic equations is then obtained which are solved on a mesh within the device simulation domain. Application of the boundary conditions and the constituent relations in order to complete the solution are discussed. Sucient detail is provided to enable programming a hydrodynamic simulation code. Several numerical pitfalls are pointed out so the reader can avoid numerical instabilities. To illustrate the importance of the hydrodynamic method, two di€erent devices are simulated. The ®rst is the standard test device for hydrodynamic simulations, the ballistic diode. Both silicon and GaAs based ballistic diodes are investigated using the drift± di€usion and hydrodynamic simulations. Comparison of the calculated results using these two di€erent approaches reveals the importance of the hydrodynamic method. It is shown that the hydrodynamic approach captures physical phenomena that the drift±di€usion technique cannot, speci®cally non-stationary transport e€ects. Speci®cally, the hydrodynamic simulation shows the importance of velocity overshoot within the structures. The ballistic diode device displays carrier heating in electric ®eld regions and the di€erence between the peaks in the ®eld region and the carrier temperature. The e€ects of band non-parabolicity within the hydrodynamic model are also investigated. The band structure is modeled using the Kane non-parabolicity formulation and the resulting transport equations are modi®ed accordingly. Comparison of the Kane nonparabolicity formulation to the standard parabolic formulation reveals some important di€erences. The current calculated using the non-parabolic model is less than that

Hydrodynamic simulation of semiconductor devices


determined based on the parabolic model. Additionally, velocity overshoot, as predicted by the hydrodynamic model with the non-parabolic band formulation, is not as pronounced as for the parabolic case. However, some caution should be exercised when employing the Kane non-parabolicity formulation to hydrodynamic simulation. It is found that at modest carrier energies, roughly 00.4 eV for GaAs device simulation using the parameter set here, the validity of the simulation equations becomes suspect. The binomial approximation was used in the derivation of the closed form coecients of the model. The binomial expansion of the di€usion and ®eld terms becomes inaccurate at modest carrier energies ultimately invalidating the results. The second class of devices simulated is the heterojunction bipolar transistor (HBT). This device highlights the need for including the energy balance equations for both carriers. The simulation of this device also displays the need for a complete two-dimensional simulation of this particular device. The most interesting result from the HBT analysis is that the electron temperature remains substantially hotter in a thin base device than in a thick base structure. As a result, the electrons thermalize in the wide base device but not in the thin base structure. This results in important di€erences in the terminal characteristics, speci®cally the current voltage characteristic, when comparing the hydrodynamic to drift± di€usion simulations. REFERENCES 1. 2. 3. 4. 5. 6. 7. 8. 9. 10. 11. 12. 13. 14. 15. 16. 17. 18. 19. 20. 21. 22. 23. 24. 25. 26. 27. 28. 29. 30. 31. 32. 33. 34. 35.

Pantankar, S. V., Numerical Heat Transfer and Fluid Flow, Taylor and Francis, New York, 1980. Cook, R. K. and Frey, J., COMPEL, 1982, 1, 65±87. Selberherr, S., Analysis and Simulation of Semiconductor Devices, Springer-Verlag, New York (1984). Hansch, W. and Miura-Mattausch, M., J. Appl. Phys., 1986, 60, 650±656. Widiger, D. J., Kizilyalli, I. C., Hess, K. and Coleman, J. J., IEEE Trans. Elec. Dev., 1985, ED-32, 1092± 1102. Hess, K., Advanced Theory of Semiconductor Devices, Prentice Hall, Englewood Cli€s, N.J., 1988. McKelvey, J. P. Solid State and Semiconductor Physics, Krieger, Malabar Florida (1982). Reif, F., Fundamentals of Statistical and Thermal Physics, McGraw-Hill, New York (1965). Dabrowski, W., Progress Quantum Electronics, 1989, 13, 233±266. Kane, E. O., J. Phys. Chem. Solids, 1957, 1, 249. Rogalski, A. and Piotrowski, J., Progress Quantum Electronics, 1988, 12, 87±289. Mains, R. K., Haddad, G. I. and Blakey, P. A., IEEE Trans. Elec. Dev., 1983, ED-30, 1327±1338. Scharfetter, D. L. and Gummel, H. K., IEEE Trans. Electron Dev., 1969, ED-16, 64. Tang, T.-W., IEEE Trans. Electron Dev., 1984, ED-31, 1912. Basore, P. A., PC-1D Installation Manual and User's Guide Version 2.1, Iowa State University Research Foundation, Ames, IA, 1989. Miller, J. J. H., Eng. Comput., 1985, 2, 71±73. Miller, J. J. H., Mathematics and Computers in Simulation, 1986, 28, 1±11. Leonard, B. P. and Mokhtari, S., Inter. J. Num. Meth. Eng., 1990, 30, 729±766. Dennis, Jr, J. E., Schnabel, R. B., Numerical Methods for Unconstrained Optimization and Nonlinear Equations, Prentice-Hall Inc., Englewood Cli€s, N.J., 1983. Stratton, R., Phys. Rev., 1962, 126, 2002±2014. Blotekjaer, K., IEEE Trans. Elec. Dev., 1970, ED-17, 38±47. Higman, J. and Hess, K., Solid-State Elec., 1986, 29, 915±918. Gaur, S. P. and Navon, D. H., IEEE Trans. Elec. Dev., 1976, ED-23, 50±57. Alwin, V. C., Navon, D. H. and Turgeon, L. J., IEEE Trans. Elec. Dev., 1977, ED-24, 1297±1304. Adler, M. S., IEEE Trans. Elec. Dev., 1978, ED-25, 16±22. Cook, R. K. and Frey, J., IEEE Trans. Elec. Dev., 1981, ED-28, 951±953. Cook, R. K. and Frey, J., IEEE Trans. Elec. Dev., 1982, ED-29, 970±977. Cook, R. K., IEEE Trans. Elec. Dev., 1983, ED-30, 1103±1110. McAndrew, C. C., Singhal, K. and Heasell, E. L., IEEE Elec. Dev. Letts., 1985, EDL-6, 446±447. Szeto, S. and Reif, R., IEEE Elec. Dev. Letts., 1987, EDL-8, 336±337. Wang, C. T., Solid-State Elec., 1985, 28, 783±788. Rudan, M., Odeh, F. and White, J., COMPEL, 1987, 6, 151±170. Forghieri, A., Guerrieri, R., Ciampolini, P., Gnudi, A., Rudan, M. and Baccarani, G., IEEE Trans. Comp. Aided Design, 1988, 7, 231±242. Snowden, C. M. and Loret, D., IEEE Trans. Elec. Dev., 1987, ED-34, 212±223. McAndrew, C. C., Heasell, E. L. and Singhal, K., Semicond. Sci. Tech., 1987, 2, 643±648.

360 36. 37. 38. 39. 40. 41. 42. 43. 44. 45. 46. 47. 48. 49. 50. 51. 52. 53. 54. 55. 56. 57. 58. 59. 60. 61. 62. 63. 64. 65. 66.

A. W. Smith and K. F. Brennan

Azo€, E. M., Solid-State Elec., 1987, 30, 913±917. Azo€, E. M., COMPEL, 1987, 6, 25±30. Azo€, E. M., J. Appl. Phys., 1988, 64, 2439±2446. Ou, H. H. and Tang, T. W., IEEE Trans. Elec. Dev., 1987, ED-34, 1533±1539. McAndrew, C. C., Heasell, E. L. and Singhal, K., Semicond. Sci. Tech., 1988, 3, 886±894. Feng, Y. K. and Hintz, A., IEEE Trans. Elec. Dev., 1988, ED-35, 1419±1431. Gardner, C. L., Jerome, J. W. and Rose, D. J., IEEE Trans. Comp. Aided Design, 1989, 5, 501±507. Quade, W., Rudan, M. and Scholl, E., IEEE Trans. Comp. Aided Design, 1988, 10, 1287±1294. Roberts, J. W. and Chamberlain, S. G., COMPEL, 1990, 9, 1±22. Shawki, T., Salmer, G. and El-Sayed, O., IEEE Trans. Elec. Dev., 1988, ED-37, 21±30. Shawki, T., Salmer, G. and El-Sayed, O., IEEE Trans. Comp. Aided Design, 1990, 9, 1150±1163. Zhao, Z. Y., Zhang, Q. M., Tan, G. L. and Xu, J. M., IEEE Trans. Comp. Aided Design, 1991, 10, 1432± 1440. Hrivnak, L., Progress Quantum Electronics, 1993, 17, 235±271. Kishore, R., Physica A, 1993, 196, 133. McAndrew, C. C., Heasell, E. L. and Singhal, K., Semicond. Sci. Technol., 1987, 2, 643. Woolard, D. L., Tain, H., Trew, R. J., Littlejohn, M. A. and Kim, K. W., Phys. Rev. B, 1991, 44(20), 11119. Ashcroft, N. W. and Merman, N. D., Solid State Physics, Holt, Rinehart and Winston, New York (1976). Higman, J. and Hess, K., Solid State Elec., 1986, 29(9), 915. Smith, A. W. and Rohatgi, A., IEEE Trans. Computer Aided Design and Integated Circuits, 1993, 12(10), 1515. Dresselhaus, G., Phys. Rev., 1955, 100, 580. Jain, S. C., Heasell, E. L. and Roulston, D. J., Progress Quantum Electronics, 1987, 11, 105±204. Parmenter, R. H., Phys. Rev., 1955, 100, 573. Smith, A. W. and Brennan, K. F., Solid State Elec., 1996, 39, 1055±1063. Smith, A. W. and Brennan, K. F., Solid State Elec., 1996, 39, 1659±1668. Schweps, R., Progress Quantum Electronics, 1996, 20, 271±358. Landsberg, P. T., Recombination in Semiconductors, Cambridge Univ. Press, Cambridge, U.K., 1992. Landsberg, P. T. and Beattie, A. R., J. Phys. Chem. Solids, 1959, 8, 73±75. Landsberg, P. T. and Robbins, D. J., J. Phys. C: Solid State Phys., 1977, 10, 2717±2739. Cassi, D. and Riccb, B., IEEE Trans Elec. Dev., 1990, 37(6), 1514. Stewart, R. A. and Churchill, J. N., Solid State Electronics, 1990, 33(7), 819. Yokoyama, K., Tomizawa, M. and Yoshii, A., IEEE Trans. Elec. Dev., 1984, ED-31, 1222±1229.