Formulation of Macroscopic Transport Models for Numerical Simulation of Semiconductor Devices

According to different assumptions in deriving carrier and energy flux equations, macroscopic semiconductor transport models from the moments of the Boltzmann transport equation (BTE) can be divided into two main categories: the hydrodynamic (HD) model which basically follows Blotekjer's approach [1, 2], and the Energy Transport (ET) model which originates from Strattton's approximation [3, 4]. The formulation, discretization, parametrization and numerical properties of the HD and ET models are carefully examined and compared. The well-known spurious velocity spike of the HD model in simple nin structures can then be understood from its formulation and parametrization of the thermoelectric current components. Recent progress in treating negative differential resistances with the ET model and extending the model to thermoelectric simulation is summarized. Finally, we propose a new model denoted by DUET (Dual ET)which accounts for all thermoelectric effects in most modern devices and demonstrates very good numerical properties. The new advances in applicability and computational efficiency of the ET model, as well as its easy implementation by modifying the conventional drift-diffusion (DD) model, indicate its attractiveness for numerical simulation of advanced semiconductor devices


I. INTRODUCTION
he semi-classical Boltzmann Transport Equaion (BTE), together with the Poisson equation, has been regarded basic and sufficient to model the electrical behaviors of modern semiconductor devices in spite of its inability to account for quantum interference and Heisenberg's uncertainty relations [5].However, it is extraordinarily difficult to solve the BTE in its explicit form, an integro-differential equation in the seven-dimensional phase space (r,k, t) OfqF (Of) Ot --cot, where r is the spatial coordinates, k is the wave vector, is time, f is the carrier distribution func- tion, q is the carrier charge, F is the electric field and the right hand side denotes the collision inte- gral, which contains the probability function of transfer between different states in k space.An alternative approach is to apply the Monte Carlo (MC) method [6]- [10] to obtain a stochastic solution for f(r,k, t).The ability to incorporate full band structures (and hence the hot-carrier effects at high energies) and detailed quantum mechanical scatter- ing rates makes the MC method very powerful in device modeling.Nevertheless, the MC method has many practical problems which limit its applications, such as, the need for extremely large computational resources, slow convergence in seeking self-con- sistent solutions with the Poisson equation in low- field and barrier regions, lack of rigorous calibration and verification of its large parameter set with ex- periments, and difficulty in treatment of generation- recombination processes and parameter sensitivity analysis and optimization [6, 9].By using the method of moments in the k space, we can simplify the BTE to a set of macroscopic transport equations.Theoretically, an infinite num- ber of moments will be necessary to regenerate all transport properties in the BTE.Yet it is practical to use only the first few moments to be computa- tionally efficient and physically meaningful.The ze- roth-to third-order moment equations can be generated from the moment variables: unity, the carrier velocity v, the carrier energy E and the carrier energy flow E-v.The corresponding macroscopic moment quantities are concentration n (or p), elec- tric current (carrier flux) J, total energy Wtota and energy flux S.
There is relatively little disagreement in the even-order moment equations, namely the continu- ity equations and the energy balance equations, which conserve the particle and energy respectively.For the odd-order moment equations of particle and energy fluxes, different assumptions and formula- tions have led to two main categories of models.The first model uses the same procedure as the even- order moments for the first moment, i.e., integration over k space is done with respect to v and an additional relaxation time is introduced for ((9f/t)col ).The moment generation is usually closed at the third moment with S expressed through some phenomenological laws.We denote this ap- proach as the hydrodynamic (HD) model [1, 2].
While various simplifications and formulations of the HD model exist in the literature (see, for exam- ple, [11][12][13][14][15][16][17][18][19][20][21][22][23]), the main assumptions are the momen- tum relaxation time and its parametrization and the closure of the third moment.The second model originates from Stratton [3, 4] and involves a very differ- ent approach.The basic formulation has been im- plemented for device simulation [24][25][26], although its assumptions and derivations have not been inves- tigated in detail.It is sometime treated as a varia- tion of the HD model with the same limitations [27].The fluxes, as odd moments of the BTE, only explic- itly depend on the odd part of the distribution function, fodd, (the parts with the even part, feven, vanish in integration due to symmetry).It is possible to make a microscopic relaxation time approximation to foaa as (0foaa/ 0 t )coll foaa / roaa and ex- press foaa in terms of feoe, and Zoaa [5].If we know the approximate functional dependence of feoen and Zoaa, closed forms of fluxes can be obtained.We denote this approach as the energy transport (ET) model [28].
After formulation of transport equations, there are coefficients to be parametrized in both models.We will look into the implications and influence of parametrization.The numerical properties of both models are then presented.The spurious velocity spike of the HD model in nin structures will be analyzed as a case study of influence from inappro- priate parametrization.Finally, we summarize the recent development and applications of the ET model on multiple band minima and nonuniform lattice temperature from Joule heating.We propose a new transport model, from which the solution pr6files of potential, electron and hole concentra- tions, and electron, hole and lattice temperatures can be obtained self-consistently to describe the interaction among the three subsystems: electron, hole and lattice.We denote this model as the DUal Energy Transport (DUET) model.

