Convergence Properties of the Bi-CGSTAB Method for the Solution of the 3 D Poisson and 3 D Electron Current Continuity Equations for Scaled Si MOSFETs

As semiconductor technology continues to evolve, numerical modeling of semiconductor 
devices becomes an indispensible tool for the prediction of device characteristics. 
The simple drift-diffusion model is still widely used, especially in the study of subthreshold 
behavior in MOSFETs. The numerical solution of these two equations offers 
difficulties in small devices and special methods are required for the case when dealing 
with 3D problems that demand large CPU times. In this work we investigate the 
convergence properties of the Bi-CGSTAB method. We find that this method shows 
superior convergence properties when compared to more commonly used ILU and SIP 
methods.


INTRODUCTION
The demand of having more transistors on a single chip continues to push MOSFETs toward shorter channel lengths.In order to maintain the long- channel behavior and suppress undesirable short- channel effects, channel doping levels must increase, whereas oxide thickness, junction depths and isolation spacing decrease.To make final threshold voltage adjustments and to eliminate the punchthrough effect, ion implantation is used at several *Corresponding author.steps in the fabrication process of these state-of- the-art devices.The stochastic nature of ion implantation across the wafer leads to significant deviation in both the number of dopants and the arrangement of these dopant atoms from device to device [1,2].While a continuum doping model may be a valid assumption for large devices, it fails in the case of deep-submicrometer devices.For example, in MOSFETs with effective channel lengths and widths on the order of 0.1 lam and channel doping of 51017cm-3, the number of impurity atoms in the inversion layer is less than 100.Thus, differences in the arrangement of dopant atoms, as well as fluctuations in the number of dopants under the gate, lead to significant variation of the off-state leakage current in these state-of-the-art devices.
The high doping levels used in submicrometer devices lead to large gradients in the potential and charge within the device, which requires a very fine mesh.This makes some of the commonly used iterative methods for the solution of both 3D Poisson and 3D electron current continuity equations rather inefficient because of poor con- vergence properties.To overcome this, we use a variant of the Conjugate Gradient Squared meth- od (CGS), the so-called Bi-CGS Stabilized (Bi- CGSTAB) method [3].The convergence properties of this method are found to be superior to those of the Incomplete Lower-Upper (ILU) factorization method [4].A variant of the ILU method, Stone's [5] Strongly Implicit Proced.ure(SIP) was also investigated, but found to be comparable to ILU.Finally, we discuss our 3D simulations of the subthreshold characteristics of 0.1 lam gate-length MOSFETs with various gate-widths.

APPROACH
The geometrical and structural parameters of the device used in our numerical analysis are illu- strated in Figure 1.The depth of the discrete doping region is 0.14m which, for substrate doping N, =51017cm -3 and oxide thickness to= 3nm, corresponds to more than twice the depth of the depletion region.The 3D Poisson and electron current continuity equation are solved assuming Fermi-Dirac (degenerate) statistics with Dn and # constant.
The boundary conditions for , n and p are as follows: Contact regions are assumed to be charge neutral so the electrostatic potential and the carrier concentrations n and p at the contacts are prescribed by the usual Dirichlet conditions.At the free surfaces and artificial boundaries, the perpendicular field and the perpendicular current are set to zero.We use the decoupled scheme [6] in which, the linearized Poisson equation is repeat- edly solved for the electrostatic potential p.The resulting values of are then used to obtain the improved values for n and p.This loop corre- sponds to one Gummel cycle.The above process is repeated until self-consistent values of , n and p are obtained.
The application of finite-difference techniques to both 3D Poisson and 3D electron current con- tinuity equations leads to algebraic equations having a well-defined structure, represented by Ax b.The most suitable methods for solution of these two equations are direct methods, but the computational cost becomes prohibitive as the number of equations increases.This has led to iterative procedures that utilize the structure of the coefficient matrix.The ILU methods, and the use of preconditioning together with conjugate gradi- ents, provide a significant increase in the power of iterative methods.
Within incomplete factorization schemes [4], the original matrix A is decomposed into a product of lower and upper triangular matrices L and U respectively.This is achieved by modifying matrix A through the addition of a small matrix N. Thus, one solves the modified system LUx=b+Nx by solving successively the matrix equations LV b + Nx and Ux, where is an auxiliary matrix.
The superfluous terms of N affect the convergence of the ILU method.Stone [5] suggested the introduction of partial cancellation which mini- mizes the influence of these additional terms and accelerates the rate of convergence.This partial cancellation procedure is implicitly included in the calculation of L and U. To provide rapid convergence, the fraction that is canceled is varied on successive pairs of iterations.
The basic CG algorithm loses its applicability when the resulting system of equations is not SPD, as it is the case with the electron current continuity equation.In such circumstances, the best alter- natives are Lanczos-type algorithms.In recent years, the CGS method due to Sonneveld [7] has been recognized as an attractive transpose-free variant of the Bi-Conjugate Gradient (Bi-CG) iterative method.The Bi-CGSTAB method due to Van der Vorst [3] is a variant of the CGS algorithm which avoids squaring of the residual polynomial.The convergence behavior of this method is smoother because it produces more accurate residual vectors and, therefore, more accurate solutions.To reduce the spectral conditioning number and improve the distribution of the smallest eigenvalues, one applies preconditioning techniques together with conjugate gradient solvers.Successful preconditioning matrix can be obtained by using ILU factorization [8]: If L and U are the strictly-lower and the strictly-upper triangular parts of A, then the preconditioning matrix is KILU (L + l)I-i (U --]), (1)   where diag(KiLu)=diag(A).Once the diagonal is computed, one scales the original matrix A / I-l/2A] -1/2 diag(/) + L + .(2) The preconditioning matrix for this symmetrically scaled matrix is of the form where I is the identity matrix.The Bi-CGSTAB method is now applied to the preconditioned system (1 + 1) -1 (I + [)-1 (] _]_ i)-l, (4) where --l/2b and x--l-l/2(l + )-1.In the calculation of the product (] + I) -(I + I])-, significant amount of extra work is avoided by using Eisenstat's trick [9].
In [1], the implantation process was simulated in the following way: n impurity atoms were placed in a region Vtot which was 8000 times larger than the discrete doping region Vdisc.Only those atoms (k)   that fell within the discrete doping region were retained.In this way, one generates a binomial distribution [10], which in the limit n Na Vtot, p= Vdisc Vtot--+O and np---a (where a=Na I/disc is the mean value of the process), approaches a Poisson distribution k a E, (5) This observation suggests that one can mimic ion implantation by drawing a random number k from the Poisson distribution itself and uniformly distribute these k impurity atoms within the discrete doping region of the device by using triplets of independent uniformly distributed ran- dom numbers.This simplification is essential when dealing with devices with wider gates and very high substrate doping.

