Research Article Current Issue Versions 5 Vol 3 (4) : 20030403 2020
On the History of the Numerical Methods Solving the Drift Diffusion Model
: 2020 - 10 - 08
: 2020 - 11 - 02
: 2020 - 12 - 30
120 15 0
Abstract & Keywords
Abstract: In 1964 Hermann Gummel published the first numerical solution method for the one-dimensional Drift Diffusion model. In his seminal paper [1] already the nonlinear iteration method and the basics of the discretization method named after him are outlined. Soon after this paper appeared many research groups worldwide tried to solve the Drift Diffusion equations in two and more dimensions applying predominantly general finite element discretization methods which were very popular at these days. Due to this a large variety of different codes solving the multidimensional Drift Diffusion equations based on many different space discretization schemes existed in the seventies. However already in the nineties all Drift Diffusion simulators being of importance for semiconductor device design in industry and academia still used Gummel’s nonlinear iteration method but were entirely based on just one specialized space discretization method, which incorporates the basic ideas of the Scharfetter-Gummel discretization scheme [2]. All other codes which were not based on this special space discretization method had nearly vanished already in the nineties and this is still the case today. This paper tries to shed some light on the hidden reasons for this astonishing development.
1.   Introduction
In the following the development of the numerical solution methods for the Drift Diffusion (DD) equations will be reviewed, since much can be learned from this development for comparable tasks in the future. Aspects covered well in the literature will be only shortly mentioned by citing the appropriate references. Other aspects that are very important as well but rarely mentioned in literature or even nearly forgotten today will be highlighted. Especially the aspect of preserving the inherent "stability" of the underlying differential equations in terms of monotonicity properties during the discretization and solution process will be carefully examined. The general flow of arguments presented here follows closely those outlined in Reference[3]. Several details have been previously published as well in References [4, 5]. Some hints concerning the Hydrodynamic (HD) model [6] will be given at the end.
2.   History of the Numerical Models
Since the quasi Fermi potentials Φn and Φp (imrefs) for electrons or holes, respectively, cover a much smaller range of numerical values compared to the electron or hole densities n, p, intermediate and final solutions of the drift diffusion system of equations are typically saved by storing the electrostatic potential Ψ, and the imrefs Φn and Φp instead of Ψ, n and p. Moreover, formulating the DD equations in Ψ, Φn, and Φp makes it much easier to introduce the solution algorithms that are typically applied. The stationary drift diffusion equations for a homogeneous semiconductor (e.g., silicon) device formulated in these variables have the following form (see References [7, 4, 8, 9] for details and a derivation),
where e is the elementary charge and κ is the permittivity of the different materials for which Poisson’s Equation (1) is solved. Moreover ni is the intrinsic density, ND , NA are the ionized donor and acceptor concentrations and µn , µp the electron and hole mobilities of the semiconductor material with homogeneous band gap within which the electron and hole continuity Equations (3, 2) are solved. In addition, VT is the thermal voltage and G the generation density within the semiconductor. For simplicity it is assumed that at all contacts the Dirichlet boundary conditions of the ideal Ohmic contact model (see References [4, 9] for details) are valid for all three potentials Ψ, Φn and Φp and that at all other boundaries homogeneous Neuman type boundary conditions can be applied.
It can be shown that the above system of differential equations has a unique solution provided a number of reasonable assumptions is fulfilled especially for the generation term G. Please refer to Reference [10] for a general theory and to References [11, 12] for special results concerning the DD set of equations.
There is an important property, that deserves special attention.
I. All above operators TP - TE are in divergence form and result in important conservation laws if integrated over a finite volume and after the application of the divergence theorem.
For example integrating TH+TE over the simulation domain results in Kirchhoff’s law for the stationary terminal currents of the device under consideration. If a numerical device model is used inside a circuit simulator it is absolutely mandatory that this law is exactly reproduced by the numerical model. Therefore, it is very important to maintain the validity of such conservation laws in some sense during the discretization process.
The most important solution algorithm for the DD system is an iterative method often addressed as Gummel’s nonlinear relaxation method [1]. Assuming the result Ψk, Φp,k, Φn,k, after iteration k as known, this method evaluates the new approximate solution after iteration k+1 by solving the three boundary value problems (1) - (3) successively as follows:
For the above partiell differential equations it is always assumed that the variable with the highest number of derivatives in the equation is updated and the other variables are kept unchanged. The individual nonlinear equations in (4) are typically solved by Newton’s method, which converges very fast and robust even for bad initial solutions, if the underlying equation is nearly linear. Therefore, Newton’s method is typically applied for the operators and defined in (2) and (3) and not for the operators TH and TE, since and are nearly linear in the new variables
provided carrier generation G has no dominant influence. These new variables are often addressed as Slotboom variables because they were first introduced in Reference [13] and the advantages they have for generating “stable” discretization schemes were probably mentioned in Reference [7] for the first time. Nevertheless it is still possible to calculate the Newton updates based on the numerically more convenient Jacobians of TH and TE and the variables Φp and Φn. The only modification necessary for performing the Newton iterations for the more linear operators is to modify the Newton update itself as shown below for the hole continuity equation and the solution function before (Φp,b) and after (Φp,a) one Newton step.
is the Newton update calculated using the Jacobian of TH for the variable . For a general report on using alternative solution variables for enhancing convergence please refer to Reference [14].
In order to understand which criteria in addition to the mandatory consistency criterion should be considered for the discretization of the continuous operators TP - TE it is very useful to look at the Jacobians of , and that are necessary for performing Gummel’s nonlinear relaxation method based on Newton’s method. These Jacobians are typically evaluated for some existing intermediate solution Ψb, Φp,b, Φn,b and operate on the functions , , . If in addition for the calculation of the Jacobians the dependence of the mobilities on the solution variables is neglected,
these Jacobians have the following form.:
All above partial derivatives with respect to solution variables are assumed to be Frechet derivatives on suitable function spaces [15]. If only direct recombination and Shockley-Read-Hall (SRH) recombination [9] is considered for the Frechet derivative of the carrier generation G
holds. With this additional condition all Jacobians (7)-(9) have important properties that again deserve special attention. They are
II. self adjoint,
III. positive definite,
IV. of monotone typ,
if appropriate boundary conditions are assumed for the function spaces considered [15]. Especially property IV is very important, since it means that for all these Jacobians monotonicity theorems hold (Reference [15], Chapter 23.5) which restrict for example the possible form of the update functions , , during the Newton iterations required for Gummel’s nonlinear relaxation method very much and enhance the robustness and convergence properties of this solution method decisively.
The mathematical properties of the model after discretization should be as similar as possible to the properties of the continuous model! If this is fulfilled the discrete model is an analogon of the continuous model, even if only coarse grids can be afforded, which is typically the case. This property is very important for the discretization error control on coarse grids. Therefore discretization methods are preferred which are able to conserve conditions I-IV in some discrete sense.
The first 2D DD simulations used a rectangular solution domain and tensor product grids [16, 7] so that standard finite difference discretization methods for tensor product grids could be applied, that conserved most of the conditions mentioned above. However many device cross sections were not rectangular, so that at least at the beginning of multi-dimensional numerical semiconductor device modeling many groups developed discretization methods (see References [17, 18, 19, 20, 5, 21, 22] and citations therein) for non rectangular solution domains based on the finite element approach [23], but it turned out that general finite element methods typically have problems to conserve conditions I-IV [24, 18, 5, 25]. Finally, even for solution domains with complicated polygonal boundaries which cannot be discretized efficiently by a tensor product grid, the method of choice used today by the vast majority of DD and HD simulators is a straight forward generalization of the integration method published in Reference [26], Chapter 6 for 2D tensor product grids. This method is able to conserve conditions I-IV as will be shown below. This generalization for 2D problems was already briefly mentioned in Reference [26] using the early work of Reference [27] as guideline. This general method is addressed as box integration method in Reference [4], as box method in Reference [5] and today often named finite volume method [28]. An especially well documented example of this historical development from the application of general finite element methods to the final exclusive application of the finite volume method is the development of numerical device modeling codes at IBM research. There two groups independently developed two general numerical device modeling codes. One group with a clear focus on the application of general finite element methods [24, 19, 20] and the other group shifting more and more from hybrid finite element/ finite volume discretization methods to the exclusive application of the finite volume method [29, 18]. Both codes were developed over a decade until the beginning of the eighties. Ten years later the finite element code development is not mentioned at all any more in a comprehensive review paper about the TCAD development at IBM with more than 200 citations [30].
Possibly the first mathematical analysis of the general box integration method for the DD model in 3D is due to Reference [31]. The box integration method can be interpreted as a finite element method [5], but the grid elements for this method (e.g. triangles) cannot be considered as the basis from which the discretization method proceeds like it is the case for a general finite element method [23] but instead the elements are constructed in a unique way based on the predefined grid points. The method is based on the construction of two dual grids the Voronoi diagram and the Delaunay tessellation. The early work of the two Russian mathematicions M.G. Voronoi [32] and B.N. Delone [33] is typically cited in this context but the method in 2D is even much older and has been rediscovered various times. Figure 1 below is used to explain the basic principles of this method. The n grid points Pk (n = 9 in this example) with the coordinates are considered as given. In a first step for each point with index k its Voronoi volume Vk is constructed as the set of all points in space that are closer to Pk than to any other grid point.
Moreover for each grid point Pk its environment Sk is defined by
M indicates the 1D Lebesgue measure M1 for 2D problems and the 2D Lebesgue measure M2 for 3D problems. With this environment the boundary of Vk is given by