II. FORMULATIONS OF MACROSCOPIC EQUATIONS A. Basic Equations
Without further approximations, the even-order (zeroth and second) moment equations directly re- sult in the continuity equations and the energy bal- ance equations for electrons and holes: on q V. Jn -qU; (2) Ot Op + V-Jp -qU; (3)

O(p<Ep>)
Ot + V" St, F.Jp W (5)   where n and p are the electron and hole concentra- tions, which are also used as subscripts for other terms to denote the electron and hole parts, U is the net carrier recombination rate, and Wn and Wp are the net energy loss rates for electrons and holes.
The angular brackets represent averages over the distribution function and follow the convention in [3].
The terms U, IV, and W need to be carefully evaluated, considering all the important carrier and energy exchanging mechanisms, such as SRH re- combination, Auger recombination, impact ioniza- tion, and carrier scattering with phonons, ionized impurities, etc.According to different electron-hole interaction mechanisms, the dark net recombination rate can be written as Aug ii ii U--esrh -k-RAn ug + Rp a ap (6)   where Rsrh is the SRH recombination rate, ln ug and /pUg are the Auger recombination rates related to electrons and holes, and Gin and G are the impact ionization rates due to electrons and holes, respectively.The conventional expressions for these terms can be found in [29] and will not be discussed in this work.
The energy loss for electrons, W n, includes the following terms: 1. electron energy loss through scattering with phonons and others: Wla n(<En> Eo)lT-wn nkB(T ZL)/Twn (7) sidered the reverse procedure of the impact ionizations, the electron energy obtained can be estimated as Waug Egap RAn ug (11) The total energy loss for electrons can be as- sembled as Wn--Wla 4r-Wsrh --Wii-Waug ---Egap (Gff < ug ) (12)   The total energy loss for holes can be derived in like manner 3 3 Wp p +-kBZpRsrh rwp 2 --Egap(Gi RAp ug ) (13)   where Twn(Tn, TL) is the energy relaxation time for electrons; 2. electron energy loss due to SRH recombina- tions: Finally, the fluxes Jn and S are defined by an -qn(v) -q f d3kvf (14)   (8) 3. electron energy loss due to impact ionizations." For each electron related impact ionization, the electron energy loss is approximately Euap + 6nKnTp, while for each hole related impact ionization, electrons obtain an energy about 6pKnTn, where Ega v is the energy gap and the factors 6 n and 6p depend on electron and hole distributions, respectively.Therefore the elec- tron energy loss due to the total impact ioniza- tion is l/Vii (Egap + 6,,KnTp)Gin i-6pKBTnGj (9) S n n(Ev) fd3kEvf. (5)  Jv and Sp can be defined similarly.To obtain closed forms of J and S, there are mainly two approaches, namely, the HD and ET models.

