Investigation of numerical algorithms in semiconductor device simulation

Investigation of numerical algorithms in semiconductor device simulation

Solid-Srare Electronics Vol. 30, No. 8, pp. 813-820, 1987 Printed in Great Britain. All rights reserved 0038-I101/87$3.00+ 0.00 Copyright 0 1987Perga...

767KB Sizes 18 Downloads 157 Views

Solid-Srare Electronics Vol. 30, No. 8, pp. 813-820, 1987 Printed in Great Britain. All rights reserved

0038-I101/87$3.00+ 0.00 Copyright 0 1987Pergamon Journals Ltd

INVESTIGATION OF NUMERICAL ALGORITHMS IN SEMICONDUCTOR DEVICE SIMULATION AKIRA YOSHII, MASAAKI TOMIZAWA and KIYOWKI YOKOYAMA NlT, Atsugi Electrical Communications Laboratories, Atsugi-shi, Kanagawa, 243-01, Japan (Received 21 February 1986; in revised form 4 December 1986) Abstract-The

algorithms used in semiconductor device simulation are investigated. Inversion algorithms, such as SOR, SI, generalized ICCG and Crout methods are compared in terms of convergency and required computer resources for various devices and bias conditions. For linearization of the basic equations, a quasi-coupled method is also compared with Gummel’s conventional decoupled method. Numerical experimentation shows that even the SOR method, which has t’ie slowest convergency among these algorithms, efficiently provides good results when used properly. Moreover, the quasi-coupled method is also effective in linearizing the basic equations for transient analysis or high bias conditions without a significant increase in the required memory. Consequently, a properly used simulator having several algorithms is shown to be necessary for two- and three-dimensional analysis. Furthermore, guidelines for applying these numerical algorithms effectively are described in detail.

1. INTRODUCTION

