Accounting for bandstructure effects in the hydrodynamic model: A first-order approach for silicon device simulation

Accounting for bandstructure effects in the hydrodynamic model: A first-order approach for silicon device simulation

Solid-State ElectronicsVol. 35, No. 2, pp. 131-139, 1992 Printed in Great Britain. All rights rcscmd 0038-l 101/92 S5.00 + 0.00 Copyright 0 1992 Pcrg...

808KB Sizes 0 Downloads 7 Views

Solid-State ElectronicsVol. 35, No. 2, pp. 131-139, 1992 Printed in Great Britain. All rights rcscmd

0038-l 101/92 S5.00 + 0.00 Copyright 0 1992 Pcrgamon Pmw plc

ACCOUNTING FOR BANDSTRUCTURE EFFECTS IN THE HYDRODYNAMIC MODEL: A FIRST-ORDER APPROACH FOR SILICON DEVICE SIMULATION T. J. BORDELON, X.-L. WANG, C. M. MAZIAR and A. F. TAKEI Department of Electrical and Computer Engineering, The University of Texas at Austin, Austin, TX 78712, U.S.A. (Received 7 May 1991; in revised form 12 August 1991) Aimtract-This paper presents a non-parabolic formulation of the hydrodynamic model for Si device simulation. By adopting a first-order approach, non-parabolic band effects are included as simple corrections to the conventional hydrodynamic model. These corrections are shown to improve the accuracy of the hydrodynamic model under submicron device conditions. The iirst-order approach captures the important effects of non-parabolicity in the conservation relations by introducing a minimal number of approximations. Unlike previously reported non-parabolic models, only a single empirical constant is required in addition to the relaxation times. An advantage of the formulation presented here is that the average electron energy is retained as a state variable, and an electron temperature need not be defined. Advanced hot-carrier models based on average electron energy will benefit directly from the improved accuracy of the model.



The hydrodynamic charge transport model[l] is an attractive alternative to the drift-diffusion model for routine analysis of submicron Si devices, offering increased physical accuracy without the computational burden characteristic of Monte-Carlo methods. The hydrodynamic model consists of a set of conservation equations for particle number, average momentum and average energy derived from the Boltxmann transport equation[2]. Unlike the driftdiffusion model, the hydrodynamic model is able to describe non-stationary transport phenomena, such as velocity overshoot. In addition, the average energy can readily be used to model hot carrier processes with greater accuracy than the local electric field[3]. The purpose of this paper is to present an improved formulation of the hydrodynamic model which more accurately models transport under submicron device conditions by accounting for the non-parabolic nature of the energy band structure. Commonly, a parabolic band assumption is invoked in order to facilitate the derivation of the conservation relations which comprise the hydrodynamic model. The conservation relations are, however, a&ted by nonparabolicity. For example (as will be shown), a consequence of first-order non-parabolicity is that diffusion is lower then that predicted by assuming a parabolic band. The approach adopted here introduces non-parabolic corrections to the conventional hydrodynamic model in a first-order sense. The flrstorder approach capture the salient features of nonparabolicity in the conservation relations with a minimum number of empirical parameters,