B. Carrier and Energy Flux Equations
in the HD Model The first-order moment equation is obtained from multiplying both sides of (1) by v and then integrat- ing over k-space, For simplicity but without losing much accu- racy, W/i can be approximated by l/Vii EgapG (10) 4. electron energy gain due to Auger recombina- tions: Since Auger recombination can be con- qn(VkV)F + V(n(w)) -n v Here we will use the case for electrons as an illus- tration.The subscript n will be omitted when there is no confusion.The momentum relaxation time rp is defined by which is usually assumed to be a function of average energy (E), i.e., Zp rp((E)).Two additional as- sumptions are generally employed (for exceptions, see [21, 22]): 1.The conduction effective mass is constant ex- cept in the collision terms: E h2k2/(2m*) (18) where m* is the carrier conduction effective mass.
Equation (16) can now be written as S qrpV(nuu) q/xnF +/xV(nknT) (20)   where u -J/qn (v) is the mean electron ve- locity and /x is defined as qrp/m*.Notice that (20)   is nonlinear with respect to J due to the convective term.Also S can be expressed as , S= -knT + Ea +Q (21)   q where Ear m'u2/2 is the drift energy and the heat .flowQ is defined as m* Q --n((v u)(v u)2). (22)  In order to obtain a closed set of transport equa- tions, the third moment Q is usually identified with the thermodynamic heat flux and expressed phenomenologically using the Fourier's law where the thermal conductivity is estimated by the Wiedmann-Franz law k rc r -An/x (24) q where a is a proportionality factor.Equations ( 20), (21) and (22) represent the flux equations in the common HD models.
C. Carrier and Energy Fluxes in the Energy Transport (ET) Model The carrier distribution function f can be divided uniquely into its even and odd parts, denoted by feo and foaa, respectively, i.e., f feo + foad (25) and where leo(r, k) leo(r, -k) and fodd(r, k) --load(r, --k).Note that it is not necessary to assume foaa << feoen" Imposing separate microscopic relaxation time approximations reo(r, k) and roa(r, k) for feo and foaa, respectively, we have coil coil O coil leo fo foaa %v (27) since the equilibrium distribution function f0 is an even function.The relaxation time approximation in [3] is only a special case of ( 27), when rev roa z.
If both Zev and %a are assumed to be even functions of k, the odd terms in (27) give eL a ) Ot con roa(r,k). ( By substituting (25) and ( 28) into (1) and using the symmetry, we can obtain the relation between fev and fodd q%d F fdd= h Vkfe--'rdVVrfe" After substitution of ( 25) and ( 29) into ( 14) and (15), we can in general write J and S as n/2F + D. Vn + n---VT e (30)   0b ] nher + b e. Vn + VTe (31)  Ore where all transport coefficients in (30) and ( 31) are tensors and their components are given by e De 1 fd3kzoEvivfe < zoaEviv) (35)   Here, the subscripts i, j represent any combination of the coordinates x, y, z and ( ) denotes an aver- age over re, instead of f.The divergence of the tensors in (30) and (31) means that the divergence of each row of the tensor is taken as a vector element.
The equivalent carrier temperature T in (30) and (31) is defined to characterize the energy dependence of re," Note that T e has a different definition from T c in (19) of the HD model.The detailed derivation will be published elsewhere [39].An im- portant conclusion is that D/OT e is a proper driv- ing force for current, while OD/Ox is not.Note that, for a spatially independent relaxation time od(k) and thus spatially independent )(T e) and De(Te), (30) and (31) are equivalent to Eqs. ( 4) and (5) in [28].
The field-induced anisotropy of feo shows a re- markable influence for materials with small effective mass and needs to be included to achieve accurate device modeling.In the HD model, this effect is accounted for by the convective term qzpV. (nuu)   in (20).In the ET approach, this effect can be included using an anisotropic feo" In our previous work, we have developed a phenomenological method to extend an arbitrary isotropic distribution fi,o(E) to its corresponding anisotropic form fe.(k) [36] feo(k) 1/2[f/so(E(k + k0)) + fiso(E(k-ko))] ( where k 0 m* (v)/h.We can derive a nonisotropic ET model from (52) in a straightforward way.The detailed derivation will not be repeated here [36].
IIl. CRITICAL COMPARISONS The main assumptions of the HD model are the momentum relaxation time (17) and its parametrization, and the closure of the third moment by the Fourier and Wiedmann-Franz laws in (23) and (24).
The main assumptions of the ET model are the functional forms of fe. and rodd.These different assumptions have their respective advantages and limitations.Yet it is possible to express (40) and (41)  of the ET model in forms similar to [3] or the HD model through some reorganization and change of variables.The flux equations of both HD and ET models consist of contributions from the gradients of the potential (drift), the concentration (diffusion), and the carrier temperature (thermal diffusion).However, the thermal diffusion terms have different forms and the carrier temperatures have different definitions.The three flux terms can be understood from a phenomenological thermodynamic point of view [58], which is derived rigorously under the local equilibrium assumption.This assumption dictates that at a fixed position r and time t, any accessible state of electrons and holes can be uniquely charac- terized by a local temperature, T and Tp, respec- tively, and a quasi-Fermi level bn and bp.Related comparisons of the formulation of the current equa- tion regarding the diffusion term have been performed in [43].
We will present the parametrization and numeri- cal properties of the HD and ET models in the following.The well-known velocity spurious spike in submicron nin structures of the HD model will be used as a case study.

Parametrization
Although different assumptions and formulations of J and S are used, both HD and ET models assume specific functional dependencies of their transport coefficients.From the continuity and energy balance equations, the same set of transport coefficients, U, 'r w, Rsrh, Gii and R Aug, are required in the HD and ET models.From the flux equations, rp (or/x) and A are required in the HD model, while/x and C are necessary in the ET model.Sometimes these coeffi- cients cannot be directly obtained from experimental measurement, such as, when they are parametrized through the average energy or the carrier temperature.The most popular way of func- tional extraction is from bulk MC calculations, which theoretically can give a stochastic, yet exact solution of the BTE.As an illustration, in the HD model -p is first assumed to be a function of (E) only.Then from bulk MC calculations at different field strength, relations between (E) vs. F and between -p vs. F are obtained.Finally 'p can be expressed as a func- tion of (E) or T [11, 31].There are four question- able areas in this methodology: 1. MC calculations require a large set of parameters (usually over 60) such as the scattering rates and pseudopotentials [6].However, this large set of parameters is obtained mostly from theoretical calculation and is often merely cali- brated against experimental velocity-filed curves for bulk at different doping levels.Dif- ferent calibration methods for high-energy ef- fects have been proposed [44,45], but are still far away from unanimous acceptance.2. Monte Carlo calculations have built-in statisti- cal noise, which can be qualitatively physical, since current fluctuation for extremely low car- rier concentration is experimentally observable [5].In high carrier concentration regions, this is an artifact from small sampling numbers.This makes functional extraction from derivatives very inaccurate.Alternative approaches such as the scattering matrix method [32], which also gives an exact solution of the BTE, may be- come helpful in this respect.3. It is hard to justify the accuracy if parameters such as zp an r w, which are related to nonlocal and transient effects, are initially extracted from the bulk MC calculations and are directly applied to simulate devices with large field variations.Notice that the nonlocal transport characteristics built in the distribution function from MC calculations may be never utilized.Prescribed-field MC calculations, where field variations in space or time are either similar to those in the intended device [46, 47] or are from theoretical analysis and transformation [50,51], are required for the MC information of nonlocal effects to be really used in the macroscopic models.Unfortunately, even though smart prescribed-field MC methods are indeed used, we are still faced with the dilemma that the MC parameters related to the nonlo- cal effects are hard to be calibrated with experiments.
4.Last but not least, some of the parametrization or the functional formulation of transport coef- ficients is actually an ad hoc assumption and is subject to questionable accuracy and scaling.Again using 'p parametrized by (E) as the example, the parametrization can not be justi- fied for all scattering processes [37].We cannot expect at all that the physical effects included in bulk MC calculations will be reflected cor- rectly in the macroscopic models through this kind of ca4ibration.The prescribed-field MC calculations [47] can provide a necessary check on the parametrization.Although still not suf- ficient, it indeed shows the problems of the assumptions in such parametrization.More physical parametrization of Zp and z has been attempted [17], but it is based on the theoreti- cal work in [16], where it is assumed that the carrier distribution function f can be con- structed accurately by using the linear combi- nation of moments, n, J, W and S, i.e., f(r, k) n(r)g0(k) + J(r)gl(k) + Tn(r)g2(k) + S(r)g3(k).This assumption removes one of the most important advantages of the HD model: no prescriptions of f [1].It also lumps the streaming properties of tail and origin of the distribution which are known to be different.
Moreover, judging from the Maxwellian distri- bution functions where f is inversely propor- tional to Tn 3/2 exp(1/Tn), the assumption of f as a linear function of T n is rather poor.
The methodology of parametrization and calibra- tion of both macroscopic models needs more careful thought.There are good examples of parametriza- tion such as the low-field mobili.ty of MOSFET in the drift-diffusion (DD) model.The scaling relations [52] and the universal mobility [53] provide guide- lines for physical parametrization [54, 55].Similar philosophy should be applied to other transport coefficients.The MC calculations can still be used as a validation and calibration tool [49].However, while waiting for a better methodology, users have to keep in mind the assumptions and limitations behind present abstraction steps.

B. Numerical Properties and Discretization
For the HD model, the convective term in (20) will make the current equation become hyperbolic [57] and for critical cases, shock waves in carrier velocity may be observed [59, 60].Furthermore, due to the nonlinearity introduced by the convective term, it is not possible anymore to substitute J of the HD model into the continuity equations ( 2) and (3) as in the DD and ET models, where the current equation is linear with respect to J or u.This indicates that the current equation has to be solved explicitly by a hyperbolic partial differential equation (PDE) solver.
Many authors have ignored the convective term for the HD model, and did not solve the current equa- tion explicitly.However, there is no strong evidence that the convective term can be ignored under general device conditions.Matching a few experimental results by tuning parameters in the reduced HD model does not imply that the convective term is negligible.The field-induced anisotropy corresponding to part of the convective term in the HD model is actually quite significant in some situations [36].If the current equation is solved explicitly, an increase in the total number of equations also implies a severe decline in computational efficiency for simi- lar numerical schemes.Fortunately, the hyperbolic PDE solver is usually easier to parallelize for mas- sive computation.
For the HD model, proper upwinding [59] or sophisticated Essentially-Non-Oscillatory (ENO) discretization schemes of the flux equations may be necessary to resolve the numerical problems from possible shock-wave formation.For the reduced HD model (neglecting the convective term) and the ET model, the discretization scheme plays an impor- tasnt role in the stability of Newton-like methods.It is well-known that the .simplecentral difference discretization scheme for the current equation will generate oscillations or wiggles in the solution when the potential difference between grid points exceeds twice the thermal voltage [61], which will deteriorate both numerical accuracy and stability.The popular Scharfetter-Gummel [61], proper upwinding [62] or ENO [23, 40] discretization scheme is necessary.The numerical (artifical)diffusion introduced by Scharfetter-Gummel discretization does not affect the solution severely provided that the grid is rela- tively dense in space-charge regions [62].Moreover, special attention should be paid to the energy flux equation because of its complexity.For reduced HD models, since the parameter K is proportional to n and hence may change exponentially between mesh points, the discretization schemes are usually very complex and sometime even show numerical insta- bility owing to incorrect upwinding schemes [63][64][65][66][67]. On the other hand, the energy flux equation of the ET model has an elegant parallel to the Scharfetter-Gummel discretization scheme of the current equation [28].This implies that the ET model will be stable for Newton-like methods [42] and may be easily implemented by modifying exist- ing DD-based software.

C. The Spurious Velocity Spike---A Case Study
A spurious velocity spike, in addition to the physical velocity overshoot, was commonly observed in the HD model for both silicon nin and npn structures close to the anode junction [13, 23].It is interesting that for the same structure there is much less spurious spike in the ET model [28].This does not imply that the ET model is superior owing to other limita- tions.Nevertheless, based on the similarity between the HD and ET models, additional clues may be provided for the spurious spike from their compari- son.The spurious spike is not caused by neglect of physical mechanisms in the HD model such as the nonparabolic band structure or the carrier-carrier scatterings, since for the self-consistent ensemble MC calculations with these effects purposedly omit- ted, the spurious velocity spike is not observed [40].The Wiedmann-Franz law was suspected to be re- sponsible for the spurious spike since it was re- ported that the spike can be eliminated by phe- nomenologically reducing A in (24)  [13].However, it is found that different structures require various values of A to eliminate the spike [34].In [35,48], components from drift (V), particle diffusion (Vn)   and thermal diffusion (VT n) are compared, where is the electric potential.The spike is identified with the inbalance between particle and thermal diffu- sion terms.Reduction of A will only push the ther- mal diffusion velocity peak toward the anode junc- tion.Hence, the optimal A to avoid spurious spike is rather structurally dependent.In addition, phe- nomenological adaptation of a certain parameter to fit specific curves of different devices is risky, since more serious errors many be introduced in other device behaviors.By investigating the 1-D profiles of Q [30, 33], the Fourier law in (23) is shown to be very inaccurate, which implies that tuning the pro- portionality factor of (24) is not very meaningful.Explicit treatment of Q is suggested [33] for more accurate modelling.
From the MC simulation results, the formation of the spurious spike can be identified by two carrier populations with very different average energies [34].Close to the anode junction, a pack of hot carriers is injected into the sea of cold carriers in the n / region.The average energy of all carriers will drop sharply since the number of cold carriers grows almost exponentially.However, it takes a larger distance for the hot carriers to lose energy from scattering and hence two populations of hot and cold carriers coexit.The magnitude of the spurious spike can actually be estimated using the two-popu- lation model [40].This means that the distribution function and hence the transport properties can not be well characterized by (E) as postulated in [58].The work in [17], [33] and [48] supports this micro- scopic point of view by stating that either explicit treatment of Q or inclusion of J and S (possibly VT c) in the parametrization of zp will enhance accu- racy close to the anode junction.The ET model also fails to explain conditions with two carrier populations.Although the ET model does not show a spurious spike due to different assumptions that lead to a negligible thermal diffusion velocity, it does not mean that it is more accurate in a general way.The estimation of the spurious spike in [40] based on the difference between the ET and HD models is still valid simply because the velocity profile of the ET model is much closer to the exact solution from MC calculations around the anode junction.
The case of the spurious spike demonstrates the importan.ce of parametrization of transport coeffi- cients.A more well-known case is the impact ioniza- tion, whose transport contributions are mostly from the high energy tails of the distribution function [41,44].It has been shown that impact ionization is inaccurately modeled by F or (E) [56].tional efficiency.It is interesting to note that r m of /z in the HD model cannot be summed up directly as in the DD approach to account for the Gunn effects [71] unless r m or/x can be considered spa- tially slow-varying, which is not generally applicable for submicron devices [72].The composite mobility in the DD model, averaged over multiple valleys, will show negative differential resistance (NDR) if the valley with higher energy has smaller mobility.For example, the popular GaAs electron mobility in the DD model is expressed as [73]   Id, o "Ji-Usa ( ,F, F ;icrit ) /x(/zo, FII )

(53t
where/x 0 is the low-field mobility, Usa is the satura- tion velocity and Fcn is the critical field.However, the inclusion of NDR will require a very fine grid to resolve the Gunn-domain region where a large elec- tric field exists, and will contribute to numerical instability in typical DD simulations of low-doped regions.Nonphysical ripples in simulated I-V curves caused by Gunn domains hopping between the grid points are typically observed, although the magni- tude of the ripples can be reduced with finer grid [71, 74].On the other hand, it is possible to use a composite mobility for different valleys in the ET model similar to the DD model since in the dis- cretized flux equations only the product of n in- stead of n and /x separately appears in each term [38].In [74], it has been proposed to use (1 A)/.zt + (1 + A)/., H t( 0, Ze) For III-V materials with multiple conduction band minima and heterojunction applications, the HD model has been formulated and applied from the very beginning [1, 68, 69, 70].Each valley will con- tribute a new set of moment equations and interval- Icy scattering rates need to be parametrized.These additional equations will seriously degrade computa- where A tanh(aA(T -TA)) is a fitting function and a A and T A are fitting parameters./zc and /x H are defined by m( L) 110 + q id, O'wFc2rit -L 1 1+ q ld, oT.w Fc2r --1 3/2 (55) re) t'tO if" q Osat Tw Fc -L 1 1 + -qvs.t.,Fcri ---1 (56) By ignoring the diffusion current contribution in the energy balance equation, /XL is extracted from the constant mobility approximation and /zn from the constant drift velocity approximation of ( 53).This mobility model has been implemented in the ET model on a 1-/zm gate GaAs MESFET and has shown good numerical convergence.The ripples in the I-V curves are eliminated on a relatively coarse grid, on which the DD model with (53) demon- strates large ripples.This is possibly due to the smoother profile of T than that of F. The Gunn domain, a localized accumulation and depletion of electrons, can be clearly observed in the ET simula- tion near the gate-drain region [74].
For heterojunction applications, the implementa- tion of the ET model is very similar to that of the HD model [24, 26, 70].Both models will resolve the mobility degradation problems from built-in fields in the DD model since the mobility is now a function of the average energy instead of the local electric field [75].

B. Lattice Temperature Variations
Lattice temperature may significantly differ from the room temperature due to Joule heating for devices operating at high currents, especially for III-V materials and silicon-on-insulator (SO1) tech- nology due to low thermal conductivity.Exclusion of Joule heating will seriously affect the simulation accuracy, particularly for breakdown effects.In the HD model, the thermal diffusion equation has been adopted for simulation for years [19], since in the computational fluid dynamic regime (the predecessor of the HD model) the heat source has been taken into account for a long time.In the ET model, the thermal diffusion equation can be easily added CL Ot V" (KLVT) H (57) where C., x. and H are the heat capacity, the thermal conductivity and the net heat source for the lattice, respectively.Considering the total energy conservation for the carrier and lattice system, we can express H as 3n kn(T T) 3p kn(Tp-T) It= 2 'wn 2 Twp 3 + Tp) + Egap)Rsr h + F.J n + F.Jp (58) Equations (2-5), ( 30), ( 31), ( 57) and (58) constitute a new transport model, DUET (DUal Energy Trans- port), which accounts for all thermoelectric effects in most devices.
Simpler models can be derived from DUET by using limiting assumptions.For example, assuming that T L is constant, Dual ET becomes the previous ET model.With approximations T, Tp T L, the thermodynamic (TD) model [58] can be derived from DUET, which can handle cases involving lattice heating and thermal diffusions.The DUET model has been implemented in the 2-D device simulator PISCES [73].By using the generalized Scharfetter- Gummel scheme [28, 38] for discretization and the fully-coupled Newton method, excellent conver- gency behavior was observed in almost all cases.
Simulation results using the DUET model for a p n diode are shown below.The doping profile is plotted in Fig. 1.Profiles of T, Tp and T L for a forward bias of 1V are shown in Fig. 2. It can be seen that T n = Tp = T L under this forward bias and the hot-carrier effect is negligible.Hence, both TD and DUET models are valid.Since the tempera- tures are not uniform in space due to Joule heating from large forward current, the ET model is less accurate.The reverse-bias breakdown curves for this diode model at a reverse bias where the reverse current is 3 A/cm.using DD, ET and DUET models are shown in Fig. 3, where the field-dependent impact ionization model is used for DD, and the carrier temperature-which is critical to many designs, differs signifi- dependent impact ionization model is used for ET cantly.For the DD model, where neither carrier and DUET [77].The temperature model will still energy transport nor the lattice energy transport not be accurate, since impact ionization is affected (thermal diffusion) is included, Vbo is about 208 V. by carriers above .somethreshold voltage and has For the ET model, Vbo is reduced to 125 V by little dependence on the average energy, as we have including the carrier energy transport.For the mentioned before.Nevertheless, the following anal-DUET model, Vbo is only about 40 V.The distribu- ysis still demonstrates how the carrier and lattice tion of T n, Tp and TL under reverse biases are system interact through generation-recombination shown in Figures 4 and 5.At a low or intermediate terms.Breakdowns due to impact ionization occur current level, where T L is close to the ambient at reverse bias around 18-20 V.Moreover, a sharp temperature, thermal diffusion can be omitted.snap-back (burn-out) occurs at high current levels Therefore, ET as well as DUET can be applied.for all models.However, the burn-out voltage, Vbo, However, when the reverse current increases to a very high level, the lattice is heated by currents and T begins to increase (see Fig. 5).This results in a dramatic increase in the intrinsic carrier concentra- tion and hence the total current.This is similar to a positive feedback mechanism and the lattice tem- perature runs away.In this case, only the DUET model is applicable, since both the hot-carrier effect and lattice heating effect need to be taken into account.

V. CONCLUSIONS
We have presented and compared the derivation and assumptions of both HD and ET models.Impli- cations of the parametrization of transport coeffi- cient should be considered carefully for accurate device simulation.The numerical properties and the applicable domains of both models have also been examined.In view of the numerical efficiency, physi- cal applicability, and easy implementation by modi- fying well-developed drift-diffusion software, the DUET model is considered to be an attractive tradeoff for advanced semiconductor device simula- tion.

FIGURE 4 FIGURE 3 I
FIGURE 4 Temperature profiles generated by the DUET model at a reverse bias where the reverse current is mA/cm.

FIGURE 5
FIGURE 5 Temperature profiles generated by the DUET