Monte Carlo Analysis of Anisotropy in the Transport Relaxation Times the Hydrodynamic Model

This paper investigates the anisotropy properties of the relaxation times used in the hydrodynamic model. To this purpose a calculation of the collision term of the Boltzmann Transport Equation is performed by means of a Monte Carlo code accounting for the proper band structure and scattering features. Numerical results for electrons in silicon are shown at 77 and 300 K for E < 00> and E (111).


INTRODUCTION
The hydrodynamic model is a widespread tool for the analysis and design of advanced semiconductor devices.The model is made of the moments of ranks 0 through 3 of the Boltzmann Transport Equation (BTE), constituting a set of two second-order non lin- ear PDEs in the real space [1 ].Its solution provides, in addition to the concentration of the carriers, their mean velocity, energy, and energy flux.In the deriva- tion of the model it is possible to incorporate the effects of the collisions into a set of generalized relax- ation times for the moments of ranks 1, 2, and 3, this yielding three tensors of rank 2, 4, and 6, respectively.The rank is then decreased by taking the trace of the corresponding PDEs [1 ].Up to this point the deriva- tion is carried out with no loss of generality but, in the subsequent solution scheme, some approximations are introduced.In particular, the above tensors are supposed to reduce to scalars.This assumption is in principle not justified under non-equilibrium condi- tions.In fact, two types of anisotropy effects can be relevant.The existence of an electric field disrupting the cubic symmetry of the crystal is expected to reflect into some degree of anisotropy of the distribu- tion function, whose moments over k then acquire a tensor nature.Furthermore, the electron populations around different minima of the band, equivalent by crystal symmetry, show non-equivalent dynamical and scattering properties due to the different effect of the electric field.In order to investigate anisotropy effects it is necessary to determine a full solution of the BTE, incorporating those features of the band structure -such as the presence of six nonspherical, nonparabolic valleys in silicon, and scattering mecha- nisms-that are relevant for a correct description of the problem.To this purpose, a Monte Carlo (Mc) simula- tor for electron transport in S has been used, account- ing for the full 3D electron dynamics in k space for a homogeneous system [2] to obtain results in a number of physical conditions that are most relevant for the applications.In particular, the intermediate range of the electric field (5 < E < 100 kV/cm along (100) and 111 ), at 77 and 300 K) has been investigated, where anisotropy effects are known to be relevant [2].When transport is analysed at much higher electric fields and electron energies (typically above 100 kV/ cm and eV, respectively), the anisotropy effects become progressively less relevant because the carri- ers spread more uniformly over larger and larger por- tions of the Brillouin zone [3].

THEORY AND IMPLEMENTATION
The collision terms for the moments of rank 1, 2, and 3, after calculating the trace as indicated in the Introduction, are expressed by means of three generalized relaxation-time tensors.The definitions of the latter are e(k) f ",f d3k z;ln(w weq), fe;,,e. (2)  d3k In the above, f is the distribution function and n-pfd3k the electron concentration.The inverse total scattering rate z is defined by 1/'c BH(k, k') d3k', the scattering-in and -out terms of the collision integral by f/'c IB f(k')H(k', k) d3k and f/'c, respectively.Quantities w, v, and P are the mean energy, velocity, and energy flux of the electrons, respectively, and are the averages on f of the electron energy e, group velocity u, and energy flux ue.Finally, "c w is the scalar energy-relaxation time, while the diagonal components of the rank-2 tensors p, q are the relaxation times for momentum and energy flow.If the electric field is given, "p, ;Cq "c w, are actu- ally the only "true" coefficients of the equations of rank 1, 2, and 3. On the other hand, unless heavy sim- plifications are made in their structure, a realistic analysis of such coefficients brings back to (1) and (2) anyhow, this calling for a full solution of the BTE.In particular, the analysis of (2) is of interest to deter- mine the anisotropy effects.It is worth observing that a direct evaluation of (2) provides a statistically sound information only in the direction parallel to the elec- tric field, that is, where the distribution function is drifted with respect to the equilibrium one.As a con- sequence, the generalized relaxation times "c w, "Cp33, and q33 were calculated directly from (1) and ( 2), since the average quantities therein are different from zero.Other components of ;Cp and ;Cq can be evalu- ated (see below) using a set of microscopic correla- tion functions obtained as described in [4].The MC program used for the evaluation of the integrals of ( 1) and ( 2) is a fully three-dimensional simulator of the electron dynamics in k space for a homogeneous sys- tem, and includes all six ellipsoidal nonparabolic val- leys as described in [2].Six intervalley-phonon scattering mechanisms (three of f type and three of g type) have been considered, together with acoustic intravalley scattering with the correct energy dissipa- tion and phonon equilibrium distribution.The set of phonon parameters, along with the other silicon con- stants used in the calculations, are listed in [5].
According to the theoretical description of the phonon-scattering mechanisms included in the simulation (described by the usual deformation-potential cou- pling), the inverse total scattering rate "c is a function of the carrier energy and can be evaluated analytically at the beginning of the MC procedure.As for the eval- uation of the "in" and "out" terms of the collision integral, they were recorded during the simulation over a fixed grid in k space, before and after each scattering event.

