An Efficient Solution Scheme for the Spherical-Harmonics Expansion of the Boltzmann Transport Equation Applied to Two-Dimensional Devices

This paper presents a novel numerical scheme applicable to the solution of the Boltzmann transport equation by means of a spherical-harmonics expansion. This scheme improves the solution at low energies, keeping the desired accuracy in the calculation of the mean quantities while saving a significant amount of CPU time. This is important in view of the applications of the method, since the typical number of nodes to be used in the combined space-energy domain is in the range of 104105.

In recent years the spherical-harmonics expansion (SHE) of the Boltzmann transport equation (BTE) has successfully been tested in a wide range of problems in the field of electron transport simulation 1,2].The method transforms the differential equation in energy and space into a differential equation in space and a difference equation in energy [3].This allows one to connect to each point in space only those nodes in energy coupled via the scattering operator.To achieve a detailed knowledge of the distribution function in energy, a rather fine discretization of the energy domain is necessary.Thanks to the structure of the SHE equation, the energy discretization can be decomposed in several grids, say L, displaced in the energy domain where the solution has to be com- puted.In order to correctly compute the averages 239 (namely concentration, mean energy, mean velocity, etc.), the calculation of the solution on several grids is especially important in the low-energy range.We pro- pose here a new scheme that improves the solution at low energies, keeping the desired accuracy in the cal- culation of the mean quantities while saving a signifi- cant amount of CPU-time.
The SHE method in steady state provides the sec- ond-order difference-differential equation [3]   -g-[N+opfo Nopfo } -0, (1) in the unknown distribution function f0.The equation is solved in a three-dimensional domain (x, y, H) with boundaries defined in space by the boundaries of the device and in energy by Hmin =-qtp(x,y) and Hmax Emax -qtp(x,y), where Emax is the maximum energy of the band system.A prismatic mesh with triangular elements in the (x,y) plane is adopted.The nodes at constant total energy H are uniformly spaced in energy by intervals AH hoop with n integer.
Thanks to the absence in (1) of partial derivatives with respect to H each node is connected along the H direction only to the nodes at (x,y, H + hO3op) via the phonon scattering operator.The resulting algebraic system is thus decomposed in L decoupled subsys- tems.In the solution scheme we present here, eqn.(1)   is solved in the full energy domain [Hmin, Hmax] only once, while several additional solutions (i.e.L 1) are computed in the reduced energy domain [Hmin, rn Hop ], where Hop hcoop q(p(x,y).The terms in eqn.(1) which are external to [Hmin, rn Hop are expressed by interpolating the values previously com- puted on the full domain.By letting rn 10 one can express the values of the distribution function outside the domain by an interpolation scheme that maintains the linearity of the system; conversely, to compute the local solution with rn one has to resort to a more sophisticated interpolation scheme.Both schemes are based on a Lagrange polynomial interpolation per- formed on the logarithm of the distribution function but, while the first one only involves previously com- puted values of the distribution function, the second one involves also the unknowns of the local system.It should be added that, in addition to the loss of linear- ity of the system, a local solution scheme with rn does not ensure the conservation of the electron cur- rent density because of the early truncation of the solution domain.This can be understood by consider- ing the fluxes in a simplified one-dimensional case (Fig. 1).The x axis is in the horizontal direction, hence node and the two adjacent nodes i-1, + have the same total energy, say H i. The H axis is in the vertical direction, hence the energy of node + is Hi+ ho)op.In particular H is assumed here to repre- sent the lowest energy at location i, namely, there is no other node below node i.The incoming and outgo- ing contributions to the fluxes of f0 for node are also represented in the figure.The dotted arrows indicates the fact that the second term in brackets of ( 1) is miss- ing from the system matrix when the scheme with rn is used, because no equation such as (1) is writ- ten for node +.Still considering the one-interval scheme, it is worth adding that the value of fi used in the calculation of the flux i not fully consistent with foi.This leads to a non-conservation of the car- rier distribution flux which yields a non-conservation of the current density at the device contacts.In theory, this phenomenon also takes place when the scheme with rn 10 is adopted; however, it occurs at rela- tively higher energies, hence its effect on the distribu- tion function and its derivative is not relevant.

RESULTS
We present here, for the case rn 10, the comparison of the results of this solution scheme (consisting of full solution and 9 local solutions) with those of a full solution computed on 10 grids in [Hmin, Hmax].The test device is a 2D-MOS transistor with gate length 0.33 lttm.The device has a bulk doping concentration of 1017 cm-3, a Gaussian p-type channel implant with a 2 1017 cm -3 peak located at the surface and a 0.04 m standard deviation.The source and drain junction depth is 0.11 m and the oxide thickness is 7 nm.The simulation has been performed with VGS VDS--3 V, using a non-uniform mesh with tri- angular elements.Figs. 2 and 3 present the compari- son of the electron concentrations and the normalized electron mean energies along a section of the device parallel to the Si-SiO 2 interface, and Fig. 4   FIGURE 2 Electron concentration along a section parallel to the Si-SiO interface computed with the SHE full-domain and the local schemes the drain currents.For comparison the current com- puted by HYDRO-HFIELDS [4], using a mobility model consistent with SHE, is also reported in Fig. 4. In conclusion, the scheme based on the solution of L subsystems displaced in the full-energy domain can be replaced without loss of accuracy by a scheme which considers a single complete solution and (L-1) to its structure, the new method is intrinsically more efficient: in the example shown here, a reduction fac- tor of about 5 in the CPU-time has been achieved.
Output characteristic computed HYDRO-HFIELDS and two solution schemes of the SHE with Biographies Maria Cristina Vecchi 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 emphasis 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, Ziirich, working on the deterministic solution of the Boltzmann transport equation.
Jan Mohring received the degree in industrial mathematics in 1992 from the University of Kaisers- lautern, Germany.Since then he has been working towards the Ph.D. degree in mathematics at the Uni- versity of Kaiserslautern and, for four months, at the Department of Electronics at the University of Bolo- gna, Italy.His thesis will deal with numerical improvements and a mathematical consolidation of a deterministic method to solve the Boltzmann equation in device simulation.
Massimo 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

FIGURE
FIGURE Sketch of the fluxes relative to node in the case of the one-interval local solution compares

FIGURE 3
FIGURE 3Normalized electron energy along a section parallel to the Si-SiO 2 interface computed with the SHE full-domain and the local schemes