Figure 1.   Voronoi diagram (dashed lines) and Delaunay tessellation (thick drawn lines) for a 2D example with 9 grid points.
The union of all Voronoi volume boundaries defines the Voronoi diagram. Based on this diagram the dual grid structure of the Delaunay grid is defined by defining the set of edges D for this grid by the straight lines li ,k between the points Pi and Pk for which and i and k vary between 1 and n with . Please note that the li ,k are oriented curves starting at Pi and ending at Pk. If only the edges of the Delaunay grid without orientation are important li ,k and lk ,i can be considered as equal. The typical elements in such a Delaunay grid for 2D problems are triangles and tetrahedra for 3D problems. But 2D and 3D rectangular tensor product grids are special cases for this method and fit perfectly into this framework. The DD and HD problems are typically formulated as boundary value problems on some finite region and in order to incorporate Newman type boundary conditions into the framework of box integration in a natural manner the boundary is typically assumed to be composed of edges (faces in 3D) of the delaunay grid elements. In the given example in Figure 1, is the closed polygonal line composed of ,,,,, and and is the interior of this closed polygonal line. The Voronoi volumes for the grid points on are typically unbounded but for box integration bounded boxes are mandatory. This leads to the following definitions:
Clearly and is only relevant if is part of the Delaunay grid. If no additional constraints are fulfilled it can happen that for some edge of the Delaunay grid. This is typically not good for the consistency of the discretization method and mostly avoided by constructing the grid in such a manner that for example for 2D problems like the given example the Delaunay triangles that have some common boundary edge with have all interior angles smaller than 90 degrees. Such triangles are called acute or nonobtuse. See for example Reference [34] for algorithms generating grids in such a manner. If this additional condition is fulfilled for all edges of the Delaunay grid holds and only the boxes of grid points on the boundary share some boundary with like for the given example. Similar additional conditions with comparable consequences are considered as well for 3D problems [31]. For the box integration process presented here it is not necessary that all triangles are nonobtuse for preserving condition IV but for finite element discretization schemes this condition must be fulfilled [34, 5].
The boundary value problems TP, and defined in Equations (1) - (3) and considered as individual problems that are solved separately have the following common form:
u (r) is the solution variable. Therefore u is either or or . Morover holds always. The dependence of a and f on and u considers the typical physical models for the mobilities and the generation rate and their dependence on the solution variables [9]. The first step always performed for the box method is to integrate Equation (16) over the box for each grid point , which is not determined by a Dirichlet boundary condition. Moreover the divergence theorem is used to transform the suitable parts of the integral over into integrals over . This yields:
is a unit vector having the same orientation as . The above formulation is for the 3D case, where the first integrals integrate fluxes over an area and the remaining part is a volume integral. In the 2D case the first integrals are line integrals and the second part is an area integral. If it can typically be assumed that for this part of the boundary a Newman type boundary condition holds such that the flux through this part of the boundary is zero. This applies for instance for in the example, if for no Dirichlet boundary condition is given. It is clear that
where the first integral is related to and the box integration over , whereas the second integral is related to the point and the box integration over . During the discretization process the integrals in Equation (17) are typically approximated independently by difference approximations and quadrature rules such that a consistent discrete approximation is generated. For conserving property I it is important that Equation (18) holds exactly after discretization. If this is the case the discretized formulas for two points and with (like and in the example) can be summed leading to an expression where the sum of the discretized integrals over and is represented by the discretized flux integrals over the boundary of which is . The resulting expression can be interpreted as a discrete version of the divergence theorem for and its boundary. If the discretized problem is solved exactly this discretized version of the integral theorem holds exactly and does not depend on any discretization error. Relations of this kind are very helpful for checking the global numerical accuracy and the consistent calculation of for instance terminal currents.
In order to study how to preserve properties II-IV during the discretization process equation (16) is simplified again by neglecting the dependence of a on the solution variable u and considering only the dependence of direct and generation on the solution variable. Thus (16) simplifies to
and holds always. If the discretization preserves property I the discretized integral over Fk,j in Equation (17) can be considered as a function . This implies that G depends as well on everything clearly connected to like Pk and Pj . Since Equation (18) shall be preserved must hold. Lets assume that the discrete solution is represented by a vector u with n entries uj and each uj is the discrete approximation of the function u at the point Pj . The Jacobian of the equation system after discretization should be self adjoint (symmetric). This requires the condition . There are not too many alternatives left if these above two conditions must be fulfilled simultaneously. One discretization formula, for which both conditions hold, is
Here is a suitable mean value of a(r) that can be considered as a function of and for which must be satisfied. If in addition the integral over the box in Equation (17) is discretized using the simplest quadrature formula
the strict diagonal dominance of the Jacobian matrix of the discrete system is guaranteed as well. indicates the 2D Lebesgue measure for 2D problems and the 3D Lebesgue measure for 3D problems. So far only the discretization for all points which are not given by a Dirichlet boundary condition has been studied. The set of indices of these points should be given by SB . Moreover SD contains all indices of points that are given by a Dirichlet condition. The Jacobian entries for the latter points are simply 1 for the main diagonal and 0 for the other entries. In summery the Jacobian of the discretized boundary value problem (19) is strictly diagonal dominant, all main diagonal elements are strictly positive and all other elements are negative or zero. Such matrices are positive definit and M-matrices as well, which means that their inverse matrix has only elements that are positive or zero [26]. These properties are very beneficial for a large number of solution algorithms solving linear equations involving matrices. The convergence of iterative methods like the Jacobi or Gauss-Seidel methods is guaranteed [26], semi-iterative methods like conjugated gradient algorithms work well [35] and even Gaussian elimination profits because a pivot element search is not necessary and the accumulation of the rounding error during elimination is well controlled. Finally and probably most important the discrete system introduced above is of monotone type, which means that important stability inequalities can even be derived for the maximum norm, which is the most important norm for practical applications. For example in Reference [36], Chapter 4.2.2, the following is proven. If for all kSB
and for all , then the following stability inequality holds for the discrete solution:
This inequality is directly applicable for the estimation of the maximum Newton correction if the Jacobian on the left hand side of Equation (22) is used to solve the discretized nonlinear problem (19) by Newton’s method. Moreover such stability inequalities are very useful for evaluating upper bounds for the discretization error even for the nonlinear problem (19). In Reference [25] an excellent example is given demonstrating clearly how bad the discretization error is controlled if the off diagonal elements of the Jacobian of the discretized electron continuity equation have both signs. Moreover it is shown as well in Reference [25] for the same problem that the discrete solution gets much more accurate and very well controlled by the applied voltages if all off diagonal elements of the Jacobian become always negative or zero after a modification of the grid. The underlying reason is not as falsely stated in Reference [21] that the original grid had one obtuse triangle but that the original grid was not the Delaunay grid constructed on the basis of the Voronoi diagram, whereas the modified grid is the Delaunay grid. As pointed out already earlier, a Delaunay grid may contain obtuse triangles!
Another advantage of the discretization scheme presented above is that it allows a straight forward incorporation of the Scharfetter-Gummel discretization formula for the balance equations, which was originally developed for rectangular grids [2, 4, 5]. The application of this discretization formula is mandatory for achieving accurate simulation results on coarse grids. For general finite element schemes it is typically very difficult to incorporate this formula but for the scheme described above it is easily done by choosing as follows:
For holes:
For electrons:
are appropriate mean values of the mobilities that can be regarded as a function of the edge . Results concerning the consistency and convergence of the discretization scheme described above for the balance equations can be found in Reference [31].
Gummel’s nonlinear relaxation method (4) performed in such a manner that the relevant discrete Jacobian matrices fulfill conditions II-IV, which typically means that derivatives with respect to the solution variables are neglected for the mobilities and impact ionization, converges nearly always even for very bad initial solutions. One of the rare counter examples is given in Reference [14]. The convergence is typically slow for high current applications but is very predictable so that the solution accuracy can be estimated very reliably during the iteration process [37]. This allows to switch to more coupled methods like a simultaneous Newton method not before the solution accuracy is so high that the simultaneous Newton method is in the range where it converges quadratically. Of course in this case all derivatives should be considered in the Jacobian matrix of the fully simultaneous Newton method. The availability of this combination of solution methods for the DD model featuring high robustness even for bad initial solution and high accuracy at the same time is possibly one reason why the DD model is still the numerical device model that is applied by far most even for nanoscale devices, where its physical accuracy is certainly questionable [38]. The above comments concerning the beneficial effect on robustness of neglecting the derivatives of the mobilities and impact ionization apply as well to other nonlinear relaxation methods [39, 14] and are even valid for the fully coupled Newton method outside the range of quadratic convergence. It is rather straight forward to extend the box integration method including a Scharfetter-Gummel type discretization for the energy flux densities to the HD model [6]. A nonlinear relaxation method with comparable convergence properties to Gummel’s method has been published in Reference [40] and evaluated in Reference [41] for the HD model. For this method and the numerical solution algorithms of the HD model in general convergence robustness increases as well decisively if certain derivatives with respect to mobilities and impact ionization are turned off in order to enhance the “stability” of the discretized equations.
3.   Conclusion
Based on the historic development of the space discretization and solution methods for the Drift Diffusion model it is shown how important it is for the error control on coarse grids to preserve especially the monotonicity properties of the underlying partial differential equations in the final discretized model. The author believes that this observation should serve as a guideline for the development of discretization methods for transport models in future.
[1] H. K. Gummel: A self-consistent iterative scheme for one-dimensional steady state transistor calculations, IEEE Transactions on Electron Devices, 455–465 (October 1964).
[2] D. L. Scharfetter, H. K. Gummel: Large-Signal Analysis of a Silicon Read Diode Oscillator, IEEE Transactions on Electron Devices 16(1), 64–77 (1969)
[3] B. Meinerzhagen: Effiziente Rechenstrategien zur Bauelementsimulation, Doctoral Thesis (RWTH-Aachen, 1985)
[4] W. L. Engl, H. K. Dirks, B. Meinerzhagen: Device Modeling, Proc. of the IEEE 71, 10–33 (1983)
[5] R. E. Bank, D. J. Rose, W. Fichtner: Numerical Methods for Semiconductor Device Simulation, IEEE Trans. Electron Devices 30(9), 1031–1041 (1983)
[6] B. Meinerzhagen, W. L. Engl: The Influence of the Thermal Equilibrium Approximation on the Accuracy of Classical Two-Dimensional Numerical Modeling of Silicon Submicrometer MOS Transistors, IEEE Trans. Electron Devices 35(5), 689–697 (1988)
[7] J. W. Slotboom: Computer-Aided Two-Dimensional Analysis of Bipolar Transistors, IEEE Trans. Electron Devices 20, 669–679 (1973)
[8] C. Jungemann, B. Meinerzhagen: Hierarchical Device Simulation: The Monte-Carlo Perspective, Computational Microelectronics, ed. by S. Selberherr (Springer, Wien, New York 2003)
[9] S. Selberherr: Analysis and Simulation of Semiconductor Devices (Springer, Wien 1984)
[10] O. A. Ladyshenskaya, N. N. Uraltseva: Linear and Quasi-linear Elliptic Equations (Academic Press, New York 1968)
[11] M. S. Mock: Analysis of Mathematical Models of Semiconductor Devices (Boole Press, Dublin 1983)
[12] T. I. Seidmann: Steady State Solution of Diffusion Reaction Systems with Electrostatic Convection, Nonlinear Analysis Theory, Methods & Applications 4, 623–637 (1979)
[13] J. W. Slotboom: Iterative scheme for 1 and 2-dimensional DC transistor simulations, Electron. Lett. 5, 677–678 (1969)
[14] K. H. Bach, H. K. Dirks, B. Meinerzhagen, W. L. Engl: A New Nonlinear Relaxation Scheme for Solving Semiconductor Device Equations, IEEE Trans. Computer–Aided Des. 10(9), 1175–1186 (1991)
[15] L. Collartz: Funktionalanalysis und numerische Mathematik (Springer-Verlag, Berlin Heidelberg New-York 1968)
[16] H. H. Heimeier: A Two-Dimensional Numerical Analysis of a Silicon N-P-N Transistor, IEEE Trans. Electron Devices 20, 708–714 (1973)
[17] J. J. Barnes, R. J. Lomax: Finite-Element Methods in Semiconductor Device Simulation, IEEE Trans. Electron Devices 24, 1082–1089 (1977)
[18] E. Buturla, P. E. Cottrell, B. Grossmann, K. Salsburg: Finite-Element Analysis of Semiconductor Devices: The Fielday Program, IBM J. Res. Develop. 25, 218–231 (1981)
[19] G. D. Hachtel, M. H. Mack, R. O‘Brian, B. Speelpenning: Semiconductor Analysis using Finite- Elements - Part I: Computational Aspects, IBM J. Res. Develop. 25, 232–245 (1981)
[20] G. D. Hachtel, M. H. Mack, R. O‘Brian: Semiconductor Analysis using Finite-Elements - Part I: IGFET and BJT Case Studies, IBM J. Res. Develop. 25, 246–260 (1981)
[21] J. F. Burgler, R. E. Bank, W. Fichtner, R. K. Smith: A new Discretization Scheme for the Semiconductor Current Continuity Equation, IEEE Trans. on Computer-Aided Design 8(5), 479–489 (1989)
[22] P. J. Mole: Discretization of the semiconductor current continuity equation for finite element solvers in 2 and 3 dimensions. In: Proc. Fouth Conf. Numerical Analysis of Semiconductor Devices and Integrated Circuits-NASECODE IV, ed. by J. J. H. Miller (Boole Press, Trinity College, Dublin, Irland 1985) pp. 429–435
[23] W. G. Strang, G. J. Fix: An analysis of the finite element method (Prentice-Hall, Englewood Cliffs 1973)
[24] G. Hachtel, M. Mack, R. O‘Brien: Semiconductor Device Analysis Via Finite Elements. In: Proc. of the Eighth Asilomar Conference on Circuits and Systems (1974) pp. 332–338
[25] C. S. Rafferty, M. R. Pinto, R. W. Dutton: Iterative methods in semiconductor device simulation, IEEE Transactions on Electron Devices 32(10), 2018–2027 (1985)
[26] R. S. Varga: Matrix Iterative Analysis, Series in Automatic Computation (Prentice-Hall, Englewood Cliffs, New Jersey 1962)
[27] R. H. MacNeal: An Asymmetrical Finite Difference Network, Quart. Apll. Math. 11, 295–310 (1953)
[28] C. Großmann, H. G. Roos: Numerische Behandlung partieller Differentialgleichungen (B. G. Teubner Verlag / GWV Fachverlage GmbH, Wiesbaden 2005)
[29] P. Cottrell, E. Buturla: Steady State Analysis of Field Effect Transistors via the Finite Element Method. In: IEDM Tech. Digest (IEEE, 1975) pp. 51–54
[30] E. B. D. C. Cole, S. Furkay, K. Varahramyan, J. Slinkman, J. Mandelman, D. Fotyi, O. Bula, A. Strong, J. Park, T. L. Jr, J. Johnson, M. Fischetti, S. Laux, P. Cottrell, H. Lustig, F. Pileggi, D. Katcoff: The Use of Simulation in Semiconductor Technology Development, Solid State Electronics 33, 591–623 (1990)
[31] M. Sever: Delaunay Partitioning In Three Dimensions and Semiconductor Models, Compel 5, 75–93 (1986)
[32] G. Voronoi: Nouvelles application des parametres continus a la theorie des forme quadratiques. Deuxieme memoire. Recherche sur les parallelloedres primitifs, Journal für reine und angewandte Mathematik 134, 198–287 (1908)
[33] B. N. Delone: Sur la sphère vide, Bulletin of the Academy of Sciences of the U.S.S.R. 7, 793–800 (1934)
[34] B. S. Baker, E. Grosse, C. S. Rafferty: Nonobtuse triangulation of polygons, Discrete & Computational Geometry 3, 147–168 (1988)
[35] J. A. Meijerink, H. A. vanderVorst: An iterative solution method for linear systems of which the coefficient matrix is a symmetric M-matrix, Math. Comp. 31, 148–162 (1977)
[36] A. A. Samarskij: Theorie der Differenzenverfahren (Akademische Verlagsgesellschaft Geest & Portig K.-G., Leibzig 1984)
[37] K. Bach: Nichtlineare Relaxationsverfahren zur Lösung von Differentialgleichungssystemen am Beispiel der Halbleitergleichungen, Doctoral Thesis (RWTH-Aachen, 1991)
[38] C. Jungemann, T. Grasser, B. Neinhuis, B. Meinerzhagen: Failure of Moment-Based Transport Models in Nanoscale Devices near Equilibrium, IEEE Trans. Electron Devices 52(11), 2404– 2408 (2005)
[39] B. Meinerzhagen, H. K. Dirks, W. L. Engl: Quasi-simultaneous solution method: A new highly efficient strategy for numerical MOST simulations, IEEE Trans. Computer–Aided Des. 4, 575– 582 (1985)
[40] B. Meinerzhagen, K. Bach, I. Borg, W. L. Engl: A New Highly Efficient Nonlinear Relaxation Scheme for Hydrodynamic MOS Simulations. In: NUPAD IV Technical Digest (Seattle 1992) pp. 91–96
[41] Y. Apanovich, E. Lyumskis, B. Polsky, P. Blakey: An Investigation of Coupled and Decoupled Iterative Algorithms for Energy Balance Calculations. In: Simulation of Semiconductors and Processes, ed. by S. Selberherr, H. Stippel, E. Strasser (Springer Verlag, Wien 1993) pp. 233–236
Article and author information
Bernd Meinerzhagen
Bernd Meinerzhagen is Full Professor and Head of the Institute for Electron Devices and Circuits at the Technical University Braunschweig, Braunschweig, Germany. He received the Dipl.-Ing. degree in electrical engineering, the Dipl.-Math. degree in mathematics, the Dr.-Ing. degree in electrical engineering, and the “venia legendi,” all from the RWTH Aachen (Aachen University of Technology), Aachen, Germany, in 1977, 1981, 1985, and 1995, respectively. From 1978 to 1986, he worked mainly on the development of numerical device modeling codes as a Research and Teaching Assistant at the RWTH Aachen. In 1986, he joint AT&T Bell Laboratories, Allentown, PA, as a Member of Technical Staff, where he developed advanced numerical models for MOS substrate and gate currents. From 1988 to 1995 he was the Head of the Research and Development Group for Silicon technology modeling and simulation (TCAD) at the RWTH Aachen and he continued this research as Professor at the University of Bremen, Bremen, Germany between 1995 and 2003. His research interests are focused on the physics, characterization, modeling and design of Si/SiGe integrated devices and circuits and on the mathematical foundations of electrical engineering. Prof. Meinerzhagen has coauthored a book and more than 200 papers published in international journals and conference proceedings and has been a Technical Program Committee member of IEDM, ESSDERC, SISPAD, IWCE and other conferences.
Publication records
Published: Dec. 30, 2020 (Versions5
Journal of Microelectronic Manufacturing