The Response of the Hydrodynamic Model to Heat Conduction , Mobility , and Relaxation Expressions

We study simulations of the n +-n-n + diode, by means of higher moment models, derived from the Boltzmann equation. We emply such realistic assumptions as energy dependent mobility functions, with doping dependent low field mobility. It is known that a critical role is played in the hydrodynamic model by the heat conduction term. When the standard choiCe of the Wiedemann-Franz law is made for the conductivity, and constant low field mobility values are used, spurious overshoot is observed. Agreement with Monte-Carlo simulation in this regime has in the past been achieved by empirical modification of this law. In this paper, we consider the effect of representing the heat flux by the sum of two terms. It is found that the effect is negligible with respect to overshoot in comparison to that achieved by employing a doping dependent low field mobility. We also compare the hydrodynamic model to recently introduced energy transport models. Finally, in low temperature regimes, we study the dependence of shock formation on the momentum relaxation time representations and on the heat conduction term. The algorithms employed for both models are the essentially nonoscillatory (ENO) shock capturing algorithms.


INTRODUCTION
n previous work, we have demonstrated the ro- bustness of an algorithm (ENO" Essentially Non-Oscillatory) designed for conservation systems ad- mitting steep fronts and possible shock like behav- ior.The algorithm captures steady states by explicit time stepping and identification of characteristic directions.It is an apt choice for the simulation of the hydrodynamic model for semiconductors over a wide range of parameters.In [7] and [8], n/-n-n / diodes in one dimension and MESFETS in two dimensions were simulated, for both the hydrodynamic model and an energy transport model (ET), introduced in [4].Both of these models are charac- terized as self-consistent higher order moment models, expressing various conservation laws.The most significant difference between the hydrodynamic and ET models is that the former is based upon macroscopic and the latter upon microscopic relaxation time approximations.Also, the hydrodynamic model contains hyperbolic modes not found in the ET model.
An issue of continuing debate in the literature is the role of heat conduction in the hydrodynamic model.The spurious velocity overshoot at the chan- nel/drain junction observed in simulations has been addressed by various authors; Gnudi, Odeh, and Rudan [6] appear to have been the first to present the problem clearly, and to indicate the appropriate empirical choice of the thermal conductivity con- stant required to establish consistency of the data 132 JOSEPH W. JEROME and CHI-WANG SHU with respect to Monte-Carlo simulations.In [13], the traditional choice of the heat flux representation, in terms of the Wiedemann-Franz law, was critically examined.A first principles' study was employed to confirm the appropriateness of adding a transport type term; this had also been discussed by Lee and Tang [9].The authors of [13] convincingly demon- strated that the new flux gave far more reliable agreement with the rigorous flux as extracted from Monte-Carlo simulations.
The original goal of this work was to detail the effect of this use of the heat flux in the hydrodynamic model, with particular attention paid to alter- ation of the spurious velocity overshoot at the right junction.We felt that an extremely careful study was in order, and carried out simulations which included the effects of doping dependent low field mobility as a multiplier of a temperature dependent mobility.It was a surprising discovery that it is this low field mobility quantity, not the general structure of the flux, which most influences the overshoot.We were also interested in determining whether the flux con- tributed in any way to shock prediction.Although this seems unlikely from first principles, since it is the convective term in the momentum subsystem which, together with supersonic/subsonic transi- tions, determine shocks, we felt that such transitions might be inhibited if heat flux were correctly cali- brated.Our simulations show that this is possibly the case.In the process, we also discovered the decisive part played by the momentum relaxation term in possible shock formation.Certain relaxation expressions do, in fact, inhibit the supersonic/sub- sonic transitions.Moreover, these are precisely the expressions employed in Monte-Carlo simulations.Our shock studies were carried out at 77 K, whereas the other simulations assume an ambient temperature of 300 K.All simulations are for Silicon.We coordinated our case studies at 77 K with the work of Gardner [5].
In setting the goals of the paper, we felt it to be instructive to provide comparisons between models.
Therefore we have included simulations of the en- ergy transport model mentioned earlier.These are included at 300 K only, due to the lack of parameter data at 77 K. Summaries of the models and algorithm are presented in the second section, whereas the third section contains the graphs of extensive simulations.We confirm, for example, the finding of [13] that the Wiedemann-Franz law overstates the (signed) flux in significant portions of the device.This is perhaps the principal reason for making use of the modified heat flux.Concluding remarks are given in the last section.

The Hydrodynamic Model and Parameters
The equations as presented here are discussed in references [3] and [10].They are derived as the first three moments of the Boltzmann equation,  Here, f f(x, u, t) is the numerical distribution function of the electron species, u is the species' group velocity vector, E E(x,t) is the electric field, e is the electron charge modulus, m is the effective electron mass, and C is the time rate of change of f due to collisions.The equations are given by: tgn TIvI + mne, + V. v Tlvl + mne, -env.E-V.(vP)-V.Q+C w. (2.4) The Poisson equation for the electric potential must be adjoined; each species contributes a corresponding moment subsystem, with appropriately signed charge.We begin with standard definitions and as- sumptions.The concentration is given by n .'=ffdu; the average velocity by v (1/n)fufdu; the mo- mentum by p := mnv; the random velocity by c u v; the pressure tensor by Pij '= mfcic,fdu; and the internal energy density by et := (1/2n)flcl2fdu.
This function represents energy/unit mass/unit concentration.The heat flux Q is given by Qi (m/2)fcilcl2fdu.Finally, for reference in subse- quent subsections, the electron current density is given by J -env, and the energy flux is given by S fu{(m/2)lul2}fdu.The assumptions on f are now stated.The function f is assumed to decrease sufficiently rapidly at infinity: In addition to these transport equations, we have Poisson's equation for the electric field, where n d doping and e dielectric: V'(oV) .,eini n d. (2.6) Here, we have used the convention that there are different species, each of concentration n and charge e i.The entire system consists of equations (2.2), (2.3), (2.4), repeated according to species, and (2.5), (2.6).By moment closure is meant the selection of com- patible relations among the variables, 4, n, v, P, et, and Q, where attention is confined now to a single species, so that the number of equations is equal in number to the remaining primitive variables se- lected.We begin by introducing a new tensor vari- able T, the carrier temperature, defined by Pij nkTi/, where k is Boltzmann's constant, and a scalar vari- able W, the total carrier "energy.Reduction to a set of basic variables, n, v, W, and b, or a set equiva- lent to these, can be implemented by the following assumptions: 1.The pressure tensor is isotropic, with diagonal entries P and off-diagonal entries zero, for a suitable scalar function, P. P is related to e 3 via mnet Ps.

It follows from the previous assumption that
temperature may be represented by a scalar quantity T, and that the internal energy is represented in terms of T by me -kT.
3. The total energy density (per unit concentra- tion) w is given by combining internal energy and parabolic energy bands with m assumed constant.
2 w me I + -mlvl.. and the total energy (per unit volume) W is the product, W nw. The heat flux is obtained by a differential expression involving the temperature (cf.(2.15), (2.17) to follow).
It is possible to rewrite the system (2.2, 2.3, 2.4) with the closure assumptions incorporated.We have the following. (vnkT)-V.Q+C w. (2.9) The form of Q is pivotal and is discussed below.The final step deals with the replacement of the collision moments.For a one carrier system, we define C n 0, and the momentum and energy re- laxation times, "rp and z, respectively, in terms of averaged collision moments as follows.
2. The energy relaxation time r w is given via m 2 lul2Cdu:=C w.
Here, W0 denotes the rest energy, 3 kTo, where T O is the lattice temperature.The forms for the relaxation times selected at 300 K while at 77 K the form of z w is given by the modified Baccarani-Wordeman model, " 2 r P ( l}kT)2 (212) rw=-:-1+ mv Here, rpo (m/e)t,o and two 3kT o t,o/2ev2s.Typ- ically, the low field mobility, /z0, is taken to be constant, but in this paper we also investigate a doping dependent version given in units of Ixm2/V/ps by 14.2273 /Zo(X) 0.0088. 1.0 + nd (2.13 t 1.0 + 143200 Here, n d is expressed in units per /zm 3.This for- mula is used only for the 300 K case.We note that it does not directly take into account the fact that the influence of ionized impurity scattering is known to decrease for high carrier energies, and perhaps overestimates such influence at the right contact. Note, however, that the momentum relaxation time is proportional to the mobility, which means that the mobility functions are energy dependent in their complete representation.Moreover, in the simula- tions performed at 77 K, we also compare with replacing rp by a version motivated by Monte-Carlo simulations; this has the effect of introducing a divisor, 6,, (1 + (nd/nrel,)), into (2.10):rp rvo 6. (2.14) Here, the values a .659and nref= 1440//,m 3 are employed.
The traditional form of Q (denoted standard Q in the simulations) has been, Q -rVT. (2.15) Here, K is the thermal conductivity governed by the Wiedemann-Franz law (cf.[2]).This may be described by, TnTe oo (2.16)   where r is typically taken to be .6,giving a value of 1.5 for the dimensionless coefficient.The new form of Q, suggested in [13], is 5(1 r) kT o J. (2.17 The value of r .72 is recommended by the authors of [13].Thus, the kinetic heat flux has both convec- tive and diffusive components.In order to obtain comparative uniformity, we select r .72 in both (2.15) and (2.17)when T O 300 K.For T O 77 K, we select r .02, to conform to the selections made for the thermal conductivity in [5].Strictly speaking, the authors of [13] not only found that Fourier's law Was violated, but also found that the mobility was a nonlocal function of the electron temperature.It is therefore in some sense inconsistent to use the heat flux model of [13] with a temperature dependent mobility.Our model should therefore be viewed as a first step toward understanding these effects.

The Energy Transport Model
In this section, we shall employ a microscopic as- sumption upon the momentum relaxation time, viz., we shall assume that the collision term C in (2.1) rp where fl is the odd part of f.Note that this con- trasts with tee macroscopic assumption on rp, em- ployed in the hydrodynamic model.In this case, we may obtain an expression for the energy flux S: where /x w and D w are tensor expressions for mobil- ity and diffusion, defined in terms of moments, and w represents average energy per unit concentration.The details are furnished in [4].The following con- stitutive relation holds for w: w }kT(1 + -flkT).
(2.20) Equation (2.20) allows for nonparabolic energy bands as well as non-Maxwellian distributions.Alto- gether, the model may be written in terms of the Poisson equation, (2.6), in conjunction with the sys- tem, schemes, consider the scalar one dimensional problem, and a conservative approximation of the spatial operator given by V-J= 0, (2.21) V S J. E nC. (2.22) Here, J is the drift-diffusion current, J -elznV + eV(nD), (2.23) and S has been described previously, in (2.19).In the expressions for J and S, the assumptions made for the model lead to scalar representations for the mobility and diffusion coefficients.For example, the choice made in [4] leads to ix po Tof T /x -1 kr (2.24) tzkr. (2.25) The diffusion coefficients are defined by Einstein's relations.The collision term in (2.22) is a specified quadratic function of T. One significant advantage of the microscopic relaxation time assumption is that certain key parameters may be fitted via Monte-Carlo preprocessing, ensuring reliability of their values.

Algorithm
We shall only briefly describe the algorithm used in this paper, namely the ENO scheme developed in [11] and [12].The ENO scheme is designed for a system of hyperbolic conservation laws of the form, ut + f(U)x g(u, x, t), (2.26)   where u (Ul,... Um )T, and the hyperbolicity con- dition, is diagonalizable with real eigenvalues oU holds.An initial condition is adjoined to (2.26).
For systems of conservation laws, local field by field decomposition is used, to resolve waves in different characteristic directions.For this purpose, analytical expressions are needed for the eigenvalues and eigenvectors of the Jacobian matrix.This reduces the determination of the scheme to the case of a single conservation law.Thus, to describe the L(u)g Ax ('+1/2 -'-1/2)" (2.27) Here, the numerical flux f is assumed consistent: +1/2 f (U_l,...,U+k); f (u,...,u)=f(u). (2.28) The conservative scheme (2.27), which character- izes the f divided difference as an approximation to f(u),, suggests that f can be identified with an appropriate function h satisfying f(u(x)) Jx-(ax/2) h( ) ds .(2.29)If H is any primitive of h, then h can be computed from H'. H itself can be approximated by polyno- mial interpolations using Newton's divided differ- ence method, beginning with differences of order one, since the constant term is arbitrary.The neces- sary divided differences of H, of a given order, are expressed as constant multiples of those of f of order one lower.The main ingredientof the ENO method is the adaptive choice of stencil: it begins with a starting point to the left or right of the current "cell"by means of upwinding, as deter- mined by the sign of the derivative of a selected flux (or the eigenvalue of the Jacobian in the system case); as the order of the divided differences is increased, the divided differences themselves deter- mine the stencil: the "smaller" divided difference is chosen from two possible choices at each stage, ensuring a smoothest fit.In all the calculations for next section, we use the third order ENO scheme with 300 uniformly spaced grid points.We have verified the resolution of the scheme by overlapping the results obtained by using 300 grid points with those obtained by using 600 grid points (not shown).
Convergence towards the steady state is obtained in a time asymptotic fashion.A continuation with vbias is used, with the steady state solution from a lower vbias used as the initial condition for a higher one.

SIMULATION RESULTS
We first simulate a standard n /-n-n / silicon diode with a length of 0.6/zm.The high and low dopings FIGURE 1 The doping profile nd of the 300 K case are n d 5 X 105/zm -3 and n d 2 X 103/m -3, re- spectively, with slightly smoothed out junctions at x 0.1 and x 0.5 (Fig. 1).We point out that the purpose of smoothing out the junction is to have it completely resolved by at least 5-8 grid points.
The "shock capturing" ENO algorithm we use is capable of handling abrupt junctions without any smoothing, although a slight numerical kink is some- times observed at the junction.We will not show results with abrupt junctions in this paper.The lattice temperature T O is taken as T O 300 K.The voltage bias applied is vbias 1.5V.We first compute with the standard Q given by (2.15) and a constant /x 0 -0.14.The result, as is well known, shows a large spurious velocity over- shoot at the right junction (Fig. 2).We plot this "standard Q, constant /x0" case (SQCM) as a solid line background for all the pictures from Fig. 2 to Fig. 5, in which we draw the velocity v in Fig. 2, the temperature T in Fig. 3, the total energy W in Fig. 4, and the heat flux Q in Fig. 5.The concentration n and the electric fieldnth are not very sensitive to the changes we make, hence their graphs are not shown.Five cases are compared with the back- ground case of "standard Q, constant/z0" (SQCM).These are: standard Q with a variable /z0, which depends on the doping (2.13) (SQVM), in Fig. 2(a) through 5(a); the new Q as defined by (2.17) with a constant /z 0 (NQCM), in Fig. 2(b) through 5(b); the new Q as defined by (2.17) with the variable /x0, defined by (2.13) (NQVM), in Fig. 2(c) through 5(c); the ET model (2.21, 2.22)with a constant /x 0 (ETCM), in Fig. 2(d) through 5(d); and the ET model with variable /x 0 (2.13) (ETVM), in Fig. 2(e) through 5(e).
From Figures 2 through 5 we can observe: (1) The variable/x 0 has the most significant effect upon reducing the spurious velocity overshoot at the right junction (Fig. 2(a)).The new Q even increases the spurious overshoot (Fig. 2(b)), although combined with a variable/x0, it does show a decrease of the spurious over- shoot together with a decrease in velocity in the entire middle low doping region (Fig. 2(c)).
The velocity obtained by the ET model, on the other hand, does not seem to be too sensitive to the choice of /x 0 (Fig. 2(d) and 2(e)); (2) The differences in the concentration n (not shown), the temperature T (in Fig. 3), and the electric fieldb' (not shown) are relatively small; (3) The total energy shows a significant change near the right junction for all cases (Fig. 4); (4) The heat flux Q is significantly reduced at the right junction in all cases (Fig. 5), which is consistent with the discussion of [13].Next we compute a shorter silicon diode with a length of 0.45/xm.The high and low dopings are n d 106m -3 and n d 103/zm -3, respectively, with slightly smoothed out junctions at x 0.1 and x 0.35.The lattice temperature T O is taken as T O 77 K.The voltage bias applied is vbias 1V.This is the case where shock-like behavior is found in [5].
We again start with computing with the standard Q given by (2.15) and momentum relaxation time Zp given by (2.10).The result is the same as those shown by Gardner in [5] and we can see shocks especially clearly in the velocity (Fig. 6).We plot this "standard Q, standard Zp" (SQST) case as a solid line background for all the pictures from Fig. 6 to Fig. 10, in which we draw the velocity v in Fig. 6, the temperature T in Fig. 7, the total energy W in Fig. 8, the electric field-'-th' in Fig. 9, and the heat flux Q in Fig. 10.Three cases are compared with the background case of standard Q and standard Zp (SQST): new Q (2.17)with standard Zp (NQST), in the plus symbol; standard Q with the nd-dependent Zp as defined by (2.14) FIGURE 7 The temperature T of the 77 K case--with the standard Q, standard 'p result (SQST) as the solid line back- ground.Plus signs: new Q, standard ', (NQST); circle signs: standard Q, hal-dependent -p (SQNT);-star signs: new Q, r d- dependent -p (NQNT) FIGURE 8 The total energy W of the 77 K case--with the standard Q, standard -, result (SQST) as the solid line back- ground.Plus signs: new Q, standard -, (NQST); circle signs: standard Q, hal-dependent -, (SQNT);-star signs: new Q, n d- dependent -p (NQNT)  FIGURE 9 The electric field--b' of the 77 K case--with the standard Q, standard -, result (SQST) as the solid line back- ground.Plus signs: new Q, standard -p (NQST); circle signs: standard Q, nd-dependent -p (SQNT); star signs: new Q, n d dependent -p (NQNT) FIGURE 10 The heat flux Q of the 77 K case--with the standard Q, standard -p result (SQST) as the solid line back- ground.Plus signs: new Q, standard -, (NQST); circle signs: standard Q, nd-dependent 't, (SQNT);-star signs: new Q, n d- dependent -p (NQNT) and the new Q with the nd-dependent p as defined by (2.14) (NQNT), in the star symbol.
From Figures 6 through 10 we can observe that the changes are now much more significant than in the 300 K case.In particular, the shock seems to disappear in all the three cases.

CONCLUSIONS
We have demonstrated the robustness of a shock capturing algorithm (ENO) to simulate the hydrody- namic model, with convection, at 300 K and 77 K.The role of a new heat conduction expression, con- taining both convective and diffusive components, is analyzed for its effects upon the model.It is deter- mined that the heat flux is accordingly more accu- rately represented throughout the device than with the Fourier law.However, this refined heat conduc- tion expression has negligible impact upon the so- called spurious velocity overshoot at the drain junc- tion which is most impacted at 300 K by the doping dependent low field mobility.The energy transport model, on the other hand, is more robust to such changes in low field mobility, and does not en- counter the problem of spurious overshoot.At 77 K, the heat conduction and the momentum relaxation appear to influence decisively the formation of shocks in the hydrodynamic model.We did not employ the energy transport model in this case due to lack of parametric data.
oW ot +V.(vW) -env.E-V are the Baccarani-Wordeman models [1], device, in /,m The moment equa- tions are expressed in terms of certain dependent variables, where n is the electron concentration, o is the average velocity, p is the momentum density, P is the symmetric pressure tensor, Q is the heat flux, et is the internal energy, and C., Up, and Cw (SQNT)in the circle symbol; The velocity v of the 77 K case--with the standard Q, standard % result (SQST) as the solid line background.Plus signs: new Q, standard '. (NQST); circle signs: standard Q,