Three-Dimensional Hydrodynamic Modeling of MOSFET Devices

The hydrodynamic (HD) model of semiconductor devices is solved numerically in three-dimensions (3-D) for the MOSFET device. The numerical instabilities of the HD model are analyzed to develop a stable discretization. The formulation is stabilized by using a new, higher-order discretization for the relaxation-time approximation (RTA) term of the energy-balance (EB) equation. The developed formulation is used to model the MOSFET.


INTRODUCTION
Although the hydrodynamic model is not new, its numerical solution by computer simulation for multi- dimensional problems remains difficult.The problem is nonconvergence, due to oscillations in the iterations or divergence.The problem can be divided into two parts, which are non-convergence of the EB equation itself and of the overall system.Furthermore, a suc- cessful solution mesh in unusually difficult to con- struct.More work on numerical methods for the HD model needs to be directed to making it stable and routinely convergent.

HYDRODYNAMIC MODEL
The equations of the HD model consist of five cou- pled partial differential equations (only equations for electrons are shown; hole equations are analogous), q (n-p-D), V" Sn Jn E n Rwn, where eq. ( 1) is the Poisson equation, (2) is the elec- tron current-continuity (CC)equation, and (3) is the electron energy-balance (EB) equation.In the EB equation, the first and second terms on the right-hand side (RHS) are called the input energy and relaxa- tion-time approximation (RTA) terms.The unknowns in the above equations are the electric potential , car- rier concentrations n and p, and (scalar) carrier tem- perature T n.The material parameters are dielectric permittivity s and net doping concentration D NNTt.The net recombination R, average car- rier energy w n, energy relaxation time Xnw, and aver- age thermal equilibrium carrier energy w 0 (3/2)kBTL, with Boltzmann's constant k B and the lattice temperature Tt, can be expressed in terms of the unknowns.
The vectorial quantities electric field E, current density Jn, and energy flux S n appearing in the model equations above are related to the unknowns as fol- lows: Jn --qlunnV + qDnVn + npnkBVTn, ( S, J--(w, + kBTn) 2(kB/q)2qnt,TnVTn, ( -q where the mobility [.t n and diffusivity D n have been used.

Discretization
The device geometry is discretized by introducing the computational mesh, which is a set of nodes and edges connecting them.Given certain mesh proper- ties, a Voronoi cell, with volume V i, associated with node can be constructed by taking the smallest vol- ume enclosed by the planar perpendicular bisectors of the edges leading from the node to its neighboring nodes j [1 ].The equations at a node are discretized by integrating each equation over the Voronoi cell asso- ciated with that node.
The oscillation problem can be controlled some- what by locating the problem spot and increasing the mesh density in its vicinity.This has the effect of removing or diminishing the oscillations.This tech- nique has been adequate for simulations using 2-D meshes, where mesh sizes of the order 3,000-5,000 nodes can be used.However, the mesh design proce- dure is tedious, time-consuming, and unsuitable for 3-D simulation.
Efforts to obtain robust, general-purpose solution algorithms have focused primarily on new discretiza- tions of the vectorial fluxes J and S. The stability of the EB equation has been enhanced by the application of the Scharfetter-Gummel (SG) discretization of the differential operator [2][3][4].This type of discretization is known to produce nearly diagonally-dominant coefficient matrices, which in practice converge for a wide range of meshes and initial guesses.The discretization of Ref. [4] was implemented and tested in the 3-D device simulator SIMASTER [5].The equations were solved within the context of the block-Gummel loop.It was found that oscillations in the iterations, illustrated in Fig. 1,developed between the electron temperature at a small number of nodes oscillating between high and low temperatures; diver- gence sometimes occurred if the range of variation was not restricted.Frequently the problem was exac- erbated by an oscillation in space, involving both car- rier concentration and temperature, as illustrated by Fig. 2.These are primarily multidimensional effects, where a low n-density region near a high density region carries the oscillations.Thus, the discretization of Ref. [4] resolves the non-convergence issues of the CC and EB equations individually; however, difficul- ties remain in the convergence of the overall system when using a block-Gummel solution method.Recognizing that the convergence of the entire system of equations was the problem, some researchers sug- gested that the input energy term is the root of the nonconvergence problem [6, 7].In particular, the numerical problems were thought to be caused by the implicit dependence of the current density on the car- rier temperature [7].The implicit dependence can be made explicit by expanding the current into its three components, as in eq. ( 5).In Ref. [7], a special quad- rature rule for the term involving thermal diffusion was suggested, which resulted in improved stability and convergence speed.However, by not consistently using the SG-discretized expression for current den- sity, this scheme underestimates the input energy.Here, a vector identity is used in the discretization, following [4, 6, 8].In this work, the effect of the dis- cretization of the RTA term on convergence is investi- Electron FIGURE 2 Correlated spatial oscillations of carrier concentration and temperature in a 3-D simulation.The quantities are plotted along a line in the z (width) direction at the location of maximum temperature gated.The RTA term appears on the RHS of the EB equation, as n(w n Wo)/wn.In the box-integration scheme, its discretization proceeds from n .dVnvndV Twn mn 3 kB fv,.n(Tn--TL)dV (7)   + The second term on the RHS of the above equation is the dominant term; no special treatment of the first term is necessary.The simplest quadrature rule for the second term approximates the integrand as a constant over V i.The convergence of the EB equation is rou- tine; however, the convergence of the block-Gummel loop is still a problem.It turns out that approximating the electron concentration as a constant is a bad choice, since it varies exponentially along mesh lines.
A higher-order quadrature rule should be used which preserves the stability of the EB equation and stabi- lizes the block-Gummel iterations.
A more accurate approach would be to partition the integration over the Voronoi cell into sections, where the subvolumes Vij are.defined by the polygonal pyra- Within each sub-volume, the approximation is made that the integrand only varies along the mesh edge, which is perpendicular to the base of the pyramid.The variation in T n is taken to be linear, and the varia- tion of n exponential; after simplification the integration results in nW"-WdV ---(Tni-TL) Z [, ndV The above equation can be split and included in the LHS coefficients.In the first term on the RHS of ( 9), since the integral is always positive, the first term contributes to the diagonal dominance of the coeffi- cients.However, the second term will always degrade this property.Fortunately, it can be shown that the second term is (h2) and smaller than the first term, which is of ((h).The integrals in eq. ( 9) can be com- puted analytically [5].
The new discretization of the RTA term stabilizes the block-Gummel iterations.Extensive simulations using this discretizadon of the RTA term have demonstrated that this quadrature rule is stable in the EB equation, and also stabilizes the block-Gummel iterations.
This alternative discretization of the RTA term requires changes to the discretization of the other equations to avoid inconsistencies.In the Poisson equation, there appears an integration of charge den- sity over the Voronoi cell.The higher-order quadrature rule should also be applied to this integration.In the SG-discretization of the CC equation, a functional form of charge density along mesh lines is calculated which deviates from the exponential form [4].The impact of these potential inconsistencies may be small, as the order of the truncation error is the same.The complete analysis of this matter will be published elsewhere.
The system of discrete equations is solved by using the block-Gummel algorithm.The block-Gummel algorithm decouples the model equations into two groups, which are the drift-diffusion (DD) block, con- taining the Poisson and CC equations, and the HD block, containing the EB equations.Within each block, the variables from the other block are assumed known from previous iterations.In addition, the Gummel method is used to decouple the individual equations within each block.The discretized model equations are solved by using the fixed-point iteration method [9, 10]. 3. APPLICATIONS A 0.5 }.tm LOCOS MOSFET was simulated in 3-D using the HD model.The device geometry is shown in Fig. 3.In the figure, note that the semiconduc- tor-oxide interface is nonplanar; several others figures will plot various quantities at this surface.Fig. 4 shows a typical plot of the maximum relative error of three equations vs. the block Gummel iteration number.At first, the n-EB error is small, since the ini- tial guess uses a 2-D HD simulation which is spread across the 3-D domain.Overall, the convergence is smooth and stable.Typical simulation times on a Sun Sparcstation-10 workstation were 8 and 81 min for 2-D and 3-D HD simulations, respectively.
The electron concentration and temperature profiles at the nonplanar interface are shown in Fig. 5.In these figures, the active device area lies where z > 0, with the drain at x > 0.7, source at x < 0.1.The elec- tron concentration spreads out along the interface near the drain junction due to the elevated electron temperature, shown in Fig. 5(b).The elevated elec- tron concentration beneath the isolation oxide harms the device isolation.One technique to control the problem introduces a channel-stop doping, which is a high p-type doping beneath the isolation oxide.The Block Gummel iteration number FIGURE 4 The maximum relative error of all three equations vs. block Gummel number showing successful convergence for the bias condition V D 4 V, V G 2 V.The solution time was 81 min subsurface charge, illustrated in Fig. 6, is somewhat damped, especially near the source.However, the electron temperature still drives the charge under the isolation near the drain.Finally, drift-diffusion model simulations of the same structure, shown in Fig. 7, cannot predict these effects.
In conclusion, this paper has shown that SG-fitting alone is inadequate for stable HD simulations using the block-Gummel method.Adjusting the discretiza- tion of the input energy term does not improve con- vergence.However, the block-Gummel iteration method for the HD model can be stabilized by adjust- ing the discretization of the RTA term.A higher-order discretization was implemented in the simulator SIMAsTER.Applications have shown that the charge Electror charge (a) Electron concentration.

FIGURE
FIGURE Oscillations in the block-Gummel iterations the DD and HD blocks.Often the problem involved mids from the central mesh point to Voronoi cell faces,

FIGURE 3
FIGURE 3  The exterior of the 3-D semi-recessed LOCOS MOS- FET.The mesh contains 24,403 nodes