A New Self-Consistent 2D Device Simulator Based on Deterministic Solution of the Boltzmann, Poisson and Hole-Continuity Equations

LDD MOSFET simulation is performed by directly solving the Boltzmann Transport Equation for electrons, the Hole-Current Continuity Equation and the Poisson Equation self-consistently. The spherical harmonic expansion method is employed along with a new Scharfetter-Gummel like discretization of the Boltzmann equation. The solution efficiently provides the distribution function, electrostatic potential, and the hole concentration for the entire 2-D MOSFET.


I. INTRODUCTION
The spherical harmonic (SH) method for solving the BTE represents a new middle option for device simu- lation which gives more information than the Hydrodynamic model and is more efficient than the Monte Carlo appraoch.It can provide the momentum distri- bution function for spherical band models [1,2], and is free from statistical noise.It has been shown to agree with Monte Carlo simulations which employ the same transport model [3,4,5], while requiring a minute fraction of the CPU time to evaluate.It has been used as a post-processor for 2-D MOSFET sim- ulation [6,7].In addition, the approach has been gen- eralized in 1-D and 2-D to include an arbitrary number of spherical harmonics [7,8].
In this paper we extend the SH approach by self-consistently solving the BTE, Poisson and Hole-Current Continuity Equations for a realistic sub- micron 2-D LDD MOSFET structure.To solve the BTE, we first use the SH approach to transform it into a 3-D equation Which is tractable for numerical solu- tion.To facilitate stability, we have developed a new 251 Scharfetter-Gummel type discretization for the BTE.To solve the overall coupled nonlinear system, while simultaneously overcoming computational restric- tions imposed by the 3-D BTE, we have developed an algorithm which combines iterative and direct matrix techniques.We use the simulator to calculate the dis- tribution function for the entire MOSFET, from which we determine channel current, substrate currents and average quantities.Our results are in qualitative agreement with experiment and/or other numerical techniques.
The remaining parameters have their standard defini- tions.)Our collision integral in the BTE accounts for phonon and ionized impurity scattering, impact ioni- zation and recombination.For ionized impurity scat- tering we use the Brooks-Herring formulation; while impact ionization is accounted for with a random-k approach.For recombination we have derived an expression which is analogous to the SHR approach.We use the acoustic, optical and intervalley phonon deformation potentials and band structure (up to 3eV) which were originally developed for MC simulations of [1,2].

III. METHOD OF SOLUTIONS
The Poisson, Boltzmann and Hole-Continuity Equations give rise to a coupled nonlinear system.One nonlinearity arises due to the coupling in the field-term of the BTE where there is a product between the gradients of the potential and the distri- bution function.Another nonlinearity arises due to the product of the hole concentration and the gradient of potential in the Hole-Current Continuity Equation.To treat these nonlinearities, we solve the BTE, Poisson and Hole-Continuity Equations self-consistently using a decoupled or Gummel-like iteration method.Within each Gummel loop, the Poisson Equation, the BTE and the Hole Equation are solved independently.The flow chart in Figure shows the overall algorithm employed, below we describe the method of solution of the individual equations.

Poisson and Hole-Continuity Equations
To solve the Poisson equation, we first discretize with finite differences, and then directly solve using sparse matrix algebra.
To solve the Hole-Current Continuity Equation, we first transform it into self-adjoint form using Slot- boom variables.We then apply the Scharfetter-Gum- mel discretization to the transformed equation.The resulting diagonally dominant system is then scaled, and solved iteratively using the fixed point approach [9].