RESULTS
Numerical results were obtained at 77 and 300 K as functions of electric field strength and direction.For illustrative purposes Fig.
shows the result for the collision term of the BTE in the isotropic Herring- Vogt space of the MC simulation [2] as a function of the component k z of the k vector along the field direc- tion and the modulus of the orthogonal component k+/-.Due to the cylindrical symmetry imposed by the electric field, this grid is the most convenient for the numerical calculations.The physical case T 77 K, T=77 K E=IO kV/cm along <100> FIGURE E 10 kV/cm along the (100) direction is shown.For this case the two valleys with principal axes along the direction of the field (cold valleys) show a dynamical behavior and electron distribution different from the four valleys with principal axes orthogonal to the % 0.5 ., (1.4 0.13 Energy-relaxation time field (hot valleys).As a consequence, electrons in cold and hot valleys show different transport proper- ties and bring different contributions to the collision term.The drift effect of the electric field can be seen in the fact that the difference between the "in" and "out" scattering terms is different from zero.If non- equilibrium transport situations are considered with an electric field along the (111) direction, different values are obtained for both the single valley contri- butions and the total collision term.
Figs. 2-4 show the energy-relaxation time and the components p33, q33, respectively, evaluated using equations (1,2), as functions of the electric field strength at 77 and 300 K.A general comment valid for the three relaxation times is that their values seem to depend very weakly on the electric field direction.
At the lowest considered temperature only "w for E (100) is higher than the corresponding value for the case E (111 ).This seems to suggest that the ani- sotropy with respect to the electric field direction observed in many microscopic properties of the system (e.g., the drift velocity [2]) is mainly related to the different effective masses along the field in the different valleys, and not to different relaxations.
Finally, Fig. 5 shows the transverse momentum-and energy-flux relaxation times :pll, 1;ql 1, as functions of the electric field strength along <111) at 300 K together with the corresponding longitudinal components.It is seen that the transverse relaxation times appear to be always larger than the longitudinal relax- ation times.
Maria Cristina Yeeehi received the degree and the PhD in Electrical Engineering from the University of Bologna, Italy, in 1990 and 1995, respectively.Since 1990 she has been working at the Department of Elec- tronics (DEIS) of the same University in the field of the numerical simulation of semiconductor devices, with emphasys on the optimization techniques for process design and on the simulation of optical sen- sors.In 1992-93 she visited the IBM T.J. Watson Research Center of Yorktown Heights, NY, working on higher-order transport models in semiconductors.
In 1994 she visited the Integrated System Laborato- ries of the ETH, Zurich, working on the deterministic solution of the Boltzmann transport equation.
Massinio Rudan was born in Bologna, Italy, in 1949.He received a degree in Electrical Engineering in 1973 and a degree in Physics in 1976, both from the University of Bologna, Italy.He joined the Department of Electronics (DEIS) of the University of Bologna in 1975, where he investigated the physi- cal properties of the MOS structures and the problems of analytical modelling of semiconductor devices.From 1978 he has been teaching an annual course of Quantum Electronics in the Faculty of Engineering of the same-University.Since 1983 he has been working in a group involved in numerical analysis of semicon- ductor devices, acting as Task leader in a number of EEC-supported Projects in the area of CAD for VLSI.In 1986 he has been a visiting scientist, on a one-year assignment, at the IBM Thomas J. Watson Research Center at Yorktown Heights, NY.In 1990, he was FIGURE 2FIGURE3