RESULTS AND CONCLUSIONS
The test device used for our convergence studies has 0.1 lan gate length and 0.2 lam gate width.The nonuniform finite-difference grid used in these sirnulations has 55 40 55 grid points along the x (channel length), y (depth) and z (width) axis.The predetermined threshold for self-consistency for the electrostatic potential is 10 -5 thermal voltages.
In Figure 2, we compare the convergence properties of various Poisson's equation solvers for Vs VD VBS 0 V and VG V. From the results shown, it follows that preconditioned Bi- CGSTAB method is more economical than the ILU and SIP methods.It gains very little when using more than two internal loops.The use of the ILU preconditioner reduced the number of itera- tions by approximately a factor of two.The slope of these curves is another indication of the performance of each solver.The convergence behavior of our continuity equation solvers is shown in Figure 3.The applied bias used in these simulations is Vs VBs=0V, VD 10 mV and VG V. Again the precondi- tioned Bi-CGSTAB method shows superior con- vergence behavior.Gummel's decoupling scheme converged in 21 (13) cycles when using Bi-CGS- TAB (ILU) method.However, the computational work within each Gummel cycle for the ILU method is more than 200 iterations, so that the total simulation time for achieving self-consistency using this method is significantly longer.
We also studied the width dependence of the threshold voltage Vth in scaled Si MOSFETs.
Here, we use 0.1 lm gate-length and 0.05, 0.1, 0.15 and 0.2 lm gate-widths.Statistical averaging of the results obtained over 25 different devices was performed.The doping profile of each device is evaluated using our simplified procedure.Using the bisection method, we determine the value of the gate voltage VG Vth for which the drain current equals ID=Ioer 10-1A.A set of 20 transfer characteristics for devices with Wg=0.2 lm which have different number and distribution of impurity atoms under the gate is shown in Figure 4.The continuum doping model results are also shown in this figure.For the devices with discrete doping distribution under the gate, the subthreshold slope does not vary significantly from device to device, and is 78.06+1.75mV/decade.Even though the gate-width of these devices is 4 times larger than that used in [1], we still observe a significant spread of the transfer  In the inset we show the width dependence of the threshold voltage for subthreshold conduction.
characteristics along the gate axis due to the nonuniformity of the potential barrier which allows for early turn-on of some parts of the channel.This observation is supported by the fact that the averaged threshold voltage is always shifted about 10 mV to lower voltages.As shown in the inset of Figure 4, this threshold voltage shift is essentially independent of the width of the gate.
In conclusion, for solving very large systems of algebraic equations derived from 3D problems, iterative methods remain the only viable solution technique.The greatly improved convergence rates achieved by preconditioned conjugate gradient methods means that the solution of such problems is more attainable now than ever before.For both 3D Poisson and 3D electron current continuity equations, the Bi-CGSTAB method showed rela- tively smooth and superior convergence behavior when compared to the ILU methods.
FIGURESchematic view of the device structure.The box underneath the gate shows the spatial extension of the discrete doping region.

FIGURE 2
FIGURE 2 Convergence properties of the 3D Poisson equation solvers.
FIGURE 3 Convergence properties of the 3D electron current continuity equation solvers.

FIGURE 4
FIGURE 4 MOSFETs subthreshold transfer characteristics.In the inset we show the width dependence of the threshold voltage for subthreshold conduction.