Formally, the derivation of the conservation equations is complicated substantially by including non-parabolicity. Higher-order moments of the distribution function, which are difficult to approximate, appear in the momentum and energy conservation equations[4]. The relationship between the carrier temperature (as defined by the variance of the velocity distribution) and the average energy also involves higher-order terms which can only be estimated with limited accuracy. In contrast, the parabolic band assumption yields a simple thermodynamic relationship between the average energy, carrier temperature and drift velocity[2]. The carrier temperature is a useful concept with a parabolic band since it results in the heat fhtx as the only higher-order unknown (which must be approximated). For a nonparabolic band, the carrier temperature is effectively an unknown quantity independent of the average energy and drift velocity, and the advantage of defining a carrier temperature is largely lost. Nonparabolicity therefore mandates additional approximations in order to close the set of transport equations with a tractable number of unknowns. Recently, several distinct approaches for including non-parabolicity in the hydrodynamic model have been presented. Thoma et a/.[51 have opted to formulate a velocity conservation equation (instead of a momentum conservation equation), and have defined the carrier temperature as a state variable. They circumvent the problem of establishing a relationship between the average energy and carrier temperature by defining the conservation equations strictly in terms of the carrier temperature. In order to accomplish this, a total of five relaxation times are needed, 131

T. J.



although their Monte-Carlo results indicate that two of these are essentially equivalent. A disadvantage of their approach is that conservation of average momentum is not insured by a velocity conservation equation with a non-parabolic band. However, it is not yet clear what distinctions exist between momentum and velocity conservation. Woolard et aL[6] have adopted a velocity conservation equation as well for modeling transport in GaAs. They introduce two temperature variables and define them as functions of the average energy via uniform field Monte-Carlo simulations. Stewart and Churchill[7] have also relied on an empirical, uniform field relationship between average energy and temperature for describing transport in GaAs. They utilize a momentum conservation equation and relate the average momentum to the drift velocity through an empirical energy-dependent effective mass derived from uniform field MonteCarlo simulations. By empirically relating the carrier temperature and average energy as well, they also preserve the parabolic band form of the momentum equation. A consequence of their approach is that the energy flux is also defined in terms of the empirical effective mass, and a heat flux contribution no longer appears in the energy conservation equation. This paper presents an alternate method for including non-parabolicity in the hydrodynamic model. By adopting a first-order approach, our objective is to minimize the uncertainty associated with the treatment of higher-order terms, as well as to reduce the effort required to implement the model by minimizing the number of empirical relationships which must be defined. By explicitly imposing a symmetry requirement on the wave-vector distribution function, a momentum conservation equation is derived which analytically accounts for non-parabolicity. The advantage of this approach is that a carrier temperature concept is not required. If the convective (or drift) component of energy is negligible, the symmetry condition is identical to the assumption of a scalar electron temperature in the parabolic band model. In addition to the energy relaxation time, one empirical parameter is required in the energy conservation equation. Since this parameter may be treated as a constant, the form of the energy equation closely resembles the conventional parabolic form.



multiplying the Boltzmann transport equation A(k) and integrating over R-space, yielding:

* (1) cdl In eqn (l), the angle brackets denote an ensemble average, E is the electric field, and u = V&k)/h is the group velocity [with r(k) as the energy]. Neglecting recombination and generation, the continuity equation is obtained from eqn (1) with A(k) = 1 as: V.(nv)=O,

2. I. Hydrodynamic model The hydrodynamic model consists of a set of equations describing conservation of particle number, momentum and energy, which are solved in conjunction with Poisson’s equation. A conservation equation for some quantity A(k) can be created by


where v =
t (nP) + V * (n(hku)) + qnE = -5, P

where p E (hk) is the average momentum. With A (k) = S(k), eqn (1) becomes the energy conservation equation: $nv)+V*(n

(4) w

where w = (8) is the average energy (w,, is the equilibrium average energy). Collisional dissipation is accounted for in eqns (3) and (4) by the momentum and energy relaxation times [r,(w) and r,(w), respectively] which, as suggested by ShutiS], can be defined as functions of average energy via uniform field Monte-Carlo simulations. Band structure effects are manifested in the hydrodynamic model through the relaxation times, the momentum and energy fluxes [n(hku) and n(uB), respectively], and through the relationship between the average momentum and the drift velocity. For a parabolic band, the momentum flux tensor takes the form: n(hk,uj) = nm*v,vj+ nk,T, 6,. (5) The electron temperature T, is related to the average energy via: w = jk, T, + $m*v*. (6) The parabolic band energy flux is: n





nkB T,v + Q,


