Computer Physics Communications 65 (1991) 133-136 North-Holland
133
Semiconductor device simulation Karl Gustafson Department of Mathematics, University of Colorado, Boulder, CO 80309-0426 USA
A short account of certain interesting problems in semiconductor physics, process, and device modelling is given. This is a branch of microelectronics which requires the most advanced computational physics and the largest supercomputers for its understanding
1. Introduction
At the IMACS 1st International Conference on Computational Physics in Boulder, CO, 11-15 June, 1990, Professor Karl Hess, Director of the N S F Center for Computational Electronics, Urbana, Illinois, gave a fascinating plenary lecture entitled " M o n t e Carlo simulation of electron transport in semiconductors at high energies", (K. Hess).
Because all other plenary presentations at the conference are represented in this volume I decided to give here, for completeness, a short account of these presentations. This is a nonexpert account, based upon very limited notes and recollections from the conference, but also supported by some background knowledge and experience about semiconductor device simulation problems which I have gained over the last two years in some research.
(1.1)
We, the organizers, were grateful for this presentation, and we understand all too well the impossibility of finding enough time to provide new original written manuscripts for all conferences to which one is invited to speak, especially in fields as active and important as semiconductor electronics. At this same IMACS conference, Dr. R. Kent Smith of the A T & T Bell Laboratories, Murray Hill, New Jersey, organized a small Special Session on Semiconductor Device Simulation. The three excellent lectures in this session were "Computational aspects of semiconductor device simulations" (R.K. Smith),
(1.2)
"Iterative methods for semiconductor device simulations" (W. Coughran),
(1.3)
"Quantum device simulations" (U. Ravaioli).
(1.4)
2. Short summaries of the lectures and references
An excellent reference for lecture (1.1) is Hess [1], in a recent issue of Physics Today devoted to ultrafast devices. The need to use Monte Carlo methods comes from the use of the Boltzmann equation when modelling electron transport in devices for which statistical effects are included but quantum effects are neglected. The supercomputer-driven demand to produce smaller and smaller chips has driven device size down to the submicrometer range. At those small sizes the relatively large supply voltages induce very large electric fields and very high electron energies. Within such devices this energy may be equivalent to temperatures above 10 000 kelvin. At such extremes it is not hard to imagine lots of nonlinear effects becoming important, to say nothing of discrete electron-hole dynamics and scattering effects.
0010-4655/91/$03.50 © 1991 - Elsevier Science Publishers B.V. (North-Holland)
134
K. Gustafson / Semiconductor device simulation
For more information, at this point, I direct the reader to the recent article by Hess [1] and the references given there, including his recent book [2]. A good reference for lecture (1.2) is Bank et al. [3]. In his lecture, Dr. Smith broke SCDS into three parts: 1. device simulation, 2. circuit simulation, 3. process simulation. Device simulation involves coupled partial differential equations (the so-called drift-diffusion equations, see section 3 below) for which the case of one space dimension may be regarded as understood but for which the cases of two and three space dimensions still pose many difficulties. Circuit simulation may involve as many a s 106 components. Process simulation, e.g. the details of impurity diffusion, oxidation, is in the infancy stage. The general goal in all cases is to model devices as small as possible. Historically, the 1970's successfully treated MOSFETS, the 1980's saw increased attention to edge effects, and the 1990's will attempt three-dimensional simulations. There is a need to better match the drift-diffusion equation results with the Monte Carlo results. The latter are also much more expensive to run. A recent reference for lecture (1.3) is Bank et al. [4], which carries exactly the same title as lecture (1.3). Dr. Coughran concentrated in his lecture on the two-dimensional steady drift-diffusion equations. The associated numerics involve the following considerations. First one needs to decide between triangular and quadrilateral grids, depending on the geometry of the device being modelled. Then there are the still unresolved discretization questions about how to best generalize the one-dimensional Scharfetter-Gummel ( S - G ) scheme to two and three dimensions. To solve the resulting nonlinear systems one cannot seem to avoid coupled Newton methods. In the associated linear solvers, iterative methods, especially those of conjugate gradient type, appear to be best. Among those being used are GSSOR, ILU(k), M I L U ( k ) ABF, G M R E S ( k ) , O R T H O M I N ( k ) , BICG, and CGS. Dr. Ravaioli (1.4) then discussed the need to
bring in quantum effects. Devices under investigation now may have layers only 5 atoms thick. At pn junctions the current/voltage characteristics will allow tunneling. The full system of equations, expanded beyond the drift-diffusion equations, now includes coupled Poisson, Boltzmann, and Schrodinger equations, boundary value problems along with eigenvalue problems. Scattering effects and quantum interfaces become important. For references see Sols et al [5] and Hess [1,2].
3. The drift-diffusion equations The basic set of partial differential equations modelling device simulation were introduced by Van Roosbroeck in 1950 [6]. These govern the electrostatic potential +(x, t), the electron concentration n( x, t ), the hole concentration p ( x, t ), and the electron and hole currents Jn(x, t) and Jp(x, t) within a device, according to: potential
- ~ A ~ = q( p - n + C ) ,
electron continuity
an at = ql div Jn -- R,
hole continuity
~p _ Ot
1 div Jp - R , q
electron current
Jn = ql~ n ( UT W n - n ~7~b) ,
hole current
Jp = - q ~ p ( Ur~Tp + p~Tq~),
(3.1) where c = permittivity of device material, q = electron charge, UT= thermal voltage, are constants, and where C = doping profile is regarded as given. Generally, the recombination-generation rate R and the electron mobility /~o and hole mobility ~p depend is a specified nonlinear manner on the principal unknowns ~, n and p. For a derivation of these equations and for further details see Markowich [7]. Putting the electron current and hole current relation into the electron continuity and hole continuity equations in (3.1) yields three basic coupled nonlinear second-order drift-diffusion equations for the unknowns qJ(x, t), n ( x , t) and p ( x , t). It turns out that the electron and hole concentration solutions have an enormous dynamic
K. Gustafson / Semiconductor device simulation
range, e.g. 1 0 20 , with very steep gradients in certain parts of the semiconductor. Part of this is due to intentionally steep gradients in the imposed doping profiles. To handle such steep and long gradients, a rescaling is often employed. Assuming Boltzmann's statistics we may write the carrier concentrations as n = e (+-O)/UT,
p : e (w-~)/UT.
(3.2)
From this change of variables and assuming a stationary solution we arrive at the following three coupled nonlinear elliptic equations in the unknowns tk, v, w A ~p =
e '~- v -
e w - '~ -
C ( x ) ,
V" (/t, e~-V~Tv) = R,
eW-
135
niques from physics to device engineering. One simulates the microscopic motion of thousands of individual carriers first, and then one averages to get the carrier densities. This is much more expensive, computationally, than the continuum driftdiffusion equations. The recent D A M O C L E S program of IBM can simulate Monte Carlo 3-D Boltzmann transport coupled to 2-D Poisson device cross sections. A description of earlier versions of PISCES and M I N I M O S can be found in Pinto et al. [8] and Selberherr et al. [9]. An account of the D A M O CLES program is given in Laux et al. [10].
5. Equations and algorithms (3.3)
w) = - R .
Here we have abused notation by absorbing UT in ~, v and w, absorbing this and q in ¢, with similar renormings absorbed in the mobilities ~n and /tp and in R. We have also taken the liberty of absorbing throughout the intrinsic concentration /'/intrinsic, which is for example about 101°/cm 3 for silicon at room temperature. The formulation (3.3) reduces the dynamic range at the expense of increasing the nonlinearity.
4. Software Because of the industrial importance of these problems, considerable software has been assembled over the last twenty years. Here in the United States, the Stanford group developed the 2-D simulator PISCES, which is widely used for the modelling of both Si and GaAs. The Vienna group under the direction of S. Selberherr has developed a stationary 3-D code called M I N I M O S 5. The M I N I M O S codes have been widely used in Europe. I donot know which are the most prevalent codes in Japan. In shrinking device sizes, the reduced validity of preaveraged mobilities and diffusion coefficients has led to a transfer of Monte Carlo tech-
Although much work has been done, relatively few mathematics, as compared to electrical engineers and physicists, have become seriously engaged in these problems. For such nonlinear partial differential equations, which may be expected to possess a number of nonunique solutions in some parameters ranges, and none in others, better understandings of the boundary conditions and initial conditions may be important for the selection of the physical from the nonphysical numerical solutions. There are also interesting analogies with combustion equations and certain fluid-flow algorithms in which similar considerations arise [11]. It is important to distinguish the time-dependent from the stationary cases. Algorithms from the latter cannot generally be expected to be very optimal for the former. Recent reductions of false diffusion by using improved Box schemes are restricted thus far to the 2-D stationary equations. Little is understood about the stability of the 2-D transient algorithms. Major simplifications are employed in recent applications of the method of characteristics approach to the time-dependent equations. Mixed finite elements still have problems in the representation of and error bounds for the exponential terms in the equations. There is the bothersome obtuse triangle problem. There are difficulties in preserving the conservation laws in adaptive gridding, and in the higher-dimensional upwinding essential to most simulation al-
136
K. Gustafson / Semiconductor deoice simulation
gorithms. Thus, in spite of much excellent work already done, there are a number of important mathematical and algorithmic considerations yet to be resolved, especially in the transient cases for 2 and 3 space dimensions.
[4] [5] [6] [7]
References [1] K. Hess, Physics Today 43 (February 1990) 34. [2] K. Hess, Advanced Theory of Semiconductor Devices, (Prentice-Hall, Englewood Cliffs, N J, 1988). [3] R. Bank, W. Coughran, W. Fichtner, E. Grosse, D. Rose
[8] [9] [10] [11]
and R.K. Smith~ IEEE Trans. Electron Devices 32 (1985) 1992. R. Bank, W. Coughran, M. Driscoll, W. Fichtner and R.K. Smith, Comput. Phys. Commun. 53 (1989) 201. F. Sols, M. Macucci, U. Ravaioli and K. Hess, Appl. Phys. Lett, 54 (4) (1989) 350. W. Van Roosbroeck, Bell System Tech. J. 29 (1950) 560. P. Markowich, The Stationary Semiconductor Device Equations (Springer, Berlin, 1986). M. Pinto et al., IEDM Tech. Digest (December 1984) 288. S. Selberherr et al., IEEE Trans. Electron. Devices 27 (1980) 1540. S. Laux, M. Fischetti and D. Frank, IBM J. Res. Develop. 34 (July 1990) 466. K. Gustafson, unpublished notes.