Analytical and Computational Advances for Hydrodynamic Models of Classical and Quantum Charge Transport

In recent years, substantial advances have been made in understanding hydrodynamic models, both from the standpoint of analytical infrastructure, as well as the parameters which play a decisive effect in the behavior of such models. Both classical and quantum hydrodynamic models have been studied in depth. In this survey paper, we describe several results of this type. We include, for example, well-posedness for both classical and quantum reduced models, and the relaxation drift-diffusion limit as examples of analytical results. As examples of computational results, we include some discussion of effective algorithms, but most importantly, some information gleaned from extensive simulation. In particular, we present our findings of the prominent role played by the mobilities in the classical models, and the role of hysteresis in the quantum models. All models are self-consistent. Included is discussion of recent analytical results on the use of Maxwell’s equations. Benchmark devices are utilized: the MESFET transistor and the n+/n/n+ diode for classical transport, and the resonant tunneling diode for quantum transport. Some comparison with the linear Boltzmann transport equation is included. Acknowledgements: The author is supported by the National Science Foundation under grant DMS-9704458. ∗Department of Mathematics, Northwestern University, Evanston, Il 60208 1

The equations as presented here are discussed in references [7, 27 and  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, rn is the effective electron mass, and C is the time rate of change of f due to collisions.The moment *e-mail: jwj @math.nwu.eduequations are expressed in terms of certain dependent variables, where n is the electron concentration, v is the average velocity, p is the momentum density, P is the symmetric pressure tensor, Q is the heat flux, ei is the internal energy, and Cn, Cp and Cw represent moments of C, taken with respect to the functions, ho(u) -= 1, ha(u) mu, h2(u) (m/2)lul .T he equations are: Ot On cot --+ V. (nv) Cn, (1.2) --+ v(V.p) + (p.V)v+ V.P -enE+ Cp, (1.3) 0 (mn ) ({mnlv[2+mnei)) 35 --Ivl= / mnei + v -+ .(vP) + .Q -env.E + Cw. (1.4) Poisson's electrostatic equation for the electric potential must be adjoined; each species contri- butes a corresponding moment subsystem, with appropriately signed charge.The concentration is given by n := ffdu; the average velocity by v := (1 In) f uf du; the momentum by p := tony.
Finally, for reference in subsequent subsections, the electron current density is given by J :=-env.In addition to these transport equations, we have Poisson's equation for the electric field, where nd := doping and e := dielectric: E -Vb, (1.5) V" (eV)) Z eini na. (1.6)Here, we have used the convention that there are different species, each of concentration ni and charge ei.The entire system consists of Eqs.(1.2)-(1.4),repeated according to species, with possible coupling terms, and (1.5), (1.6).
By moment closure is meant the selection of compatible relations among the variables, b,n, v,P, ei and Q.We begin by introducing a new variable T, the carrier temperature, defined by Po nkT6o., where k is Boltzmann's constant, and a scalar variable w, the total carrier energy.Reduction to a set of basic variables, n, v, w and qS, or a set equivalent to these, can be implemented.
The internal energy is expressed by mel-(3/2)kT.We then have: 1.The total energy density (per unit volume) w is given by combining internal energy and para- bolic energy bands with rn assumed constant: The heat flux is obtained by a differential ex- pression involving the temperature (cf.(1.9) to follow).
The final step deals with the replacement of the collision moments.For a one carrier system, we define Cn=0, and the momentum and energy relaxation times, rp and rw, respectively, in terms of averaged collision moments as follows.We have: 1.The momentum relaxation time rp is given via P-P-:= -m f uC du "= rp J 2. The energy relaxation time rw is given via w-wo.m f 2 c wlul du'-fw, Here, Wo denotes the rest energy, (3/2)kTo, where To is the lattice temperature.
The forms for the relaxation times selected at 300 K are the Baccarani-Wordeman models [5], T rw rWO T +.To +-rp. (1.8) Here, rpo=(m/e)#o and rwO=(3kTo#o/2ev2), where vs is the saturation velocity.Typically, the low field mobility #o is taken to be constant, but other choices are possible, such as doping dependent mobilities.Units are given as tmZ/V/ps.
The traditional form of Q has been (see [31] for more general expressions), Dp (kT/e)#p. (1.15) A detailed existence and approximation theory for this model was presented in [22].
Q -nVT. (1.9) Here, n is the thermal conductivity governed by the Wiedemann-Franz law.This may be de- scribed by, 5r k21zo (T) (1.10 where r is typically taken to be .6,giving a value of 1.5 for the dimensionless coefficient.The wellposedness of the reduced hydrodynamic model, for two carriers, was demonstrated in [11].

The Drift-Diffusion Model
The drift-diffusion model may be obtained by taking zeroth order moments of the BTE and adjoining the Poisson equation.Thus, one obtains the system for N carriers with recombination Ri, current density Ji, (signed) charge ei, i--1,..., N: eiOni Ot -k-V.Ji -eiRi, (1.11The model used in this paper was derived by Gardner in [20].In this section, we shall review the basic characteristics of the model as described in [14].An existence theorem for the reduced model was obtained in [33].The model is also discussed in [22]. The QHD model has exactly the same structure as the classical hydrodynamic model (electrogasdynamics), where we now permit a non-isotropic stress tensor: OV (w-(3/2)nTo) -nvi OX 7"w (1.18) together with (1.5) and (1.6).There still remains the issue of determining the constitutive current relations.Classical drift-diffusion theory gives, for N 2, nl n and nz =p, Jn -e#nnV + eDnVn, (1.12) Jp -e#ppV eDpVp. (1.13) The electronic charge modulus e is positive here, and n and p denote the electron and hole densities, respectively.The use of the Einstein relations linking the mobilities, #n,#p, and the diffusion coefficients, On, Op, is common.These relations are specified by Dn (kT/e)IZn, (1.14) in conjunction with Poisson's equation, (1.6), where V -e is the potential energy, and tem- perature here is expressed in energy units (k is set equal to 1).Spatial indices i,j equal 1, 2, 3, and repeated indices are summed over.Quantum mechanical effects appear in the stress tensor and the energy density.Gardner derived the stress tensor and the energy density based upon the O(h 2) momentum-shifted thermal equilibrium Wigner distribution function: -h2n 0 2 Po.=-nT6ij-lZmOxiOxjlog(n)+O(-h4), (1.19)   3 hZn 72 log(n) / O(h4) w -nT +-mnu 2 24--- (1.20)In one dimension, the QHD model requires eight boundary conditions.Well-posed boundary conditions for the resonant tunneling diode are n=na, On/Ox=O, and OT/Ox=O at the left and right diode boundaries x and xn, with a bias A V across the device: V(XL) T log(n/n.) and V(xn) T log(n/n.)+eAV,where n. is the in- trinsic electron concentration.

Euler-Maxwell Systems
We shall also examine a form of the Euler- Maxwell system for semiconductor transport.In particular, we shall consider the 'isentropic' mass-momentum transport system, coupled to Maxwell's equations for the electric and magnetic fields, instead of Poisson's equation for the electric field only.This Euler-Maxwell system in the isentropic case assumes the following form [3, 4]: H O, B #H, J -env, x E R3, t>0, (1.21) where p=p(n)=n'/7 is the pressure of the flow, expressed in normalized units, "y > is the adiabatic exponent, H E R 3 is the magnetic field, B E R 3 is the magnetic induction, -e(E + v B) is the Lorentz force, and/ is the permeability of the medium.

Vs
We first scale the variables (n, v, T, 4), and then show that the limit functions satisfy the drift-dif- fusion system as "tp, "tw 0 in the manner specified.
We assume that initial data (no(x),vo(x), To(x)-) are independent of "tp and "tw and are "small" as measured in appropriate norms.
The limit function u satisfies the drift-diffusion equation, + (,(, s) n())dcl, x O, x & the sense of distributions.The double integral may be identified with bx, where b is the electro- static potential in the Poisson equation.Note that physical constants have been homogenized to unity in the limiting equation.

RESULTS ON WELL-POSEDNESS FOR THE CLASSICAL MODEL
represents a perturbation of such a system, the homogeneous system is the starting point in the construction of the approximate solutions.We find it convenient to use the notation p-ran, and formulate the system in terms of p, the mass density.Consider the homogeneous system: ut+f(u)x=O, 0<x< 1, ( where u (p, p)-V and f (u) p, (p2/p) + p(p))with P(p)= P'/7, 7 > 1.
The discontinuity in the weak solution of (3.1) satisfies the Rankine-Hugoniot condition: A detailed study of the well-posedness of the initial/boundary-value problem for the reduced hydrodynamic model was carried out in [11].
The reduction refers to the assumption of a pressure-density relationship of adiabatic type, and the subsequent elimination of the energy equation from the system.The mathematical theory of (perturbed) conservation laws is not capable of handling the full system in terms of demonstrating the existence of physical solutions to the initial/ boundary-value problem.We briefly sketch a few rudiments of the theory in the following sub- sections.The interested reader may consult the references for elaboration.The basic idea is to construct approximate solutions, via so-called Riemann problems.The extraction of a conver- gent subsequential limit involves a relatively new idea of compensated compactness.The analysis of [11] actually includes the two-carrier case, as well as more complicated equations arising from multi- dimensional symmetry reduction to one dimen- sion.We now describe the basic idea of Riemann problems in the homogeneous case.

Riemann Problems as Building Blocks
In this section, we review some basic facts about the Riemann solutions for homogeneous systems.
Although the reduced hydrodynamic model where is the propagation speed of the disconti- nuity, and u0 and u are the corresponding left state and right state, respectively.The shock with speed r 0 is called the standing shock.
For the Riemann problem with data (3.3) and the Riemann initial-boundary problem of (3.1) with data: ul=o u+, Plx=o O, (3.4) we have the following facts regarding the solutions.LEMMA 3.1 There exists a piecewise smooth entropy solution u(x, t)for each of the problems (3.3) and (3.4), respectively, satisfying an invariant region condition for some region -.That is, if the Riemann data lie in ', then the Riemann solutions u(x, t) and (1 /(b a) fb u(x, t)dx ].For the Riemann initial-boundary problem of (3.1) with data: ult:o u_, Plx=l O, (3.5) we have the similar results to those for (3.4).
The entropy solutions are the physically relevant solutions.
The method by which the Riemann solutions are utilized in the reduced hydrodynamic model is by a variant of Godunov's method, well-known in gas dynamics, where typically there are no forcing terms.Piecewise constant starting values, an approximation to initial data, are used to con- struct solutions locally in time by the theory of Riemann solvers; the perturbation terms are evaluated explicitly via a fractional step, and the procedure is repeated to advance in time.The Courant-Friedrich's condition on the time step is essential for the construction.One thereby obtains a family of approximate solutions, indexed by spatial grid length.The next subsection sketches the underlying theory which permits the extraction of a limit.

H -1 Compactness of Entropy Measures
In the theory of conservation laws, the notion of entropy pairs, or weak entropy pairs, plays a decisive role.They are (usually) convex functions which may be thought of as facilitating the estimates.We have the following.A pair of map- pings (r/,q):RZ--R 2 is called an entropy-en- tropy flux pair if Vq Vr/Vf If (p, v) rl(p, pv) satisfies (0, v)= 0, for any fixed v=(p/p), then r/ is called a weak entropy.For example, the mechanical energy-energy flux pair p2 r/, -)--+ '7('7 lp2 pT-1 q* P ---+ "7 -1 (3.6) is a strictly convex weak entropy pair for (3.1).
In order to extract convergent sequences from approximate solutions, it is necessary to have a compactness criterion in function space.We need the following basic lemma (cf. [10, 16,32]).For readers not familiar with the notation of the following lemma, we are employing negatively indexed (dual) Sobolev spaces of distributions.
LEMMA 3.2 Let C R N be a bounded domain.
The approximate solutions satisfy this lemma.This is expressed in the following proposition.PROPOSITION 3.1 If {U h }, 1,2, are the approxi- mate solutions, then the measure sequence l(uhi + q( uhi )x is a compact subset ofH[o (f) for all weak entropy pairs 07, q), where f is any bounded and open set in the space-time domain Hr.
The compensated compactness framework is designed to handle the case when the approximate solutions satisfy an invariant region principle, and when the entropy pairs are compact in a suit- able topology.We have the following framework (see [10]): LEMMA 3.3 Assume that the approximate solu- tions u h (ph,ph) satisfy (1) There is a constant C > 0 such that O<_p h (X, t) <_ C, [ph(x, t)/ph(x, t)[ _< C.
(2) The measure 'l(uh)tnt-q(uh)x is compact in Ho (f), for all weak entropy pairs 07, q), where Then, for < '7 < 5/3, there ex&ts a convergent subsequence (still labeled uh) such that uh(x, t)- u(x, t)= (p(x, t), p(x, t)), a.e.This is the final component required to demon- strate that u is a weak solution of the reduced hydrodynamic system.Additional references, such as [26], may be found in [11].

EXISTENCE FOR THE QHD MODEL
The quantum hydrodynamic (QHD) model is a moment model, derived from the Wigner equation.It may be viewed as a quantum corrected version of the classical hydrodynamic equations, with the stress tensor and the energy density corrected by O(h 2) perturbations.Some applications, such as the resonant tunnel diode, also involve quantum well potentials.Ancona, Iafrate and Tiersten [2,1] derived the expression for the stress tensor, and Grubin and Kreskovsky formulated a one dimen- sional version of the model [21].The model which we study is of physical importance; in fact, a sim- plification of our model has been characterized as a pure state, single carrier transport model in [21].Note that a pure state model would have the factor (-1/4) in (4.5) whereas the mixed state model has the factor (-1/12).
Our approach is completely novel to this appli- cation area, in that we reduce the system to an integro-differential equation, with a set of bound- ary conditions, including a nonstandard second- order boundary condition, which is equivalent to specifying the quantum potential at the (current) inflow boundary.
Quantum mechanics is represented by the quan- tum potential [1]: The device domain is the x-interval, I= (0, 1).
In this section we describe the steady-state case rt (rtv)t 0.Then, after the introduction of the current density j=nv, the system (4.1)-(4.3)reduces to j(x) const., (4.6) Cxx -e(na n).

Formulation and Summary of Result
We shall present the equations for the simplified QHD model as developed in [21] and [20].

The Integral Equation
In order to transform the above equations, we first use a Green's function to solve Poisson's equation (4.8), and then we reduce the system (4.7)-(4.8) to an integro-differential equation.We shall see that the existence of a smooth solution of the original system is equivalent to that of a smooth solution of the integro-differential equation.
That is, the third boundary condition of (4.9) holds for w.Hence, the existence of a smooth solution of (4.6)-(4.10) is equivalent to that of a smooth solution of (4.16), (4.17), provided n does not vanish on L The main result was proven in [33].With minor changes, it reads as follows (due to the correction of some sign errors in [33]): Then, for each j > 0, there exists a classical solu- tion (n, 4) of (4.7)-(4.10),satisfying the properties that n E C3(I) (for compact subsets of I) and (1/K) < < K for some K > O.
5. EULER-MAXWELL WELL-POSEDNESS our approximate solution Uh(t) in the space We shall be brief in this section since we have already described the general approach to Euler- Poisson systems in 3. The approach of [12] adopts this method, but extends the ideas to the Euler- Maxwell system in one spatial variable.The critical boundary conditions, as before, are with the momentum: Plx=0 Plx=l 0. (5.1) For general large initial data in L , the solutions of (1.21) will develop singularities or shocks in finite time.Therefore, there are only global weak solutions, including shock waves, for general large initial data.The global solution of (1.21) is con- structed under the assumption of a linearized Lorentz force.In this sense, the well-posedness result is partial.If the initial-boundary conditions are bounded, the global solution will be bounded.
The global approximate solutions constructed by the Godunov method with the fractional step procedure, are then shown to be convergent to the global weak solution.

NUMERICAL ALGORITHMS
The simulations contained in [23,24 and 25], as well as the earlier simulations in [18,17], were based upon ENO schemes, or upon careful calibration between ENO and upwinding.Such schemes are actually discrete Godunov schemes.Alternative algorithms, particularly those based upon upwinding finite element methods, are per- haps less known, and we mention these briefly now.These were particularly successful in the studies carried out for the QHD model in [14].
The method is described as follows.First, we triangulate our domain $2 with triangulations "Th made solely of rectangles R such that the intersec- tion of two distinct rectangles of the triangulation Th is either an edge, a vertex, or void.Then, for each E (0, tf], we take each of the components of Vh= {pEL(f)'p[lislinear, VRTh}.( We define each of the components of u0h to be the L2-projection of the corresponding component of u0 into Vh and discretize the Eq. (3.1) in space by using the Discontinuous Galerkin (DG) method.Since the functions of the space Vh are discon- tinuous, the mass matrix of the DG method is block-diagonal.Thus, the resulting discrete equa- tions can be rewritten as the following ODE initial value problem: duh dt Lh(uh, g) + Rh(uh), E (0, T], ( Uh(t 0)--U0h, ( where Lh is the approximation of-V.F.The exact solution of the above initial value problem gives an approximation which is formally second- order accurate in space; see [15].Accordingly, a second-order accurate in time Runge-Kutta meth- od must be used to discretize our ODE; see [15, 29   and 30].Finally, a local projection AIIh is applied to the intermediate values of the Runge-Kutta dis- cretization in order to enforce nonlinear stability.
The general definition of the DG method in the case of a scalar u can be found in [15].To define the method in our case, we simply have to apply the procedure for the scalar case component by component.When the formal integration by parts is carried out after the equations are multiplied by test functions, the flux appears as a term of the boundary line integrals.The upwinding is applied in the definition of the flux representation.A local Lax-Friedrichs flux is typically employed.

SIMULATION OF THE MESFET
A benchmark MESFET is displayed in Figure 1.We present some recent results of our simula- tions, particularly with respect to symmetry and symmetry-breaking (see [11]).FIGURE 2 The 1D model with spherical symmetry assumption, in comparison with the 2D MESFET results, at vbias=0.5V, 1.0 V, 1.5 V and 2.0 V.The concentration n.Next, we show the result of using the 1D model with a spherical symmetry assumption, to approximate the 2D MESFET shown above.We take our 1D domain from r=0.025 to r=0.1, measured from the top middle point at (x,y) (0.3, 0.2) downward.The boundary conditions for the concentration n, the temperature T and the potential 4 are prescribed, using the values of the 2D simulations; the boundary condition for the velocity is floating (Neumann).In Figure 2,  we show the comparison, for the concentration n, of the 2D MESFET result with the 1D model assuming spherical symmetry, at vbias= 0.5 V, 1.0 V, 1.5 V and 2.0 V. We can clearly see a qualitatively correct agreement.This is very promising since it means that other quantities (such as T and ) which are not spherically symmetric have minimal effect on the concentra- tion through the nonlinear coupling of the equations.
The concentration n(x, t) is obtained by n(x, t) f (x, u, t)du.(8.4) Also, the electric field E(x, t) is obtained by solving the coupled potential equation, E(x, t) -x, (eCx)x e(n nd), with the boundary conditions

MOBILITY CALIBRATIONS
This work describes results presented in [9], but not previously published.We demonstrate, by comparison of the kinetic and hydrodynamic model, the critical role played by the mobility functions.The device we consider is the one di- mensional GaAs n +n n + structure of length 0.8 tm.The device used is as follows: x E [0, 0.8]; the doping is defined by nd(x)=lO6/ktm 3 in 0<x<0.175and in 0.625<x<0.8,and by ha(x) 2 103/l,tm 3 in 0.225 < x < 0.575, with a smooth intermediate transition.This is exactly the device used in Baranger and Wilkins [6], except for a smooth transition of width 0.05 l, tm at the junctions.
We rewrite the kinetic model for the linearized version we consider.The one-dimensional kinetic model can be written as follows: Of(x, u, t) Of(x, u, t) and the relaxation parameter T is computed by m# # is the mobility and we have the following char- acterizations.
More specific conditions for various models are listed below.
1.For the kinetic model (8.1): The velocity space is artificially cut at -a _< u < a (8.13)where we monitor to ensure that f(x, u, t) is always very small at the boundary u +a for the final steady state results.We learned that it is more than enough in all our runs to use a 3.5.Larger values of a are also used to verify that the results do not change in the pictures.
At x 0, take f (O, u, t) na(O) M(u) (8.14) if u >_ 0, and no boundary condition (extra- polation of the numerical solution from inside the domain to the boundary) if u < 0. Also take 4(0, t) 0. At x 0.8, take f(0.8, u, t) ha(0.8)M(u) (8.15) if u < 0, and no boundary condition (extra- polation of the numerical solution from inside the domain to the boundary) if u > 0. Also take (0, t)= vbias.At u=-a and u=a, take no boundary condition (extrapolation of the numerical solution from inside the domain to the boundary).
We compare the simulation results of the hydrodynamic (HD) model, with various assump- tions on # as given in (8.8), (8.9), (8.10) and (8.11), with the kinetic simulation results obtained with the same mobility assumptions.Figure 3 shows the results of concentration n.We can see that the results between the two models are very similar for # 4 and # =/Z(nd) as given by (8.10), but are quite different for #=0.75 and #=#(E).Alternative simulations were carried out in [19].9. SIMULATIONS OF THE QHD MODEL and the contacts (source and drain) to enhance negative differential resistance.The current-voltage curve for the resonant tunneling diode is plotted in Figure 4 for A V increasing from 0 volts to 0.22 volts (upper curve) and decreasing from 0.22 volts to 0 volts (lower curve).Note that hysteresis occurs predominantly in the region of negative differential resistance.The physical mechanism for hysteresis is that electrons "see" a different potential energy due to different accumulated electron charges in the diode when the applied voltage is decreasing than when the applied voltage is increasing.
One of the encouraging features of the QHD model is the recovery of two fundamental properties of quantum transport as observed in quantum structures such as the resonant tunneling diode.These are negative differential resistance and hysteresis.Below, we indicate in Figure 4 the current-voltage curves reproduced from [14].
To exhibit hysteresis, we simulate a GaAs resonant tunneling diode with double A10.3Ga0.7As barriers (the barrier height /3 0.209eV).The doping density na 1018 cm-3 in the n + source and drain, and na 5 x 1015 cm-3 in the n channel.The channel is 250 long, the barriers are 50 wide, and the well between the barriers is 50 wide.
22].They are derived as the first three moments of the Boltzmann transport equation, Of eE.vuf C. (1.1) 0--+ u VXf m