where the heat flux Q is formally an additional unknown. In order to close the set of equations, the heat flux can be neglected or approximated by a phenomenological Fourier law: Q = -XVT,,


where the thermal conductivity X is given by a generalized Weidemann-Franz law[9].

Hydrodynamic model for Si device simulation

With the first-order a(1




(where m* is the band-edge effective mass) the relationship between momentum and velocity is: hk = m*(l + 2ad)u.


In contrast to the parabolic band case (a = 0), the average momentum depends upon the drift velocity and the average energy flux per carrier (udp). The non-parabolic relationship between energy and velocity is: S(u2) = & [(l - 2am*uz)-“r - 11.


Expansion of the square root in eqn (11) yields terms in ascending powers of velocity which are not negligible when averaged. These quantities are additional unknowns representing higher-order moments of the velocity distribution. The average energy is no longer related to the electron temperature (as defined by the variance of the velocity distribution) and the drift velocity through eqn (6). Because of the higher-order terms appearing in eqn (1 l), only an approximate relationship between the electron temperature and the average energy can be established with a nonparabolic band. The transport equations should therefore be formulated in terms of either the electron temperature or the average electron energy, but not both, in order to minimize physical error. Thoma et al. have opted to use the electron temperature as a state variable in formulating steadystate conservation equations. Alternatively, the approach presented here adopts the average electron energy as a state variable in order to formulate a time-dependent set of conservation equations, and avoids introducing an electron temperature entirely. Since degradation phenomena depend fundamentally on energy and not temperature, this approach is appealing for the purpose of applying the hydrodynamic model in conjunction with hot-carrier models which utilize average energy. Non-local hot-carrier models can be created which may depend upon either the electron temperature or the average energy by utilizing uniformjeld relationships between the average energy or temperature and the hot-carrier processes (e.g. impact ionization). A more sophisticated hot-carrier model, which would account for nonuniform effects, is likely easier to construct using the average energy since the energy distribution of electrons is more closely characterized by the average energy than by the electron temperature (which refers to a property of the velocity distribution). A uniform field relationship between the average energy and electron temperature may not be valid under highly non-uniform conditions. In order to describe the non-parabolic momentum flux in terms of the hydrodynamic state variables n,


v and w, the following approximations duced: (k&i) N f(P) 6,j,

are intro-

((1 + 2cf#)-‘) N (1 + 2aw)-‘.



Approximation (12a) is valid as long as the standard deviation of the wave-vector distribution is much larger than the mean value of the distribution, and nearly isotropic, and as long as the wave-vector components are uncorrelated. For the case of a parabolic band, approximation (12a) can be written using Wl+((U,- V,)(U,as: Uj))





$r+(30ivj) + fk, T, 6,~ w 6,.


In general, 3vf # v2 (i.e. v, # vY#a,), but l/2m+v2 is usually small compared with 3/2k,T, in Si devices, and it is often neglected[l&l3]. Consequently, approximation (12a) is analogous to, and no more severe than, the scalar electron temperature approximation commonly invoked for the parabolic model. Approximation (12b) has been shown to be a reasonable approximation even under highly non-uniform conditions[14]l. Using eqns (9), (10) and (12), the momentum flux can be expressed as: n
$(l-(([email protected]))+iu~

=-_nw 1 8,. [

2 3

( 1+2aw >

1 a,,


For a = 0, the parabolic band result is obtained if the convective component of energy is neglected [refer to eqns (5) and (611. For a # 0, eqn (14) predicts a reduction in the momentum flux below its parabolic band value as the average energy increases. As with the average energy, the energy flux depends upon higher order terms which am not known in general. The energy flux is therefore defined as: n (ub) = @w)nwv,


where n(w) is an empirical, average-energy-dependent parameter which is determined from uniform field Monte-Carlo simulations. Under uniform field conditions, eqn (15) is a rigorous definition of the energy flux, but under non-uniform conditions, eqn (15) neglects possible contributions from spatial gradients. As will be shown, however, eqn (15) yields a better description of the energy flux under nonuniform conditions than does the parabolic band definition [refer to eqns (7) and (8)]. The energy flux coefficient n is shown in Fig. 1 as a function of average energy. Our simulations have contirmed that the energy dependence of n at very low energies is inconsequential. Therefore, Cl may be treated as a constant. Note that under uniform conditions at average energies much larger than 1/2m*uz, the parabolic band value of R given by eqn (7) (with Q = 0) is S/3, a value which is nearly 30% larger than the asymptotic value of 1.3 shown in Fig. 1. Thus, eqn

T. J.





where y = 4.2 x 10ez6W. mz K-l and R is a dimensionless parameter[lS]. The conventional parabolic band model (under the assumption of negligible convective energy) is obtained with a = 0, R = 5/3 and R = S/2. The non-parabolic model is obtained with a = 0.5, R = 1.3 and R = 0. 2.2. Monte-Carlo model 0.5 I 0.4

I 0.2




I 0.5



Fig. I. Uniform-field energy flux coefficient vs average electron energy calculated as n(ul)/now using the MonteCarlo model.

(15) insures physical consistency in the homogeneous limit, whereas eqn (7) does not. From eqns (10) and (15), the average momentum is: p = (1 + 2Ruw)m*v. (16) The second term in eqn (16) can be interpreted as accounting for the effect of a non-constant effective mass in the distribution of electrons. As a consequence of eqn (15), our model is similar in form to the one proposed by Stewart and Churchill, who define an m*(w) via a definition analogous to eqn (16). Our approaches differ in that eqn (12) is used here to define an analytic momentum flux, whereas Stewart and Churchill rely on an empirical electron temperature model to describe the momentum flux. From the above results, the momentum and energy conservation equations for the non-parabolic case are, respectively:

5 [(l + 2Gaw)nm*v] =-

+ iv ( nw =)+4nE

(1 + 2Raw)nm *v r&l(w)

-- qnv =



and &v)+V(tiwv+Q) +qnv.E=

-n 2.


Note that the mobility appearing in eqn (17), not the momentum relaxation time, is extracted from uniform field Monte-Carlo simulation data using Shur’s approach. The heat flux Q has been reinstated in eqn (18) in order to allow eqns (17) and (18) to represent both the parabolic and non-parabolic formulations of the hydrodynamic model. The heat flux is modeled by a Fourier law as:


-pRnvw, B

The Monte-Carlo method applied in this analysis takes into account acoustic and f- and g-type intervalley scattering using the scattering parameters of Jacoboni and Reggiani[l6]. An ellipsoidal, nonparabolic band is used with six valley minima located at 0.85 2z/a. Since the valley minima lie close to the Brillouin zone edge, a rejection algorithm has been implemented to account for the decrease in the density of states near the zone edge[l’l]. The simulations are performed using a fixed, positiondependent electric field established along the [Ill] crystallographic axis. 3.



3.1. Model evaluation approach In this section, the parabolic and non-parabolic formulations of the hydrodynamic model are compared with the Monte-Carlo model. For the purpose of model evaluation, simulations are performed using fixed, spatially non-uniform electric fields. A tixed electric field is an ideal approach for model comparison since it guarantees that each transport model experiences the same driving force. Self-consistent simulations involving Poisson’s equation, which are necessary for actual device analysis but not for transport model comparison, introduce variations in the electric field from one transport model to another. These variations may hinder the ability to identify the underlying transport mechanisms which are responsible for the difference in solutions. Steady-state simulation results are presented for two electric field cases which represent low and high gradient regimes. The first electric field profile has a symmetric, Gaussian-like shape (cf. Fig. 2) with a peak field of 200 kV cm-i, and a lateral extent of approx. 0.25 pm. This field will be referred to as the LG (low gradient) field. The second field profile is that of the lateral electric field along the interface of a submicron MOSFET (L,= 0.8 pm). This field, which will be called the MF field, more closely represents the degree of non-uniformity present under realistic device conditions. The steady-state from of eqns (2) and (17-19) are decoupled and solved iteratively, and discretized using standard finite difference techniques. Since the relaxation times do not depend upon the electron concentration, tixed boundary conditions for the electron concentration can be arbitrarily specified as long as they are the same on each boundary (by the symmetry of the problem)[l5]. Note that without the heat flux, eqn (18) becomes a first-order equation in

Hydrodynamic model for Si device simulation

135 -200


E 2 -100

7; 5 r .g ij z w





0.5 x












X (pm) Fig. 2. (a) Average electron energy. (b) Drift velocity for the LG field calculated by Monte-Carlo (O), and by the non-parabolic (-) and parabolic (---) formulations of the hydrodynamic model. Also shown are the results calculated using the parabolic formulation without heat flux (R = 0) (- * -. -). Note that the simulation boundaries arc at x = 0 and x = 1.

average energy (for known velocity and electron concentration). In order to maintain a stable houndary value problem when the heat flux is omitted, eqn (17) is used to explicitly define the energy flux in eqn SSE 35/2-B

(18). The energy gradient term in eqn (17) introduces a second-order energy derivative in eqn (IS) which enhances the stability of the iterative solution method when the heat flux is removed. A more comprehensive


T. J.


description of the numerical technique, including application to the transient problem, will be presented in a future publication. 3.2. Comparison of parabolic and non -parabolic for muiations Average electron energy and drift velocity are shown in Figs 2a and b, respectively, for the LG electric field case. In addition to the parabolic and non-parabolic model results, solutions obtained using the parabolic model without heat flux have been included for comparison. Previous work has indicated that improved performance is obtained with the parabolic model if the heat flux is omitted[l5,18]. Figure 2 shows that excellent agreement with the Monte-Carlo results is achieved with the nonparabolic model in both average energy and drift velocity. The parabolic formulation significantly underestimates the peak average energy, and does not predict the drift velocity profile well. The performance of the parabolic model is improved by removing the heat flux, but only the non-parabolic model accurately captures borh the velocity overshoot and average energy characteristics. Note that without heat flux, the parabolic peak energy is still too low, indicating that the nonparabolic corrections are needed for best agreement with Monte-Carlo since the heat flux will only decrease the average energy. A discrepancy is observed, however, around x = 0.65 pm where the MonteCarlo drift velocity drops more quickly than that predicted by the hydrodynamic models. By including the heat flux, the parabolic model is able to reflect the Monte-Carlo result a little more closely in this region. Figure 3 shows average energy and drift velocity generated by the MF field as calculated by the hydrodynamic and Monte-Carlo models. As with the LG field, the non-parabolic formulation gives the best agreement in drift velocity with the Monte-Carlo result (Fig. 3b). The peak energy, however, is overestimated somewhat (Fig. 3b). Residual drift velocity predicted by the hydrodynamic models is also apparent downstream of the field peak. Note that the parabolic formulation predicts a pronounced increase in drift velocity at the field peak. This overshoot behavior, which has been observed in n +-x-n + diode simulations as well[l4,19-211, is a junction phenomenon associated with the nature of the electron cooling predicted by the hydrodynamic model when the heat flux is included[l5]. 3.3. Energy flux analysis The non-parabolic formulation of the hydrodynamic model accurately models the drift velocity for both field cases except in localized regions downstream of the geld peak. There, the Monte-Carlo solution predicts a sharper decrease in the drift

et al.

velocity. A distinction between the LG and MF field results generated by the non-parabolic model is that the location of the average energy peak coincides more closely with the position of the downstream drift velocity discrepancy for the MF case than it does for the LG case. The average energy peak generated by the LG field is located upstream, away from the region of poor velocity agreement (Fig. 2), and the average energy agreement with the Mont-Carlo result is excellent. For the MF field, the average energy peak is located nearer to where the drift velocity drops (Fig. 3), and the agreement in average energy degrades. This suggests that the differences between the Monte-Carlo and hydrodynamic solutions shown in Figs 2 and 3 are associated with an effect localized to the downstream side of the field peak. Comparison of the drift velocity solutions generated by the parabolic model with and without heat flux indicates that the energy flux plays a role in the observed discrepancies between the hydrodynamic and Monte-Carlo solutions. The accuracy of the energy flux definition provided by eqn (15) can be evaluated by examining the ratio n(ud)/nuw as a function of position for each field case. The actual energy flux n (u 4’) and the convective energy flux nuw are calculated with the Monte-Carlo model. The energy flux ratios for both field cases are shown in Fig. 4 (solid lines) along with the average energy (dashed lines) for comparison. It is clear that an additional flux component not accounted for by eqn (15) exists downstream of the average energy peak in each case. With the LG field, the additional flux appears well after the average energy peaks (Fig. 4a). In the peak energy region, eqn (15) with Q = 1.3 accurately gives the energy flux, and the nonparabolic hydrodynamic solution agrees well with the Monte-Carlo solution. Where eqn (15) does not produce the correct energy flux, the hydrodynamic model overestimates both the drift velocity and the average energy. With the MF field (Fig. 4b), the additional flux is more significant. The peak energy is now affected since the flux deviates from eqn (15) over a sizeable portion of the peak energy region. The flux component is quite large for this electric field profile. Since the average energy decreases most abruptly for this case, the additional energy flux is clearly correlated with the gradient of the average energy. The Fourier heat flux law [eqn (19)] supplies an additional energy flux component where the average energy is decreasing. Thus, the parabolic hydrodynamic model is able to more closely predict the MonteCarlo drift velocity in the region downstream of the field peak where eqn (15) does not hold. The Fourier law, however, also decreases the energy flux where the average energy is increasing. Figure 4 indicates that the heat flux does not appear to follow the linear rule of eqn (19). The open circles in the figure are the uniform-field values of G (cf. Fig 1) associated with the avereage energy at the corresponding locations.


Hydrodynamic model for Si device simulation -250



z -150

p Y ;; z!


*g 4 w











X (wN




1 .o


0 0.3





X (v)

Fig. 3. (a) Average electron energy. (b) DriA velocity for the MF field calculated by Monte-Carlo (0) and by the non-parabolic (-) and parabolic (- - -) formulations of the hydrodynamic model. Also shown are the results calculated using the parabolic formulation without heat flux (R = 0) (- . -. -). Note that the simulation boundaries are at x = 0 and x = 1. The agreement between the flux ratio and R indicates

that eqn (16) holds reasonably well where the average energy gradient is positive (i.e. any gradient-related flux component is small). The Fourier law commonly used overestimates the heat flux and causes the parabolic hydrodynamic model to underestimate the average energy.


A formulation of the hydrodynamic model has been presented for device simulation which accounts for energy band non-parabolicity. The simulation results presented here indicate that the non-parabolic hydrodynamic model is more accurate than the con-

T. J.



et af.






0 0.3







X (rm) 4

’ (b)


.sJ e 5 E


% zi 15

* 1

0, 0.3





X (rm) Fig. 4. The ratio of actual energy flux to convective energy flux, n(uU)/nuw (-) for: (a) the LG field case; and (b) the MF field case. The corresponding average electron energies are also shown for comparison (---). The open circles indicate the uniform-field values of f2 corresponding to the local average energy values. All quantities are calculated with the Monte-Carlo model.

ventional parabolic band formulation of the hydrodynamic model. In contrast to other non-parabolic formulations, the first-order approach we have adopted captures the salient eflbcts of non-parabolicity while closely retaining the form of the parabolic model, and requires only one empirical constant in addition to the relaxation times. Non-parabolic

effects am manifested in the mode1 as reduced diffusion, and as a modified description of the energy flux which, in contrast to the parabolic band model, does not include a Fourier heat flux term. Simulation results have been presented for two electric field profiles: a Oaussian shaped field and a lateral electric field along the interface of a submicron

Hydrodynamic model for Si device simulation MOSFET (15~ = 0.8 pm) (representing low and high gradient conditions, respectively). By comparison with Monte-Carlo and the conventional hydrodynamic model, the non-parabolic model gives improved agreement with Monte-Carlo in both drift velocity and average energy for each case. The simulations reveal, however, that a gradient-related energy flux component exists which is not accounted for in the non-parabolic model. This component can cause the non-parabolic model to slightly overestimate the average energy peak under high gradient conditions, and can cause a localized discrepancy between the hydrodynamic and Monte-Carlo drift velocities. An important observation is that this flux component is significant only when the average energy is &creasing. The conventional treatment of the heat flux in the parabolic model overestimates the heat flux contribution, and is primarily responsible for the poor agreement between the conventional hydrodynamic model and Monte-Carlo. Under sufficiently low gradient conditions, the heat flux contribution is less significant, and the non-parabolic corrections yield excellent agreement with Monte-Carlo results. Acknowledgements-The authors would like to thank Mr Choh-Fei Yeap for helpful discussions concerning the numerical solution of the hydrodynamic model. The authors also gratefully acknowledge support by the Semiconductor Research Corporation, Intel Inc. and Motorola Inc. REFERENCES

1. K. Bletekjaer, IEEE Trans. Electron Devices ED-17, 38 (1970).


2. M. Lundstrom, Transpart FuuaivuentaIs for Device Applications. Addison-Wesley, Reading, MA (1990). 3. M. Fukuma and W. Lui, IEEE Electron. Device Lett. EDI& 214 (1987). 4. E. AaoIT,J. appl. Phys. 64, 2439 (1988). 5. R. Thoma, A. Edmunds, B. Meinerahagen, H. J. Feifer and W. L. Engl, IEEE Tram. Electron Devices ED-38, 1343 (1991). 6. D. L. Woolard, R. J. Trew and M. A. Littlejohn, _ Solid-St. Electron 31, 571 (1988). 7. R. A. Stewart and J. N. Churchill. Solid-St. Electron 33. 819 (1990). 8. M. Shur, -Electronics Lett. 12, 615 (1976). 9. R. Stratton. Phvs. Rev. 126. 2002 (1962). 10. M. Fukuma and R. Uebbing, IED* Tech. Dig., p. 621 (1984). 11. R. Cook and J. Frey, COMPEL 1, 65 (1982). 12. B. Meinerahagen and W. Engl, IEEE Trans. Electron Devices ED-36, 689 (1988). 13. M. Tomiaawa, K. Yokoyama and A. Yoshi, IEEE Tram. CAD 7; 254 (1988): 14. T. J. Bordelon. X.-L. Wana. C. M. Maziar and A. F. Tasch, Proc. ht. Electron--Devices Meeting, p. 353 (1990). 15. T. J. Bordelon, X.-L. Wang, C. M. Maziar and A. F. Tasch, Solid-St. Electron 34, 617 (1991). 16. C. Jacoboni and L. Reggiani, Rec. Mod. Phys 55, 645 (1983). 17. X.-L. Wang, T. J. Bordelon, C. M. Maziar and A. F. Tasch, Proc. ISDRS-9f (1991). 18. H.-H. Ou and T.-W. Tang, IEEE Trans. Electron Deuices ED-34, 1533 (1987).19. C. Gardner. J. Jerome and D. Rose. IEEE Tram. Electron D&ices ED-8, 501 (1989). 20. E. Fatemi, J. Jerome and S. Osher, IEEE Trans. CAD. 10, 232 (1991). 21. M. Rudan, F. Ckleh and J. White, COMPEL 6, 151 (1987).