Boltzmann Transport Equation
To solve the BTE, we begin with the generalized SH technique [7].It was shown in [7], that by substituting the SH expansion, f(}, ) _,lmf' (r, e)Y'(O, ) into the BTE, while taking advantage of the recur- rence and orthogonality relationships between spheri- cal harmonics, the BTE can be transformed into the following system of equations: i--1 at J. It is interesting to note that each equation has an iden- tical form, the system can therefore be automatically generated to arbitrarily high order and then be solved numerically.
Drift Diffusion Model _-O(x,y), n(x,y), p(x,y) Maxwellian Distbribution Initial Guess---fO(x,y,e) o (x,y), n(x,y), p(x,y) Poisson Equation----"+l(x,y) !(x,y), n" (x,y) P"(x,y) (x,y), (x,y,e) It has been demonstrated that truncating the SH expansion at 1 leads to results which agree with Monte Carlo calculations [3,4,10].We therefore trun- cate the expansion at 1.This leaves an expansion of the first four spherical harmonics (Y0 , y]-l, y0, Y1 ).Next, the H transformation is employed [11], where the total energy H is used as the independent energy variable as opposed to the kinetic energy e (H e- Next, to overcome numerical difficulties arising from rapid variations in the distribution function, we have developed a Scharfetter-Gummel (SG) like dis- cretization for the BTE.To achieve this discretization, we perform a piecewise analytical integration of the BTE between mesh points, and then discretize the resulting integrated form.After truncation and a sig-nificant amount algebra, followed by the SG discreti- zation, we arrive at the following discretized form for the BTE: 2(Fi+ljm-Fijm)

2(Fijm-Fi-ljm)]
Tli+1/2jm hi(hi-Fhi-1) -Tli-1/2Jm hi-l(hi+hi-1) rl 2(/j+ m Fijm) Where Fi,j,m F(x i, yj, mail), n the closest integer to AH' Yijm (mAH + q(x i, Yi)) and i+ l,]m--i]m T + jm In i + ljm In ijm Since equation ( 5) is a finite difference equation in 3-D, it actually represents a system of N x x Ny x N H 100, 000 discrete equations.It is usu- ally not practical to solve such a large system directly with sparse matrix solvers.Instead, we use an itera- tive technique, which minimizes memory usage and facilitates fairly rapid evaluation of the 3-D system.While the aforementioned Gummel scheme appears straight forward, actual implementation is particularly complicated because our domain is curvelinear.The lower and upper H-space domain boundaries of the BTE-Poisson system are the curved surfaces H =-q(x,y) and H e-q(x,y), respectively.These surfaces can vary very rapidly, especially in the drain region.Since the potential is updated at every Gum- mel iteration, these surfaces, and therefore the simula- tion domain, change at each Gummel step.Furthermore, in general H-space nodes will not lie on these surfaces.To help overcome these obstacles, careful monitering of the domain boundaries and interpolation is employed.In addition, an SOR-type damping term is used in the Poisson equation to ensure the simulation domain does not vary too rapidly between iterations.

IV. RESULTS
We demonstrate this device modeling method with calculations of a realistic 0.5 gm LDD MOSFET.The device has a junction depth of 0.11 gm, source and drain doping of N D 102/cm3, LDD region doping of N D lO18/cm3, and substrate doping of N A 1016 cm3, We show results for an applied bias of Vcls= 3 V and Vgs 3 V.Fig. 2 shows the energy distri- bution function along MOSFET channel, Fig. 3 shows the energy distribution function along a line which is parallel to the interface and goes through the depletion regions at the bottom of the source and drain.
To help demonstrate the accuracy of the method, we plot various moments the distribution function, and show the results have the same qualitative fea- tures as other 2-D device models and/or experiment.Fig. 4 shows the electron concentration, which is obtained by integrating the distribution function over energy-space.In Fig. 5 we show the average electron velocity determined by calculating the first moment of the distribution function.Fig. 6 shows the concen- tration of electrons or holes generated per second from impact ionization throughout the device.Fig. 7 shows the I-V characteristics for a standard MOSFET, which are calculated directly from the distribution function without empirical mobility models.Vds(V) FIGURE 7 I-V characteristics for the LDD MOSFET

( 4 )
where v() m/m', ' d/ds, repre- sents the dispersion relation, and is energy.The rais- ing and lowering operators which we developed to relate the coefficients are:

FIGURE 2
FIGURE 2 Electron distribution function in the channel of LDD MOSFET

FIGURE 3
FIGURE 3 Electron distribution function along the depletion region in source and drain

FIGURE 5 FIGURE 6
FIGURE 5 Impact ionization rate