With the progress in VLSI fabrication technology, devices with higher speed and integration performance and smaller sizes have become possible. By using Si, megabit RAMS and logic LSIs with 50k gates have been developed and much effort is continuing in this direction. Traditionally, devices with specific functions have been developed by repeating the processing-testing cycle several times. With the advance of increasingly complex integrated circuits, however, this empirical approach for device and process design has become ineffective, because of the excessive time and expense. An alternative approach to this problem utiIizes sophisticated numerical simulations[ l-31 in process and device development, and it has proven to be both reliable and cost effective. Therefore, effective process and/or device simulators have now become feasible. Numerical simulation in devices was begun by Gummel[4] who introduced a 1-D finite difference method in 1964. One- and two-dimensional programs based on Gummel’s algorithm soon followed[S-81. The rapid progress in computer capabilities and improvements in numerical techniques have made it possible to use 2-D simulators for device design. For example, the introduction of the Strongly Implicit (SI)[9] and generalized Incomplete Cholesky Conjugate Gradient (ICCG)[l&lZ] for matrix inversion has reduced computational time in some cases. It should be. noted that the ICCG method mentioned here can be applied to both symmetric and asymmetric matrices. Supercomputers, such as the CRAY, offer an order-of-magnitude speed increase in floating operation rates and have now become available in the field of device and process simulation. As a result, 2-D simulators[ 13-181 can be effectively used in analysis for short-channel effects in FETs and lateral-

current flows in bipolar transistors. Moreover, 3-D simulators[ 19-231 have also been developed to model narrow-channel effects in FETs and current spreading due to the extrinsic base region in bipolar devices. In addition, transient analyses have also been reported[ 19,24-261. Under these circumstances, device simulators would seem to be nearing their full-growth stage. However, some problems remain in the practical use of 2- and 3-D device simulators. Many algorithms have been used in this field which solve large matrices resulting from discretization of basic equations, such as the Poisson’s and current continuity equations. Unfortunately, a unified discussion on the convergency and effectiveness of these algorithms has so far not been realized. For example, it has not yet been determined that the ICCG or SI method is more effective than SOR for any device and any condition. Since general device simulators should be applicable to various kinds of devices, the most proper method should be chosen for analyzing devices and bias conditions. Therefore, a 2-D simulator has been developed which includes many inversion algorithms for large matrices, such as SOR, SI, generalized ICCG and the Crout method. Bipolar transistors, MESFETs and MOSFETs are analyzed using this simulator and the convergency and effectiveness of the algorithms are compared under various conditions. Gummel’s successive iteration is also compared with the quasi-coupled method in which basic equations are linearized using Newton’s method. In this method, the resultant Jacobian matrix is iteratively solved by block. In the next section, the proposed device simulator is briefly outlined. Section 3 is devoted to a discussion on the effectiveness of the algorithms in device simulation.

813

814

AKIRA YOSHII et al.

2. BASIC EQUATIONS AND NUMERICAL APPROACHES

combination velocity, is ignored for simplicity as pointed out by Cottrell et a1.[31]. Along the side The electrical properties of a semiconductor device boundaries, homogeneous Neumann boundary concan be specified by the equations: ditions are assumed which insure that no current can flow in or out of the device, At the Si and SiOZ V2$=-_r(P-.+N,t-N,), (1) interface in a MOSFET, the boundary conditions, where normal displacement vector component are (2) continuous on both sides, are satisfied automatically by the usual box integration method[7]. Nonlinear partial differential equations (8HlO) -q(R -G). (3) are linearized by two methods. First, using Gummel’s successive iteration algorithm, which is also called the Auxiliary equations are: decoupled method, the equations are decoupled such t = -q~,,nV$ + qD,Vn, (4) that each one can be regarded as an independent linear equation for each iteration cycle. The linearJ’ = -q&p V’i+ - qD,Vp, (5) ized differential equations are discretized by the n = nie ev[q (IL - cpJW1, (6) box integration method and in a 2-D analysis, a five-diagonal band matrix is usually obtained for P = 4 ew[q (cp, - $YW. (7) each equation. In a transient analysis, time is Here, I(/ is the electrostatic potential, n and p are the discretized by the backward Euler scheme to obtain the stiffness of the calculation. To effectively solve electron and hole densities, rp,,and ‘pPare the quasithe resulting matrix, many inversion algorithms are Fermi potentials for an electron and hole, and G and prepared such as the SOR, SI, ICCG and Crout R are the generation and recombination rates. In methods. In the procedure, both (n,p) and addition, J:, and J, are the electron and hole current [exp( - cp,), exp(cp,)] are used as variables. Under high densities, p,, and pP are electron and hole mobilities, applied-voltage conditions, the preceding variables D, and DP are electron and hole diffusion coefficients, were selected to avoid exponential overflow. An n, is the effective intrisic concentration, q is the electric charge, k is Boltzmann’s constant, and T is additional interesting method, in which the quasiFermi levels, cp, and ‘pP are directly used, was not the lattice temperature. adopted here, because it requires more complicated When Boltzmann statistics and the Einstein formulation than the above two. relationship are applied and the equations are nonUsing Newton’s method, equations (8HlO) are dimensionalized in the usual manner, the following expressions can be obtained as the basic equations for linearized simultaneously[31,32]. With this linearization, the following set of linear equations is the simulator: derived as: F, = V2@,- G2 exp(@,) + @3exp( -@,)

-qg+vJ.=q(R-G),

q$+vjp=

+N,+-N,=O, F2 =

Vbc, exp(@,)V@,] - R + G - t

F,=VbPexp(-@,)V@,]--R+G-z=O,

(8) = 0,

(9) (10)

where @,, Qi, and @, represent $, exp( -cp,) and exp(cp,), respectively. In the case of the physical models, Shockley-Read-Hall[27] and some avalanche models are used for generation-recombination. Furthermore, Scharfetter-Gummel’s empirical formula[28] and the modified formula[29] are adopted for the mobility model. The band-gap narrowing effect using the Slotboom model[30] is also taken into account. Two types of boundary conditions are introduced. For ohmic and Schottky contacts, Dirichlet boundary conditions are assumed, where the potential values are fixed to correspond with the applied voltages. The carrier densities are also fixed at the equilibrium values. For Schottky contacts, however, the effect of majority carrier density modulation at the contact resulting from the finite thermionic re-

E;

$;

Z]

[;;;I=[:;]

(11)

where F, is defined as aFJ&P,. This is called a block matrix. The Jacobian matrix F, which is constructed by the block matrices E;,, is represented in Fig. 1.

Fig. 1. Form of matrix F derived from linearization with Newton’s method and discretization with the usual box method.

815

Numerical algorithms in device simulation Since F is too large directly, an iterative explained below. The basis for this mations are made for &PI is obtained from

and complicated to be solved algorithm is adopted which is algorithm is as follows: estiS@, and &&, and a solution for the first row equation:

F,,b@, = -F,

- (F,,6@, + F,,6@,).

(12)

From the second and third row of matrix equation (1 1), && and 6@, are calculated in the same manner according to: F,,6@, = --I$ - (F2,6@++ F&@J,

(13)

F&D, = - F3 - (F,,&D, + F&ZJ~).

(14)

and

The cycle is repeated until the convergence criterion is satisfied. This solution is called a quasicoupled method in contrast with a fully coupled method. In this procedure, one parameter is introduced to accelerate convergency in the following manner: &P+’ =&J-k+ W,(lWf’

-&P).

(1) e,=max(lb (2) e,=]lb

-A&]),

-&]]/fi,

(3)

ek=max((xk-xk-ll/lxkl),

t4)

ek=maXtbk-Xk-li/(l

+Ixkl%

These definitions are classified into two categories. One includes Definitions (1) and (2), where the maximum error is estimated by the residual. Definition (1) means the maximum residual in the k th iteration and Definition (2) shows the mean square error in the iteration. The other consisting of (3) and (4) corresponds to the variation rate estimation of the solution vector. These maximum errors in the SOR iteration for one Gummel loop are indicated in Fig. 2. The calculation is accomplished for an npn transistor with V,, = 0.5 V, where the extrapolated values from the solutions of VBE0.4 and 0.45 V are used as initial values in the SOR iteration. It can be seen from the figure that Definitions (1) and (2) provide more severe

(15)

Here, a&‘+’ is the new increment vector with which the (k + 1)th solution of @ is defined as 4’+’ = &k + 66kf’. Next, 66” is the same vector in the kth iteration, 6ak + ’ is the solution vector obtained from eqns (12Hl4) in the (k + 1)th iteration, and wB is the acceleration parameter which varies between 1.0 and 1.9. This acceleration parameter is introduced on the analogy of SOR method and is determined experimentally. Since in this quasi-coupled method, three equations are linearized simultaneously, it is expected to show quadratic convergence as reported by Cottrel et a1.[31]. At each Gummel and quasi-Newton iteration, terminal currents are calculated from the obtained n, p and $ values. Current densities which are obtained from eqns (4) and (5) are integrated along the closed loop which contains the contacts. In this calculation, Gauss’s theorem is used to change the surface integral into a line integral, Next, calculated terminal currents are compared with those obtained from the previous iteration. Unless the change in currents between the two iterations is smaller than the definite value, the iteration is repeated. Moreover, the iteration is also repeated until the terminal currents satisfy Kirchhoff’s conservation law. Therefore, the accurate numerical solution of eqns (8HlO) is determined by evaluating the convergency of variables @,, Gjzand G3 and of the terminal currents which should satisfy the conservation law. 3. RESULTS

maximum error e, in the k th iteration are considered:

NUMBER OF ITERATIONS

4

0

8

& l&l

40 80 120 NUMBER OF ITERATIONS

IOO-,/’-----___ -

.

160

-____--_-_- (I)

HOLE .---------_________

(2)

AND DlSCUSSION

3.1. Error dejinitions in relaxation methods To solve the matrix equation Ax = b with a relaxation method (A is N x N matrix and x and b are N-dimensional vectors), the following definitions for

0

120 40 80 NUMBER OF ITERATIONS

160

Fig. 2. Numerical errors in SOR as a function of number of

iterations.

816

AURA YOSHII et al.

convergence criteria than (3) and (4) for the hole current equation. In relaxation calculation for a matrix inversion, the calculated result does not necessarily mean the accurate solution, even if e, by Definition (3) or (4) becomes nearly zero. It is a necessary condition but is not a suffcient condition for convergency. Therefore, in some case, these definitions are insufficient for estimation of convergency, although these are used generally. On the other hand, when e, by Definition (1) or (2) becomes less than a very small value, it is guaranteed that the results obtained are a good ,approximation of the exact solution. Hereafter, the definitions which belong to the former category, are used as the maximum error in the iteration loop.

0

20

30

NUMBER OF OUTER ITERATIONS

5

3.2. Convergency of decoupled method Bipolar devices are modeled by Poisson’s and two current continuity equations corresponding to majority and minority carriers. Especially, the accurate calculation is required for base current, because it plays an important role in device operation and is usually l/l0 and l/1000 of the collector current. In the case of an npn transistor, the base current mainly consists of holes. Therefore, we focus on the discussion of convergency of holes in bipolar transistor analysis. To investigate cost effective methods, matrix inversion techniques, such as SOR, SI and ICCG, are studied using the simulator based on Gummel’s algorithm. This simulator usually has double loops. One is the outer loop where attributes to the linearization of eqns (8HlO), and the other is the inner loop where the resulting matrix from equation discretization is solved by iterative methods. Therefore, there are two kinds of numerical errors corresponding to the loops. The error in the outer loop is defined as the maximum variation rate of the solution vector in linearization step. As the error of the inner loop, the normalized form of the definition (1) is adopted. As a total measure of convergency of the nonlinear partial differential equations, the base current is also evaluated. Figure 3 shows the calculated results for a typical npn transistor with low bias condition (V,, = 0.4 V). Here, a nonuniform rectangular mesh, which consists of 43 grids in x-direction and 53 grids in y-direction, is used. In every iterative method, the inner calculation for matrix relaxation is repeated same times (in this case 200). In Figure 3(a), the error of the outer loop in the Crout method, which is not an iterative method but a direct one, is plotted as a reference. From Figure 3(a), it can be seen that ICCG has the best convergency among these iterative methods with SI next. Especially since the convergency of ICCG is nearly equal to that of the Crout method. It is also clear from Fig. 3(c) that the final residuals for ICCG are much smaller than those for SOR and SI. Although the same tendency is observed with respect to convergency of the base current, the difference among

IO

3

0

30

NUMBER 0:OUlER

5-i

a 5

Il&ATIONS

lo-?

0

IO

20

30

NUMBER OF OUTER ITERATIONS Fig. 3. Convergency in bipolar transistor analysis for a low bias condition, (a) error of the outer loop, (b) base current and (c) maximum residual against repetition number of the outer loop.

them is not so distinguished as the inner and outer errors. Figure 4 shows the results for high bias condition ( l’BE= 0.9 V). In each iterative method, the inner calculation is also repeated 200 times. From Figures 4(a) and (b), it can be observed that there is no difference in the convergency of the outer loops among the three methods and the calculated base currents are also less affected by the methods. Moreover, the convergency of the base currents approximately corresponds to that of the outer iterations. However, the final results which are evaluated by the base currents, are not seriously influenced by the accuracy of the inner calculations. Since the coupling between Poisson’s and two current continuity equations becomes strong for high current regions, calculation for the outer iteration becomes dominant and accurate solution of each matrix in a linearization step has not such an important meaning for solving the coupled equations. To investigate dependency of the inner loop rep-

817

Numerical algorithms in device simulation

lv,fosv

-

SOR

-m-e

ICCG

0 NUMBER ‘& OUTER :iERATIONP

0.d

I

0

lo

20

30

NUMBER OF OUTER ITERATIONS 0

200

100 NUMBER

OFINNER rlERATIONS

Fig. 5. Variation of total iteration number of iterative methods and bias conditions.

the inner iteration

0

IO

20

30

NUMBER OF OUTER ITERATIONS Fig. 4. Convergency in bipolar transistor analysis for a high bias condition. Each figure corresponds to each one in Fig. 3.

etition number on the base current in detail, the number is varied from 200 to 10. In Figure 5, the total iteration numbers for iterative methods and two bias conditions are plotted against the inner loop number. The total iteration number is defined as the product of the inner and outer loops to obtain the definite base current. In this case, the value which is represented as 4 digits of significance figures of the base

current calculated at the 100th outer iteration is used. For the low bias condition, as the number of the inner loop is decreased, the total iteration number does not decrease. This means that the solution with high accuracy in matrix inversion is needed for obtaining the exact base current effectively. Moreover, the difference among iterative methods are clearly observed and a highly convergent method like ICCG, has less iterations than others. On the other hand, for the high bias condition, the reduction of the inner loop becomes effective because the total iteration numbers reduce approximately in proportion to decrease of the inner loop except for ICCG. ICCG with

less than 100 does not converge within 100 outer iterations. It seems that ICCG, especially for an asymmetric matrix, has slow convergency at the beginning of iteration steps, although it shows high convergency totally. Table 1 shows the ratio of the calculation time per iteration for the iterative methods. It is normalized by the time of SOR. The total calculation cost for iterative methods can be obtained from the Fig. 5 and Table 1. From these results, it can be concluded that for low bias region where the base current is extremely small, a highly convergent iterative method with many inner calculations is preferred. For the high bias region, however, simple algorithms with less inner loops, such as SOR and SI, become sufficient for obtaining the base current with the same accuracy as those by ICCG and Crout methods. Therefore, SOR and SI methods become more effective in that region than ICCG and Crout methods in computation time. For biases below V,, = 0.3 V, the currents calculated are not so exact even by the ICCG and Crout methods using carrier densities n and p as variables, If the result values were exact, the sum of the emitter, collector and base currents should be zero from eqns

Table 1. Ratio of computation time per iteration for iterative method

Computation per iteration

SOR

Si

ICCG

1.0

2.1

4.4

time

AKIRA YOSHII et

818

al.

Table 2. Calculated currents for a bipolar transistor at low bias voltages compared for variable selection and numerical precision. I’., = 0.0 V

current Carrier Densities V BE 0.1 v

Variables Precision IE

IB 0.2 v

fE IB IC I&B + Ic)

0.3 v

IL

Quad 1.921lOOE-11 4.318264E-12 1,482914E-I 1 I .0001872 6.6296608-10 4.1330528-l 1 6.216185E-10 1.0000256

Quasi-Fermi Double

Double 1.744113E-11 - I .432890E-09 4.556207E-11

-

6.614050E-10 I .444266E-09 6.524116E-IO

I .7660OOE-11

4.0933218-12 I .949686E- I I 0.74861655

6.610825E-10 4.1079558-l 1 6.265427E-10 0.98914367

2.647753E-08 -1.551190E-09 2.6 15044E-08

LiUa + Ic)

2.6479538-08 3,60474lE-10 2,611970E-08 0.99997567

2.6479078-08 3.605113R-10 2.612384E-08 0.99980058

0.4 v

IL IB IC IdUB + c-1

1.1226llE-06 4.7488798-09 I. I 17902E-06 0.99996447

I. 122609B06 5.2030738-09 1.1179338-06 0.99953071

1.122655B06 4.7530528-09 I. 1179068-06 0.99999639

0.5 v

L L 1, IElVB+ Id

4.876312E.05 I .340640E-07 4.863101E-05 0.99995993

4.876312E-05 1.342676E-07 4.863104E-05 0.99995514

4.876520E-05 I .341896E-07 4.863102E-05 0.99999980

IB IC

(2) and (3). This inaccuracy in numerical calculation for low bias conditions can be attributed to a fractional drop in the numerical value which is represented by finite digits in the computer. The terminal currents calculated using the exponential of quasiFermi potentials as variables, satisfy the conservation law better than those using the np method. This is because the domain of variables exp(-cp,) and exp(cp,) are smaller than that of variables n and p in the low bias region. When we use the quad-precision expression for variables, the small fractional drop in numerical values is suppressed and more accurate results are acquired even in the low bias region as shown in Table 2. As an example of FET’s, a MESFET is analyzed here. In this analysis, computation time to get the same current-voltage characteristics is compared between SOR and ICCG methods. The result is shown in Fig. 6. It can be seen that the computation time by SOR is a few times shorter than that by ICCG. This is explained as follows. In FET, the majority carrier is dominant and the minority carrier can be neglected in its normal operation condition. Moreover, the convergency of the majority carrier does not seriously influence on the calculation of the current. Therefore, SOR which is simple, is more effective than ICCG method.

moderate bias regions where log I, has linear relation with VBE,Gummel’s method is more economical than the quasi-coupled method, because its calculation is simpler. As VBEincreases, the quasi-coupled method excels in terms of computation time. For bias above 0.9 V, the computation time by the quasi-coupled method reduces to l/2 as much as that by the decoupled one. Under high bias conditions, a large amount of minority carriers are injected so that the coupling between the basic equations (lH3) becomes strong. High injection of minority carriers can be seen from the device characteristics where the curve representing the relation V,,- log I, differs from the

6.0 Z Z?l

v, =o.ov

3.3. Convergency of quasi-coupled method Device characteristics for an npn transistor are simulated by both Gummel’s decoupled method and the quasi-coupled method. The CPU time ratio between two methods is shown in Fig. 7. The collector current (I,) against the base-emitter bias voltage (V,,) is also depicted in the figure. In low and

0.0

I 0

I

1.0

0.5 V,

(VI

Fig. 6. The ratio of CPU time required by SOR and ICCG in MESFET analysis.

819

Numerical algorithms in device simulation

method becomes efficient. Based on these considerations, a mixed simulation, in which the two algorithms are used according to the conditions, is very effective for reducing computation time.

4.

0

-0.2

-0.4 -0.6 -0.0 -1.0 -1.2” v,,

(VI

Fig. 7. The CPU time between decoupled Gummel and quasi-coupled methods in bipolar transistor analysis.

straight line. The base-emitter voltages where computation by the quasi-coupled method becomes more efficient than that by Gummel’s method, agree well with high injection conditions. Another example of the transient analysis for a bipolar transistor at two VBE conditions is presented in Figs 8(a) and (b). When a voltage pulse is applied to the base electrode, the transistor state suddenly changes. In this case, the potential, electron, and hole distributions vary significantly and large currents flows. Therefore, the quasi-coupled algorithm is required for transient simulation. As time progresses, the variations of n, p and I,+gradually decrease until a steady state condition reached. When the state of the transistor approaches steady state, the decoupled

01 0

20

60 80 40 TIME lpsl

100

CONCLUSIONS

A two-dimensional simulator has been developed in which many large matrix inversion algorithms, such as SOR, SI, ICCG, and Crout, are included. Under various bias conditions, these algorithms are estimated for bipolar transistors, and FETs. The numerical experiments show that even the SOR or Sl methods efficiently provides good results with proper use, although they have slower convergency than ICCG and Crout methods. For example, in high current regions for bipolar transistors, they yield accurate base currents with half or less of the computational effort needed in the Crout method due to its simplicity. Under ordinary conditions in FETs, similar tendencies with respect to computational time are observed. In current regions which are too small and high convergency is required for carriers, the Crout or ICCG method prevails. Moreover, in these regions, it is effective to select exponentials of quasi-Fermi levels as variables instead of carrier densities n and p. The quasi-coupled method is studied through comparison with Gummel’s decoupled method. For high bias conditions of an bipolar transistor, it shows higher convergency than Gummel’s method. The efficiency between two methods can be related with device characteristics successfully.

“Ol---------’20

40 TIME

60

60

100

(psi

Fig. 8. The ratio of CPU time between decoupled and quasi-coupled methods in transient analysis (a) V,,= 0.6 + 0.8 V(lOps), (b) V,, = 0.7 + 0.9 V(lOps).

AKIRA YOSHII et al.

820

Acknowledgements-The authors would like to thank Drs Hisakazu Mukai, Tsuneta Sudo and Takayoshi Nakashima for their continuous encouragement and guidance throughout this study.

REFERENCES

5. 6. I. 8. 9. 10. 11. 12. 13. 14.

W. L. Engl, H. K. Dirks and B. Meierzhagen, Proc. IEEE 71, 10 (1983). R. W. Dutton, IEEE Trans. Electron Dev. ED-30. 968 (1983). W. Fichtner, D. J. Rose and R. E. Bank, IEEE Trans. Electron Dev. ED-SO, 1018 (1983). H. K. Gummel, IEEE Trans. Electron Deu. ED-11,455 (1964). A. DeMari, Solid-St. Electron. 11, 33 (1968). A. DeMari, Solid-St. Electron. 11, 1021 (1968). J. W. Slotboom. IEEE Trans. Electron Dev. ED-20.669 (1973). G. Manck, H. H. Heimeier and W. L. Engl, IEEE Trans. Electron Dev. ED-21. 403 (1974). H. L. Stone, SIAM J. Namer. Anal. 5; 530 (1968). J. A. Meijerink and H. A. Van der Vorst, Mafh. Comput. 31, 148 (1977). F. Y. Chang and L. F. Wagner, Electron. Letl. 18, 658 (1982). T. Wada and R. L. M. Dang, Electron. L&t. 18, 265 (1982). T. Adachi, A. Yoshii and T. Sudo, IEEE Trans. Electron Deu. ED-26, 1026 (1979). M. Tomizawa, H. Kitazawa, A. Yoshii, S. Horiguchi and T. Sudo, IEEE Trans. Electron Dev. ED-28, 1148 (1981).

15. S. P. Gauer and D. H. Navon, IEEE Trans. Electron Dev. ED-23, 50 (1976). 16. T. Toyabe and S. Asai, IEEE J. Solid-St. Circ. SC-14, 375 (1979). ._ T. Wada and K. Tanieuchi. IEEE Trans. I /. K. Kato. Electron ‘Dev. ED-32, 458 (1985). y 18. N. Kotani and S. Kawazu, Solid-St. Electron. 11, 681 (1981). and 19. E. M. Buturia, P. E. Cottrell, B. M. Grossman K. A. Salsbura. IBM. J. Res. Develop. 25, 218 (1981). M. Tomizawa, S. Horiguchi 20. A. Yoshii, H.-Kitazawa, and T. Sudo, IEEE Trans. Electron Dev. ED-29, 184 (1982). A. Yoshii and T. Sudo, IEEE 21. R. Kasai, K. Yokoyama, Trans. Electron Deu. ED-29, 870 (1982). 22. N. Shigyo and R. Dang, IEEE Trans. Electron Dev. ED-32, 441 (1985). 23. T. Toyabe, H. Masuda, Y. Aoki, H. Shukuri and T. Hagiwara, IEEE Trans. Electron Deu. ED-32, 2038 (1985). 24. M. S. Mock, Solid-St. Electron, 24, 959 (1981). II Conf., 25. W. L. Engl and H. Dirks, Proc. NASECODE p. 34 (1981). Solid-St. Electron. 26, 907 (1983). 26. K. Yamaguchi, 27. W. Shockley and W. T. Read, Phys. Rev. 87,835 (1952). and H. K. Gummel, IEEE Trans. 28. D. L. Scharfetter Electron Deu. ED-16, 64 (1969). IEEE Trans. Electron Dev. ED-26, 1068 29. K. Yamaguchi, (1974). 30. J. W. Slotboom, Solid-St. Electron 20, 279 (1977). I 31. P. E. Cottrell and E. M. Buturla, Proc. NASECODE Conf., 31 (1979). 32. R. E. Bank, D. J. Rose and W. Fichtner, IEEE Trans. Elecfron Deu. ED-SO, 1031 (1983).