Hyperbolic Hydrodynamical Model of Carrier Transport in Semiconductors

Hydrodynamical models for carrier transport in semiconductors have been widely investigated by several authors [1, 2, 3, 4, 5, 6 and al.]. Most hydrodynamical models are based on phenomenological assumptions and not derived from a consistent physical model. Here we present a hydrodynamical model based on Extended Thermodynamics which is the appropriate thermodynamic theory of irreversible processes far from thermodynamic equilibrium. The model is obtained from the set of moment equations of the Boltzmann transport equation for electrons at the level of 13 moments. The closure relations are obtained in the framework of Extended Thermodynamics [7, 8] by employing the entropy principle. The balance equations for the macroscopic variables form a hyperbolic system for values of heat flux and stress tensor which are routinely encountered in the simulation of semiconductor devices of interest for industrial purposes. In the present paper we model the productions as relaxation terms, with relaxation times functions of the moments. A suitable numerical scheme, developed in [10], has been successfully applied to the present model. It is based on the splitting method between relaxation and convection and allows an easy implementation in numerical codes. For the convective steps we used the Nessyahu-Tadmor [11] scheme for quasilinear hyperbolic systems in divergence form. It is a second order shock capturing method and does not require the analytical expressions of the eigenvalues and eigenvectors of the Jacobian matrix of fluxes.


INTRODUCTION
Hydrodynamical models for carrier transport in semiconductors have been widely investigated by several authors [1, 2, 3, 4, 5, 6 and al.].Most hydrodynamical models are based on phenomenological assumptions and not derived from a consistent physical model.Here we present a hydrodynamical model based on Extended Ther- modynamics which is the appropriate thermodynamic theory of irreversible processes far from thermodynamic equilibrium.The model is ob- tained from the set of moment equations of the Boltzmann transport equation for electrons at the level of 13 moments.The closure relations are obtained in the framework of Extended Thermo- dynamics [7, 8] by employing the entropy princi- ple.The balance equations for the macroscopic variables form a hyperbolic system for values of heat flux and stress tensor which are routinely encountered in the simulation of semiconductor devices of interest for industrial purposes.In the present paper we model the productions as relaxation terms, with relaxation times functions of the moments.
A suitable numerical scheme, developed in [10], has been successfully applied to the present model.It is based on the splitting method between relaxation and convection and allows an easy implementation in numerical codes.For the convective steps we used the Nessyahu-Tadmor [11] scheme for quasilinear hyperbolic systems in divergence form.It is a second order shock capturing method and does not require the analytical expressions of the eigenvalues and eigenvectors of the Jacobian matrix of fluxes.*Corresponding author.
The latter property is crucial because explicit expressions for the eigenvalues and eigenvectors are not available for the balance equations of Extended Thermodynamics.For the relaxation part a suitable combinations of implicit Euler steps assures a global second order accuracy.
In particular we present the simulation of a + n -n -n silicon diode.
The numerical results shows a uniform accuracy in all the fields (energy density, velocity, particle number, heat flux) within a reasonable degree of tolerance and could be employed for a systematic use as a CAD tool.

MATHEMATICAL MODEL
In the one dimensional case the evolution equa- tions become (see [9] for the complete derivation) On 0 (nv) 0 (1) O---t n +--; Ox nv3 4vp 7vo 8 q ) +-j-+5-+-- In the previous equations, n is the electron number density, v the mean velocity, p nkBT the hydrostatic pressure, q the x-component of the heat flux, cr the xx-component of the viscosity tensor, m* the effective electron mass in the parabolic band approximation and E the xcomponent of the electric field 02 Ox 2 --e(ND NA n), where 4 is the electric potential, No and NA are respectively the donor and acceptor density, e the elemetary charge and e the dielectric constant.
Numerical solutions of the system are obtained by a splitting scheme.Let us consider a system of the form Ov OF(v) (v), Ot Ox with v E R and F: Rm--+ Rm.
The splitting scheme is based on the solution of the two steps: where A dt/dx.The time step At must satisfy a stability condition .maxp(A(v(x, t))) < . (11)ere p (A(v(x,t))) is the spectral radius of the Jacobian matrix, OF This condition will ensure that the generalized Riemann problems with piecewise smooth data at time tn will not interfere during the time step At.
The values of vj./AxandF'j./Ax are a first order approximation of the space derivatives of the field and of the flux, computed from cell averages by using Uniform Non Oscillatory reconstruction (UNO), v'j.MM dj_)v +-MM(Dj_,,Dj), dj+1/2v -1MM(Dy' DJ+'

FIGURE
Average electron energy vs. distance obtained with the hydrodynamical model (continuous line) compared to Monte Carlo results (circles).A similar procedure is used for computing Fj For the relaxation step an unconditionally stable second order scheme can be obtained by analytical integration of the linearized relaxation equation, where linearization is obtained by freezing the coefficients at time tn.The electric potential is computed by a standard procedure.In order to obtain second order accuracy in time we combinine the two steps according to the following scheme where R represents the discrete operator corre- spoding to the relaxation step, C/xt is the discrete operator corresponding to NT scheme, and 79(U)gives the solution to Poisson's equation.
For the boundary conditions see [9].The stationary solution is reached within a few picosends.The numerical results obtained for the stationary case have been compared with the Monte Carlo results, obtained by the DAMO- CLES code [13].
Comparison with MC simulations shows a good overall agreement.In particular the agreement is very good everywhere in the device but in the second junction, where the discrepancy with MC simulations can reach up to 10%.

NUMERICAL RESULTS
As test problem we consider a ballistic diode n + -n + n which models a MOSFET channel.The diode is made of silicon, and the bulk temperature is supposed to be 300 K.The n + regions are 0.1 tm long and the channel is 0.3 lam long.The doping density is N-0.5 1017 cm -3 in the n + regions and ND 0.21016 cm -3 in the channel.
For the electron effective mass in the approxi- mation of parabolic band we use m*= 0.32 me when me is the electron mass.The silicon dielectric constant is given by e ere0, where er 1.7 is the relative dielecric constant and 0 8.85 10 -18 C/ Vlam is the dielectric constant of vacuum.
The initial electron temperature is the lattice temperature To 300 K and the charges are at rest.A bias voltage of Volt is applied, and this determines a charge flux in the semiconductor.

5 FIGURE 2
FIGURE 2 As in Figure fo the average drift velocity vs.distance.