A TWO-DIMENSIONAL FINITE ELEMENT ANALYSIS OF THE STATIONARY SEMICONDUCTOR DEVICE EQUATIONS by PATRICK PABLO C H A V E Z B. A . Sc., The University of British Columbia, 1995 A THESIS SUBMITTED IN PARTIAL FULFILMENT OF THE REQUIREMENTS FOR THE DEGREE OF MASTER OF APPLIED SCIENCE in THE F A C U L T Y OF G R A D U A T E STUDIES Department of Electrical and Computer Engineering . We accept this thesis as conforming to the required standard THE UNIVERSITY OF BRITISH COLUMBIA July 1997 © Patrick Pablo Chavez, 1997 In presenting this thesis in partial fulfilment of the requirements for an advanced degree at the University of British Columbia, I agree that the Library shall make it freely available for reference and study. I further agree that permission for extensive copying," of this thesis for scholarly purposes may be granted by the head of my department or by his or her representatives. It is understood that copying or publication of this thesis for financial gain shall not be allowed without my written permission. Department of f /<f c Jn*C«. I' a^J Co r^^^je r £^/ ^ The University of British Columbia Vancouver, Canada Date A y „ < f / r ; Iff?-DE-6 (2/88) Abstract The ability to model the steady-state field inside active structures, such as a transistor, is an important aspect of monolithic microwave integrated circuit (MMIC) design. This paper focuses on such an active zone of semiconductor material, and presents a finite element analysis of the classical semiconductor equations. The semiconductor equations are very nonlinear and govern the potential and carrier density distributions in semiconductor materials. A previously developed finite element method (FEM) formulation of these equations, referred to as the current conservation model, is re-introduced and re-derived with compact matrix notation. It is shown how this formulation can be solved with the Newton-Raphson iterative scheme. Then, a newly developed F E M formulation, referred to as the advection-diffusion model, of the continuity equations is described in detail. It is shown by example how this formulation solved with Gummel's iterative technique is very numerically robust. These two different solution methods of the steady-state system of coupled Poisson and continuity equations are combined into a final solution algorithm that exploits their strengths. As a specific example, GaAs MESFETs are the focus of implementation, and the resulting potential field and carrier density distributions are used, to calculate various MESFET parameters such as electrode currents, voltage gain, capacitances, and conductances. Finally, various extensions to the F E M approach involving the application of the method of moments (MoM) are briefly discussed and partially demonstrated. These extensions are intended to compensate for the assumptions and simplifications, mainly with respect to the artificial boundary conditions, used in the original stand-alone F E M formulations. i i Table of Contents Abstract ii List of Figures v 1 Introduction 1 1.1 The Big Picture 1 1.2 Objectives and Motivation 5 1.3 Thesis Structure 6 2 Problem Specifications 8 3 FEM Formulation 12 3.1 Elements 12 3.2 Discretization 15 3.2.1 Poisson Equation 16 3.2.2 Continuity Equations 17 3.2.2.1 Current Conservation Model 17 3.2.2.2 Advection-Diffusion Model 19 4 Solution Methodologies 22 4.1 Introduction to Iterative Techniques 22 4.2 Gummel's Method 25 4.3 Newton-Raphson Method 31 5 Application 37 5.1 Mesh Spacing 37 5.2 Simulations 38 5.2.1 Example 1 41 5.2.2 Example2 47 iii 6 Extensions to FEM Model 57 6.1 Method of Moments and Integral Equations 57 6.2 Parasitic Capacitances and Green's Integral Equation 58 6.2.1 Complex Images 59 6.2.2 Multipipe Model 60 6.2.3 Results 60 6.3 Artificial Boundary Condition Adjustments 62 6.3.1 Green's Integral Equation 62 6.3.2 Boundary Integral Equation 66 7 Conclusion 68 7.1 Accomplishments 68 7.2 Future Work 69 Bibliography 71 Appendix A Derivation of Shape Functions 77 Appendix B Capacitance Calculations for a Multi-Port System 80 i v List of Figures Figure 1 Layout of a monolithic microwave integrated circuit [1] 1 Figure 2 2D cross section of a MESFET device. 10 Figure 3 Six adjacent quadratic triangular elements 12 Figure 4 Local (left) and global (right) shape functions for Node i 13 Figure 5 Direct iterations; Newton-Raphson iterations 23 Figure 6 Direct iterations; Newton-Raphson iterations 24 Figure 7 Direct iterations; Newton-Raphson iterations 25 Figure 8 (a) Schottky barrier strip and (b) results for N ^ = 5 x 1 0 ^ cm"^. . . 28 Figure 9 Electron drift velocity versus electric field intensity 41 Figure 10 MESFET device for Example 1. 42 Figure 11 Mesh grid for Example 1 -959 nodes, 448 elements 42 Figure 12 Potential field 44 Figure 13 Electron carrier distribution 44 Figure 14 Current .density field 45 Figure 15 Lj vs. for Example 1. . . 46 Figure 16 Lj vs. Vg ( V d = 2.5 V) for Example 1 . 4 6 Figure 17 Mesh grid for Example 2 - 1453 nodes, 688 elements. . . . . 48 Figure 18 Electron carrier concentration and potential distributions 49-Figure 19 Current density field 51 Figure 20 1^ vs. for Example 1. 52 Figure 21 (a) Drain and (b) gate current versus gate voltage 53 Figure 22 Electrode charges for varying potentials on (a) the drain, (b) the gate, and (c) the source 54 Figure 23 Small signal model of MESFET in Example 2. 56 v Figure 24 MESFET with strip conductor feeds 57 Figure 25 Modified small signal model of MESFET in Example 2 61 Figure 26 Gain versus frequency plots for circuit models. 62 Figure 27 Modified mesh grid for Example 2 - 1677 nodes, 800 elements. . . 63 Figure 28 Drain current versus drain-source voltage 63 Figure 29 Artificial boundary adjustment results. 65 Figure B l Circuit model for two strip conductors. 81 v i 1 1 Introduction 1.1 The Big Picture As with other electrical technology, microwave technology has been following a path of smaller size and increased complexity. Monolithic microwave integrated circuits consist of a planar substrate on which passive and/or active circuit elements are constructed. Passive elements such as microstrip transmission lines, inductors, capacitors, and resistors take the form of various deposited layers of metal, dielectric, and resistive films on the surface of the substrate. In addition to this, a substrate made of semiconductor material may have active devices grown directly into it through ion implantation. For example, gallium arsenide substrates, which are the most common of the semiconductor substrates, allow for the construction of GaAs FETs which have many applications such as low-noise amplifiers, high-gain amplifiers, broadband amplifiers, mixers, oscillators, phase shifters, and switches [1]. Figure 1 shows a simple layout of a monolithic microwave integrated circuit (MMIC). Air Microstrip MIM bridge input line capacitor Inductor Ground Via GaAs Implanted GaAs plane hole FET resistor substrate metalization Figure 1. Layout of a monolithic microwave integrated circuit [1]. 2 The situation of many circuit elements operating at microwave frequencies and occupying such a compact environment that exists in MMICs leaves little room for error in the design and calls for considerable precision in the fabrication of MMICs. Unlike conventional circuit analysis in which a circuit is simply treated as the sum of its components whose operation is independent of their environments, a more thorough analysis is required. Such an analysis must take into account the electromagnetic effects of mutual field coupling among the components, discontinuities and sharp edges in the structures of the components, biasing circuitry, and packaging. Extensive use of specialized computer-aided design (CAD) software that treats the analysis as an electromagnetics problem is vital. In order to predict the voltage-current characteristics of a microwave circuit, the electromagnetic fields must be modeled. Dynamic and static fields induced by passive elements obey Helmholtz' equation and Laplace's equation, respectively, both of which are linear differential equations [1]. Dynamic and static fields induced by active elements also obey Helmholtz' equation and Poisson's equation (Laplace's equation with a source term), respectively, but are coupled to the variable charge distributions (electrons and holes) which are governed by the continuity equations [2]. The continuity equations are nonlinear. Theoretically, it is possible to deal with the problem of calculating the fields occurring in an MMIC such as the one shown in Figure 1 by simply solving the appropriate equation with specified boundary conditions. Practically, due to the complexity of such a problem and limited computing resources, this is simply not feasible. The complexity resides in the diversity of boundary shapes, types of boundary conditions, uncertainty in material parameters, and nonlinearities in active regions that exist in a typical M M I C . As a result, it becomes necessary to break down the problem into smaller ones. A smaller problem may consist of focusing on a single transmission line, capacitor, inductor, or transistor which occupies just a portion of the circuit. In order not to neglect coupling to the rest of the circuit, "appropriate" artificial boundary conditions surrounding the particular problem domain must be employed. Once a sub-problem is defined, a solution technique must be applied to the governing equation. In some cases analytic solutions are available [3], but in the vast majority of cases, 3 analytic solutions are simply not attainable. In such cases, numerical techniques have offered the only possible route to a solution. The finite element method, finite difference method, and method of moments are the three main classes of such techniques. Each has its own strengths and weaknesses. The method of moments has proven to be ideal for linear problems [4]-[8]. It requires an analytic solution of the linear differential equation at hand for a point source. This solution is known as a Green's function. Based on the equivalence theorem which states that any static or dynamic electromagnetic field in a source-free domain can be fully represented by an equivalent distribution of charge or current sources, respectively, and taking full advantage of the principle of superposition of point source solutions, a Green's integral equation that represents the solution of the linear differential equation in terms of source distributions is set up. Intimately related to the Green's integral equation is the boundary integral equation, which also requires knowledge of the Green's function and gives an expression of the field satisfying the governing equation over the entire domain in terms of the field over a closed contour in or around the domain. The method of moments involves discretizing the integral equation by way of an approximate basis function expansion of the true unknown solution and a minimization of the weighted inner product residuals reducing the problem to the determination of the unknown expansion coefficients via a matrix inversion. Finite element and finite difference methods appear similar to each other on the surface but, in reality, have completely different mathematical bases. They both solve differential equations, and they both require that the solution domain be finite. The fundamental difference between the two methods lies in the fact that finite difference methods approximate the differential equation operator and differential equation domain at hand with a difference operator and a finite number of discrete points spread over the original domain, while finite element methods approximate the solution with an expansion of known basis functions and unknown coefficients over the entire domain [9]-[l 1]. Finite element methods are formulated either directly from the original physical arguments used to derive the differential equations, known as a variational formulation, or based on minimization of weighted inner product residuals, known as a weak formulation. In this sense, 4 the method of moments and the finite element method are closely related. Generally speaking, finite element methods are more flexible than finite difference methods when dealing with inhomogeneous solution domains, irregularly shaped boundaries, and derivative boundary conditions. On the other hand, finite difference methods are more easily implemented in terms of formulation and programming since they are grid dependent, while finite element methods are coordinate-independent with mapping built into their formulation. The limitation of the method of moments is its dependency on the availability of an analytic expression of the Green's function. The Green's function satisfies boundary conditions at infinity and at other boundaries such as multi-layered substrate interfaces and the ground plane. Unfortunately, with added complexity, in the form of irregularly shaped boundaries and lossy materials, an integral equation representation of the solution quickly becomes increasingly difficult to derive [12]. Finite difference and finite element schemes do not suffer from such a limitation. On the other hand, moment method solutions of integral equations are typically more efficient than finite difference or finite element solutions to the differential equation of the same problem, and they naturally handle infinite domains. So, whenever dealing with relatively simple passive structures located in infinite domains, the method of moments has the upper hand. For more complex structures, the finite difference or finite element method must be employed [13], [14]. Also, since the principle of superposition does not apply to nonlinear problems, integral equation solutions do not exist, and the method of moments is not applicable to active structures such as transistors. In contrast, since they do not depend on integral equations, finite difference and finite element methods in conjunction with iterative solution techniques offer possibilities to reach the solution in a nonlinear domain. In referring back to the original problem of modeling an entire MMIC one subregion at a time, it becomes apparent that a choice of numerical solution is beneficial. For subregions containing simple, passive structures such as planar strip-like conductors, the method of moments along with the appropriate integral equation are appropriate. For subregions of higher complexity and/or containing active structures, the finite element method would seem to be the better choice over the finite difference method. The reason for this is that it is possible to exploit the 5 mathematical similarities of both the method of moments and the finite element method in order to directly couple the finite element solution in a finite subregion to moment method solutions in any adjacent infinite or finite subregions. Taking into account subregion interactions overcomes the problem of neglecting the mutual coupling of MMIC components. 1.2 Objectives and Motivation For amplifier and digital logic applications, MMIC layouts include transistors which take the form of doped subregions of the semiconductor substrate. This paper focuses on one such active subregion of the MMIC problem and presents a finite element analysis of the governing classical semiconductor device equations. Two different finite element method (FEM) formulations of the steady-state system of coupled Poisson and continuity equations are described, and two different iterative techniques are used to solve for the potential and carrier density distributions inside the active subregion. As a specific example, GaAs MESFETs are the focus of implementation, and the resulting field distributions are used to calculate various MESFET parameters such as electrode currents, voltage gain, capacitances, and conductances. Additionally, techniques based on the M o M that link the finite element domain to any exterior domains are briefly discussed. To this day, finding a general solution to the tightly coupled set of nonlinear elliptic semiconductor equations has presented quite a challenge to the scientific community. There are many aspects to this challenge. The first and most important one is to be able to ensure robustness of the numerical solution. Since the equations are nonlinear and tightly coupled and since the parameters of semiconductor devices, such as doping density distributions, bias conditions, and geometries, vary widely depending on the device, it is difficult to develop a general solution algorithm that is guaranteed to reach a solution in all cases. Another important aspect is the rate of convergence to the solution. Almost all solution methods implement iterative schemes, and large systems of equations must be solved simultaneously at each iteration. This task demands considerable computing power and speed, and therefore efficiency of the solution algorithm is 6 important. Equally important is the accuracy of the final solution, and this requirement often forces a trade-off with the speed of convergence. Generally speaking, almost all solutions, the vast majority of which implement either finite element or finite difference schemes, involve the discretization of the semiconductor equations over the domain of the device, which results in a grid of points or a mesh spanning and representing the domain. The trade-off occurs in the coarseness of the mesh. That is, a finer mesh implies more equations, which require longer solution times per iteration, but results in more accurate solutions, while a cruder mesh implies fewer equations, which require less time to solve per iteration, but results in less accurate solutions. Applications of finite difference techniques have been studied and fine tuned over at least the past forty years in finding solutions to nonlinear elliptic differential equations and about thirty years in finding solutions to the semiconductor equations in particular [15]-[20]. Barnes and Lomax first introduced an F E M solution in [21], although they could not claim superiority of this solution over the highly refined finite difference schemes. Since their work, relatively little has been published regarding the application of the F E M specifically to the semiconductor equations, [22]-[25]. The purpose of this thesis is to expand on the work of Barnes and Lomax in [22] and [23] by first generalizing their F E M approach with compact F E M notation carried over from Zienkiewicz [9] and Buchanan [10]. Following this, a new F E M steady-state formulation referred to as the advection-diffusion model is described. Advantages and disadvantages of both formulations are pointed out, and a final solution algorithm is presented that adopts the best features of both of them. The final algorithm is used to simulate and analyze GaAs MESFET devices similar to those in [22]. 1.3 Thesis Structure Sections 2-4 give a detailed description of the problem and a derivation of two different routes to a solution. The final solution approach offering the best performance involves a combination of two F E M models, the advection-diffusion model and the current conservation model, that are numerically solved by two different techniques, Gummel's method and the 7 Newton-Raphson method, respectively. Following this in Section 5, implementation considerations, such as stability and accuracy restrictions, of the F E M formulations are addressed. Two simple GaAs MESFET problems are introduced and fully analyzed in order to demonstrate the proposed techniques. Finally, various extensions to the F E M approach involving the application of the method of moments (MoM) are briefly discussed and partially demonstrated in Section 6. These "add-ons" attempt to compensate for the assumptions and simplifications, mainly with respect to the artificial boundary conditions, used in the original stand-alone F E M formulations. Ultimately, the goal of such extensions is to account for the effects of any exterior structures on the interior domain. 8 2 Problem Specifications The equations governing potential and carrier density distributions in the active semiconductor regions of an MMIC are the Poisson and continuity equations written as follows [23]: where V-EV<|> = -p qdp/dt = -V-Jp + qG qdn/dt = V - J n + qG (1) (2) and P = J P = Jn = E = q ( N d + p - n) q(pu.hE - D p V p ) q(nu.eE + D n V n ) -V(|> G = net generation rate N _ = donor minus acceptor ion densities <|> = potential n, p = electron and hole carrier densities J n , J p = electron and hole current densities D n , Dp = electron and hole diffusion coefficients q = electron charge, 1.602 x 10"^ C £ = semiconductor permittivity u-n = electron and hole mobilities t = time. (3) (4) (5) 9 The modeling problem is reduced to a two-dimensional cross section of the region, and steady state conditions are enforced, i.e. it is assumed that no variations of the potential and carrier density in the z coordinate and time dimension exist [19]. Also, the problem is further reduced to that of a majority-carrier device, for which the contribution of minority-carrier transport to the total current is negligible. In the case of GaAs MESFETs, electron transport dominates, and neglecting the hole continuity equation is therefore a valid approximation. Since generation and recombination phenomena depend on the existence of both electrons and holes, as is discussed in most texts on semiconductor theory such as Pulfrey and Tarr [26] and Markowich [2], the generation-recombination term is also omitted in the following analysis of a majority-carrier device. The generation-recombination term introduces a nonlinear coupling between the electron and hole continuity equations, so its omission simplifies the finite element formulation somewhat. Whether or not inclusion of this term would diminish the effectiveness of the F E M approaches described here for simulation of semiconductor devices in general has not been investigated. Regardless, how to incorporate this term into the two formulations described here is mentioned where appropriate. With these simplifications, equations (2) and (3) reduce to 0 = V - J n (6) p = q ( N d - n) (7) for electron majority carrier devices. For hole majority carrier devices, a similar expression exists. For GaAs MESFETs with submicron gate lengths, which are of interest in MMIC technology and which will be the focus of later examples, even small applied bias voltages (~ 1-2 V) generate very large electric fields (~ 10^ V/cm) thus causing electrons to reach then-saturation velocities throughout much of the device. Since the saturation velocity for electrons in GaAs is approximately 10^ cm/s [26], transit times for electrons from one end of a one-micron-long MESFET to the other is on the order of 10"^ s. Such carrier transit times give an indication that operating frequencies of a submicron gate MESFET can reach well within tens of gigahertz, the range of millimeter-wave fields. In fact, it is this high frequency responsiveness along with 10 their compact size that makes MESFETs so attractive for microwave amplifier applications in M M I C designs [1]. Also, in place of a full transient analysis of the equations, it is reasonable to use a steady-state analysis to approximate the behavior of the dynamic fields at operating frequencies which are typically well below the inverse of the transit times (in this case ~100 Ghz). Unfortunately, since this is a nonlinear problem, use of modal and frequency anlaysis in place of the transient analysis is not viable when dealing with higher operating frequencies. Source Gate Drain Figure 2. 2D cross section of a MESFET device. A MESFET device consists of a gate, drain, and source as shown in Figure 2. Under the drain and source metal electrodes, heavily doped ( > 1 0 ^ cm"^) n + regions ensure very thin Schottky barrier-induced depletion regions under the electrodes which are typically dealt with in problems by substituting them with ohmic contacts [2], [23] (See Section 5.3.1). The region under the gate is n-type with doping levels around 1 0 ^ - 1 0 ^ cm~^. Surrounding these active regions is a semi-insulating GaAs substrate which ideally has an intrinsic carrier concentration of 2.25 x 10 ^ cm~3. During operation, a current is directed from the drain to the source along a channel around the gate depletion region. By varying the gate voltage, the size of the Schottky barrier depletion region changes determining the amount of current and the operating state of the 11 device. FETs typically have two operating states, the active regime and the saturation regime. In normally-on, or depletion, devices the depletion region does not reach the buffer zone under equilibrium conditions, while in normally-off, or enhancement, devices this depletion region extends into the substrate buffer at equilibrium. For amplifier applications, a MESFET must be forced into the saturation regime, in which the electron channel is commonly referred to as being "pinched o f f between the semi-insulating substrate and the gate depletion region, by applying a negative bias voltage on the gate. This phenomenon will be demonstrated in a depletion GaAs MESFET through simulation using the F E M . 12 3 FEM Formulation The finite element method consists of approximating an unknown function, in this case potential and carrier densities, by a low order polynomial with unknown coefficients over a subregion, or element, in the solution domain. Al l the elements together span the entire domain, and qualitatively speaking, each element interacts with its neighboring elements through their common boundaries. The behavior of the function at the edges of those elements bordering the solution domain is pre-defined by the solution domain's boundary condition specifications. In solving for the polynomial coefficients by way of a matrix inversion, the boundary conditions are, in effect, propagated throughout the solution domain via element boundary interactions as dictated by the governing system equations. A brief discussion of this process is explained in the following sections, and details can be found in various books on finite element analysis such as Zienkiewicz [9], Buchanan [10], and Burnett [11]. 3.1 Elements For two-dimensional problems, an element has a polygonal shape. In practice, triangular and quadrilateral elements are the most popular because they are best suited for approximating boundaries of irregularly shaped domains. In addition to this, formulation of the polynomial approximations for these elements is straightforward as will be shown next. Figure 3. Six adjacent quadratic triangular elements. 13 Quadratic triangular elements are used in the F E M simulations that will be presented here. Figure 3 shows six such interconnected elements. Each is made up of six points, three at the vertices and the other three about half-way between vertex pairs. Over each element, a function g can be approximated with a function g a defined as where g n is the value of the function at node n, N n is a function of cartesian coordinates x and y for node n, and [N][g] is the matrix representation of the summation. N n is a basis function commonly referred to as a local shape function in F E M theory. The term "local" implies "per element". A "global" shape function for a particular node is simply the summation of all the local shape functions for that node. Figure 4 shows a linear (three-node triangular elements) shape function for a node both locally and globally. g - ga = 2 n g n N n = [N][g]- (8a) Figure 4. Local (left) and global (right) shape functions for Node i . From (8 a), the gradient of g can be written as V g - V g a = V[N][g] (8b) where V operates on each component of [N]. Throughout this paper, V[N] is a two-row matrix with the top row defined as d[N] / dx and the bottom row defined as d[N] / dy, and V [ N ] ^ is a 14 two-column matrix with the left column defined as 3[N]A / dx and the right column defined as d[N] T /3y. It follows that V [ N ] T = (V[N]) T . Derivation of N n depends on the type of element, and for triangular elements, the shape functions are often defined in terms of area coordinates rather than cartesian coordinates. A basic property of both local and global shape functions is that each N n has value 1 at node n and value 0 at node m for m*n. Also, local shape functions for nodes of each element or global shape functions for nodes of all elements must satisfy the following property: 2 n N n = 1 , (9) where N n is local, and n = 1, 2, . . . number nodes per element; or N n is global, and n = 1, 2, . . . total number of nodes. For a detailed discussion of area coordinates and the derivation of shape functions for various different elements, the reader is directed to Zienkiewicz [9] and Burnett [11]. As an extension to an exercise in Buchanan [10], a simple, yet general, method for determining shape functions for virtually any type of element is proposed in Appendix A . Although this method was not actually used in the F E M code here, it serves well as a complete and simple explanation of shape functions as an alternative to that which is given in most literature. A three node linear element with a three term expansion is the simplest triangular element. The higher the order of the element, i.e. the more terms used in (8), the better the unknown function can be represented over that element. Fewer higher order elements than lower order elements are needed to obtain the same desired accuracy in the final approximate solution. On the other hand, lower order elements require fewer nodes, and in either case, the exact solution will be approached as the elements become infinitesimally small. The choice of using quadratic triangular elements over linear triangular elements is based on the fact that they allow for the crudest possible expansion of the mobility gradient term in the advection-diffusion model of the continuity equations, while linear elements lead to a zero-valued mobility gradient term and are therefore unsuitable. A detailed explanation regarding this topic will be given in Section 3.2.2.2. 15 3.2 Discretization At this point, either a variational or a weighted residual Galerkin approach is used to complete the F E M formulation. In many cases such as this one, a variational form of the system equations is unavailable, so Galerkin's method is employed here. Most linear time-independent systems can be represented by an equation f(g,x,y) = fi(g,x,y) + f 2(x,y) = 0 and appropriate boundary conditions, where f^ is a linear function of g, and f j , f2, and g are functions of spacial coordinates, x and y. It is desired to determine the unknown function g, so it is replaced by an approximate expansion g a leaving the expansion coefficients g n as the only unknowns yet to be determined. Since f ^ is linear with respect to g, f(ga,x,y) can be written as _ n f 1 (N n ,x ,y )g n + f 2(x,y) = [f][g] + f 2(x,y), where [f] is a row vector of functions f j (N n ,x ,y) , and [g] is a column vector of the unknown coefficients. The residual is defined as R = f(ga,x,y) - f(g,x,y) = [f][g] + f2(x,y), and the set of coefficients [g] is sought so as to minimize the residual over the solution domain. This is done by setting weighted inner products of the residual to zero on an element-by-element basis which results in a matrix equation representation of the system equations for each element spanning a subdomain __e. Galerkin's method consists of using the shape functions as the weighting functions, and the result looks like this: /Q e[N] T([f][g] + f 2 (x,y))dA = [Ax)[g] + [F,] , where Q e is the region spanned by the element, [N]^ is a column vector of shape functions, matrix [Ai] is called the local stiffness matrix, and [Fi] is the local forcing vector of the element. The jth component of [g] corresponds to the approximate value of g at node j , the ith component of [Fj] corresponds to a force at node i , and the ijth component of [A]] can be interpreted as a coupling factor, or stiffness, between nodes i and j . A global stiffness matrix and forcing vector representing the whole system are formed by the superposition of the local stiffness matrix and 16 forcing vector component contributions of all the elements spanning the entire solution domain. For instance, the total, or global, stiffness between nodes 1 and 2 in Figure 3 is equal to the sum of the two local stiffnesses between these nodes which are defined in the local stiffness matrices of elements I and II. The resulting matrix equation is [A][g] = [F], where [A] and [F] are the global stiffness matrix and forcing vector, respectively, and [g] is now a column vector of all the unknown coefficients. After forcing essential boundary condition values on nodes located at domain boundaries in this equation, the solution for the coefficients is given by [g] = [A]-![F] For construction of [A], numerical quadrature is used to evaluate the inner product integrals. Gauss-Legendre quadrature is best suited for rectangular elements, while numerical formulae based on area coordinates are ideal for triangular elements (See [9], [10], and [11]). For large domains with many elements, [A] becomes very sparse, and a specialized matrix inversion routine that takes advantage of this property should be employed in finding [g]. Discretization of the semiconductor equations is described next. 3.2.1 Poisson Equation Using the notation in (8) to approximate the potential as [N][(j>] and the electron carrier density as [N][n] for each element, and inserting these expressions and (7) into (1) gives V-eV[N][<t>] = -q(N d-[N][n]) . Taking the inner products of this expression yields the matrix expression 17 fQe [N] TV-eV[N] dA [<j>] - q / Q e [ N ] T [ N ] d A [ n ] = - q / Q e [ N ] T N d dA . Using the Green-Gauss theorem on the first term, which is essentially integration by parts in more than one dimension, leads to the F E M expression of the Poisson equation: J Q e V[N] T eV[N] dA ft] + q[Qe [ N ] T [N] dA [n] = qfQe [ N ] T N d dA + J r e[N]TeV[N]-n dS [<|>] (10) where T e is the contour tracing __e, and n is the unit vector normal to T e . The last term is an integration along the outer contour of the element and represents a Neumann boundary condition on the potential, i.e. a Dirichlet boundary condition on the normal electric field component. This term is zero for insulating boundaries since, by definition, no drift current, which depends on the electric field, can exist across these boundaries. (10) is coupled to the continuity equations, and for each quadratic triangular element, it is a set of six equations in twelve unkowns, [(()] and [n]. 3.2.2 Continuity Equations Discretization of the electron continuity equation can be done in two ways yielding two different numerical models which are referred to here as the current conservation model and the advection-diffusion model. Each has its advantages and disadvantages in terms of implementation considerations, and the final solution algorithm makes use of both models. Following is a detailed derivation of each model. 3.2.2.1 Current Conservation Model The current conservation model was first introduced and demonstrated by Barnes and Lomaxin [23] and [22], respectively, and is re-derived and rewritten here with compact notation adopted from Zienkiewicz [9] and Buchanan [10]. Approximating the electron carrier density by 18 [N][n], and inserting this into (4) gives an approximate expression, J a , for the net current density: J a = q([N][n]u.eE + D nV[N][n]). Inserting this into (6) gives 0 = V - J a = qV-([N][n]u.eE + D n V[N][n]). Taking the inner product of this expression yields a matrix expression, and applying the Green-Gauss theorem results in the final form: Le [ N ] T V - J a dA = 0 L e V [ N ] T - J a d A = J r e [N] T J a -n dS (11a) q / Q e V [ N ] T - [ N ] L i e E d A [ n ] + q / Q e V [ N ] T D n V [ N ] dA [n] = J F e [N] T J a -n dS . ( l i b ) This model has the desirable property, from a circuit analysis point of view, that current is conserved regardless of mesh coarseness. This property is made possible by defining the surface integral term on the right as the net current directed out of an element through its nodes (See Section 5.2). For an interior node of a system of elements such as Node 1 in Figure 3, this term is set to zero indicating a zero net global current through the node. The same is done for exterior nodes at insulating boundaries. Thus, interior and insulator nodes obey Kirchhoff's current law, or I m = 0 where I m is the net current through each of these nodes. For an exterior node, c, at a non-insulating boundary of the domain, i.e. a contact boundary, there exists a non-zero net current, I c = / N c J a n dS, which is the sum of the local surface integral contributions and where N c is the global shape function for node c. Letting the total current through the contacts be written as 2^ I| and using the global form of (11a), an expression for the total contact current can be 19 derived in the following manner: hot = ^ i 1 ! = 2 c ! c + 2 i n I i n = / _ c N c J a _ d S + 0 = / V - j N i - J a d A . where I t o t is the total current through all the nodes, and Ij denotes net current through any node. Substituting the identity of (9) into this expression gives _ C I C = / _ c N c J a _ d S = / V - i N j - J a d A = / v ( l ) - J a d A = 0 . So the net current throughout the system is zero, and current is conserved globally. The surface integral in (11) represents a boundary condition on J a , which is the sum of both the drift and diffusion current densities. Again, for insulating boundaries, it has value zero since no current may exist across them. 3.2.2.2 Advection-Diffusion Model Another formulation of the continuity equations begins with inserting (4) into (6) and expanding the result with the chain rule: 0 = V - J n = V-q(n[i eE + D n V n ) = qnu. eV-E + qnVu, e-E + qu_Vn-E + q D n V 2 n (12) This equation is referred to as an advection-diffusion equation in Brooks and Hughes [27] and as the mass transport equation in Buchanan [10], which describes the fate of a concentration of some substance as it flows in some medium and, at the same time, diffuses into the medium. In this case the substance is the electron, and the medium is the conduction band in the semiconductor. The last term represents the diffusion, the third term represents the flow with u. eE as the electron velocity, and the first two terms represent an interaction between the electrons and their 20 surroundings. Using (1), (3), and (7) and assuming that the permittivity e remains constant over the span of an element allows substitution of V-E in the first term: 0 = qnu. e(q(N d - n)/e) + q n V u ^ E + qrA eVn-E + q D n V 2 n (13) This substitution removes some coupling between the continuity equation and the Poisson equation and is the key step in the derivation of the advection-diffusion model of the continuity equations. Now, the first term is nonlinear due to the n 2 term and the mobility u. e, which is a nonlinear function of the electric field intensity when taking velocity saturation into account. The second and third terms represent a heavy nonlinear coupling to Poisson's equation through both u.e and E , while the first term represents a nonlinear coupling to Poisson's equation only through u.£. Substituting the carrier density n with [N][n] in (13) and taking the inner products gives 0 = q/ Q e [N] T [N][n]Li e (q(N d -[N][n])/e)dA + q / Q e [N] T [N][n]V^-E dA + qL e [N] T u. e v [ N ] - E dA [n] + q D j o e [ N ] T V 2 [ N ] dA [n] ; and applying the Green-Gauss theorem to the last term leaves qL [N] T[N][n]Li eq(N d - [N][n])/e dA + q / Q e [N] T[N]Vu. £-E dA [n] + qfQe[N]TfieV[N]-EdA[n] - qD J Q e V[N] T V[N] dA [n] = - qD J r e [N] T V[N]-n dS [n] (14) The first two terms of this highly nonlinear equation offer several possibilities in the derivation of an iterative solution method as will be discussed for the application of Gummel's method (See Section 4.2). Some of these variations result in numerical convergence, and others may not. It is also important to address the mobility gradient Vu, £ in the second term of (14). In 21 building the stiffness matrix for a given element, Vu. e is calculated in the following manner: First, with the electric field intensity IEI = (V(|)-V(t>)1/'2 where Vcj) = V[N][(j)], the values of u,_(IEI) at each node of the element are approximated and are referred to as [[xe]; second, with these nodal values and the shape functions of the element, interpolation (See (8).) is used to write an approximate representation of \iQ as [N][u,e]; and finally, the expression for V[x g is given by V[N][ji e]. For a linear triangular element, the electric field, or -V(j), is a constant. It follows that Vu, g must be zero rendering linear triangular elements inadequate for approximating (15) on an element by element basis. Therefore, six node, quadratic triangular elements are applied to this problem instead, yielding a constant Vu. e for each element. Of course, the higher the order of the element, the better Vu, e can be approximated over the span of the element. The right hand term is a Neumann boundary condition on the carrier density n. Again, for insulating boundaries this term must be zero since any gradient in n would imply a diffusion current, and by definition no current can exist across such boundaries. For some boundaries it may be desired to force Dirichlet conditions on the potential (which implies a non-zero potential flux across the boundary) but also to force a zero net current condition across the boundary. In the current conservation model, this is achieved simply by setting the boundary condition on the current density to zero. To impose such a condition in the advection-diffusion model, the Neumann boundary condition on the carrier density, i.e. the diffusion current condition, must be set equal to the negative of the drift current. In other words, setting (4) to zero gives J n = q(nu.eE + D n V n ) = 0 q D n V n = - qnu,eE . So - q D j r e [ N ] T V [ N ] - n d S [ n ] = - qfTe [N]T[N]u.eV[N][c|)]-n dS [n] (15) for a zero net current across the boundary in the advection-diffusion model. 22 4 Solution Methodologies The Poisson equation and electron continuity equation together form a set of coupled, nonlinear equations to which iterative techniques must be applied in order to obtain solutions for [n] and First, a qualitative explanation of both direct iteration, of which Gummel's method [15] is a subset, and Newton-Raphson iteration in one variable is given. Following this are descriptions of Gummel's method and the Newton Raphson method in their application to the two semiconductor models derived previously: the current conservation model, (10) and (11), and the advection-diffusion model, (10) and (15). These iterative techniques may or may not successfully converge to solutions depending on the F E M model to which they are applied. Taking advantage of the best features of two combinations of model and iterative scheme, a solution algorithm is proposed that offers good robustness and performance. 4.1 Introduction to Iterative Techniques The following introduction to iterative techniques is based on a similar discussion in Zienkiewicz [9]. A single variable function, K x , for which K is a function of x is examined. Finding the values of x that satisfy Kx + f = 0, where f is assumed to be constant, can be done using direct iteration or the Newton-Raphson method. Direct iteration begins with an initial guess to the solution X Q . Next, the first iteration consists of finding a new solution Xj from x 1 = - f / K ( x 0 ) , which hopefully is closer to the actual solution than is X Q . Repeating this procedure, n iterations are carried out until x n converges to a particular solution. Gummel's method is a type of direct iteration scheme with the added simplification of decoupling the Poisson and continuity equations as will be seen in Section 4.2. Newton Raphson iteration begins with linearizing the equation 23 around some initial guess X Q by way of a truncated Taylor series expansion which gives the expression K ( x 0 ) x 0 + f + A(K(x 0 )x 0 )Ax 0 = 0 , where AXQ = X ^ - X Q and A ( K ( X Q ) X Q ) is known as the Jacobian which is, in effect, the function that is tangent to Kx at X Q . The first iteration involves first solving for AXQ and then for xj which is a new solution and guess for the next iteration. Again, repeating this procedure, n iterations are carried out until x n converges to a particular solution. Next, performance characteristics of these two schemes on three different types of functions Kx are shown graphically. Figure 5 shows how successive iterations of each method approach the solution of x for x = -fj . The initial guess is denoted as x = ag. aQ a^ a^ a^ a^ x aQ a^ a^ a^ x Figure 5. Direct iteration ; Newton-Raphson iteration. It is clear that the Newton-Raphson method has a much faster convergence rate than direct iteration for this case. Figure 6 shows the progress of a few iterations of both methods in the solution of x for K 2 x = -f 2 . 24 Figure 6. Direct iteration ; Newton-Raphson iteration. In this case, the Newton Raphson method still converges, but direct iteration does not converge to the solution regardless of the proximity of the initial guess to the exact solution. Finally, Figure 7 shows the behavior of each method for the K 3 X = -f3 case. Here, direct iteration converges to the solution, while the Newton-Raphson iterations possibly diverge unless the initial guess is close enough to the exact solution. Ideally in this situation, direct iteration could be used to reach a suitable initial guess after which Newton-Raphson iteration could be applied to rapidly reach the solution. The cases shown here will serve as analogies to the performances of Gummel's method and the Newton-Raphson method applied to the two different F E M models in the following two sections. 25 Figure 7. Direct iteration ; Newton-Raphson iteration. 4.2 Gummel' s Method Gummel first introduced his scheme in [15] as a general solution approach to semiconductor equations. It involves decoupling the equations and iteratively solving for [n] and [4>] independently. Gummel's method is actually a Gauss-Seidel method, which is described in most numerical analysis books such as Froberg [28], with specific application to the semiconductor equations. One iteration consists of two steps: 1) using an initial guess for the carrier density [n]1 in Poisson's equation, find a new potential [(J)]1+^ and 2) using [(J)]*+* in the continuity equation, find a new carrier density [n] 1 +^. The index i denotes the iteration number. These two steps are repeated until two consecutive iterations yield a pair of [n]'s and [<)>]'s with no significant changes between them. Such a change is defined by some norm, and a tolerance level is defined to be the maximum value of change for which the iterations may end (See Section 5.2). In terms of the iterative notation and leaving out the surface integral term, (10) becomes /QeV[N] T eV[N]dA[( |>] i + 1 + qfQe [ N ] T [N] dA [n]1 = q / Q e [ N ] T N d dA. 26 Assuming e is constant over each element and letting the doping density N d be represented by its nodal values as [N] [N d] allows this expression to be written in this compact form [Ai][<j>]i+1 = [B!]p([Nd] - [n]1) where [Ail = J Q e V [ N ] T V [ N ] d A ; (16) ' [Bjjp = (q/e) [B,] ; (17) [BJ = / Q e [ N ] T [ N ] d A . ' (18) It should be noted that if a continuous mathematical expression for N d is available, it should not be approximated by [N][N d] but instead should be employed directly in the integration over the element's domain. Building the global matrices [A] and [B]p from the local matrix contributions and incorporating any boundary conditions lead to a solution of [(()]^ +^ : = [ A r ^ p Q N d ] - [n]1) . (19) Formulating the continuity equation current conservation model, (11), in the iterative notation allows for three possible configurations (Again, the boundary term is not shown here.): ( [ C ^ + ^ m!!]^1 = 0 ; (20) (21) [C 1 ]i [n] i + 1 = -[Djjhi] 1 . [DjHnji+1 = -[q] i [n] i . ( 2 2 ) 27 where [D,] = q D J A J ; (23) [Ci] 1 = q / Q e V [ N ] T [ N ] M . e i E i d A = - q / Q e V f N ^ u . e M N H c } , ] 1 dA . (24) The number of possible configurations depends on the number of different stiffness matrices available. In this case, the two stiffness matrices [Ci] and [Di] can be arranged in three different combinations on the left side of the equation as shown in (20), (21), and (22) since 2^1 + 2^2 ~ 3, where n C m is the number of possible combinations of n elements in groups of m. It should be noted that D n is considered to be constant over the span of the element, and [Aj] is defined in (16). Putting (20) through (22) in global terms with [C] and [D], incorporating boundary conditions, and solving yields [ n ] i + 1 = ([Cji + tD])- 1 [F] ; (25) [ n ] i + 1 = - ( [ C ] 1 ) - 1 ^ ] [n]1 + [F]); (26) [ n ] i + 1 = - [ D r ^ W + [F]). (27) In the above equations, the forcing vector [F] accounts for the boundary condition contributions to the matrix formulation. As a quick evaluation, a simple problem consisting of a thin strip of GaAs with a metal contact having zero voltage at one end and an artificial boundary also set to zero voltage at the other end was set up to test the convergence of combinations of any one of (25)-(27) with (19). The doping density is uniform and the strip is partitioned into 100 quadratic triangular elements as shown in Figure 8(a). 28 Contact GaAs 0.01 L i m 1 [ x m ' Partitioning \ 5 \ 4 - . 3 E 2-2 >> 1 1 0 x 10 (a) C a r r i e r D i s t r i b u t i o n I t—— — •—i / / y / / / / / D e p l e t i Dn A p p r c o n - D i f f u . C o n s e r ivimatinr i 1 h- — — r A d v e c t C u r r e n s i o n M o d e l n a t i o n M o d e l 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 ,-4 x ( c m ) — > P o t e n t i a l F i e l d x 10 0.2 A 0 -0.2 To -0.4 £ o -0.6 o_ -0.8 D e p l e t i Dn A p p r c o n - D i f f u . C o n s e r iximatior A d v e c t i C u r r e n l 3ion M o d e l n a t i o n M o d e l 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 .-4 x ( c m ) — > (b) x 10 15 -3 Figure 8. (a) Schottky barrier strip and (b). results for N d = 5 x lO 1 - 3 cm The metal contact introduces a negative voltage boundary condition (28a) where (j)^ is the barrier potential of the semiconductor, § Q = kT/q ln (N c /N d ) , and V is the voltage 29 applied to the contact conductor. The barrier potential for GaAs is taken to be 0.8 V [22], and the contact conductor is at zero volts. The voltage boundary condition, in turn, also implies a carrier boundary condition n g just under the contact in the semiconductor. The relationship between ip s and n s is n s = N d e - q ^ s / k T . (28b) Explanation of (28) is given in Pulfrey and Tarr [26]. Through these boundary conditions, a Schottky barrier will be induced at the contact end of the strip. At the opposite end of the strip, the carrier density is set equal to the doping level N _ and the potential to zero. These are necessary artificial boundaries that simulate boundary conditions at infinity. In the absence of any educated guess of the carrier density, charge neutrality (n = N d ) , was used as an initial guess in the first iteration. The (19), (25) combination converged to a solution 1 C Q only when N d < 10 cm , otherwise no combinations of any one of (25)-(27) and (19) reached convergent solutions with higher doping levels, N d > 1 0 ^ cm"^. In fact, even when the initial guess of the carrier density was very close to the actual solution, all these combinations diverged from the solution. Thus, it is concluded that the current conservation model solved using Gummel's method is useless for any practical problems! Putting the advection-diffusion model, (14), which has four different stiffness matrices, into Gummel's iterative form results in fifteen (^ C^ + 4C2 + 4C3 + 4C4 = 15) different possible configurations, as opposed to the three possible configurations for the current conservation model shown in (20)-(22) and discussed earlier. After applying each of these configurations with (19) to the test problem defined earlier and shown in Figure 8(a), it was found that only two forms lead to converging iterations regardless of N d . The first is given as ([Ei] 1 + [Fi] 1 + [G]] 1 + [Dj]) [ n ] i + 1 = [0,]^] [ n ] i + 1 = ([E] 1 + [F] 1 + [G] 1 + [D])" 1 [ G ] 1 ^ ] (29a) 30 and the second is slightly different: ([Ei] 1 + [G,] 1 + [D,]) [ n ] i + 1 = [Gi]i[Nd] - [F i^n] 1 ; [ n ] i + 1 = ([E]1 + [G]1 + [D])" 1 ([G]i[Nd] - [F^n] 1) ; (29b) where the corresponding local stiffness matrices are [E]] 1 = q/„ e[N] Tu. e 1V[N]-V[N][(j)] 1 dA ; [Fi] 1 = qj*Qe [N]T[N]Vu._ i-V[N][(|)] i dA ; [Gi] 1 = (q 2 /e)/ Q e [N]T[N][n]Vei[N] dA ; and [Di] is defined in (23). Note that e is assumed to be constant over the span of the element, E has been replaced with -VfNjhtj]1, and [Ai] is defined in (16). It is also important to note that for the case of a zero net current boundary condition in (29), (15) contributes to the stiffness matrix, not the forcing vector! Both combinations of (19) and one of (29a) and (29b) are very robust, and in all cases (varying N d and contact voltages) tested on the strip, charge neutrality gave a sufficiently good initial guess of the carrier density to ensure convergence in the Gummel iterations. (29b) tended to reach a solution in slightly fewer iterations as compared to (29a). On the other hand, this scheme is more computationally intensive since [Ei] 1, [Fi]1, and [Gi] 1 must be recalculated for each iteration. Although it has not been done here, inclusion of a generation-recombination term would create an additional nonlinear term as a function of carrier densities that would also have to be recalculated for each iteration. Tan et al. [29] introduce an application of Scharfetter-Gummel (S-G) interpolation in conjunction with a seven-point Gaussian quadrature for the construction of a stiffness matrix and forcing vector for triangular elements that lead to a higher accuracy and convergence rate of the 31 solution when compared to the classical Gummel finite element formulation. The S-G interpolation is only used when applying the quadrature on the potential and carrier values of the previous iteration and therefore requires relatively minor modification to the original classical formulation in a computer program. Application of this scheme was briefly tested on the strip problem as well as the MESFET examples in Section 5.2, but there was no noticeable improvement in the accuracy or convergence rate of the solution. Nevertheless, [29] shows that the Scharfetter-Gummel scheme really shines for higher doping conditions, and due to its easy implemention in already existing F E M formulations that apply Gummel's method, it is well worth mentioning here. For later reference, the current conservation Gummel approach and the advection-diffusion Gummel approach (using (29b)) are abbreviated as C C G and A D G , respectively. 4.3 Newton-Raphson Method The Newton-Raphson method applied to the semiconductor equations consists of iteratively solving the Poisson and continuity equations simultaneously and implicitly. Its application to solving the current conservation model is derived here. The general form of a discretized nonlinear system can be written as [K([a])] = [P([a])][a] + [F([a])] = 0 (30) where ([a]) indicates that the matrix or vector is a function of the coefficients [a]. Linearizing [K([a])] by way of a truncated Taylor series expansion around some point [a]1 gives [K([a])] = 0 - [K([a]i)] + [K([a]i)]T A[ap . Defining [a]1 as the result of iteration i allows an iterative process to be described by A[a]i = -([KttaJ^Tr^Kaa] 1)] , (31) 32 [a ] i + 1 = [a]1 + A[a][ (32) where [K([a])]T = d[K([a])]/d[a] . (33) (31) and (32) describe one iteration of the Newton-Raphson procedure. Also, it is important to note that Dirichlet boundary conditions are applied to (31) by enforcing a no-change condition for a boundary node value, i.e. Aaj = 0 where j refers to a Dirichlet boundary node. The core of this method is the Jacobian or tangential matrix [K([a]1)]'p . Using (30), (33), and the chain rule leads to [K([a])]xd[a] = d[K([a])] = [P([a])] d[a] + d[P([a])] [a] + d[F([a])] (34) Rewriting the last two terms as d[P([a])][a] = [X([a])]d[a] (35) d[F([a])] = [Y([a])]d[a] (36) and substituting into (34) yields an expression for [K([a])]j : [K([a])]T = [P([a])] + [X([a])] + [Y([a])] . (37) Combining the Poisson and current conservation continuity models, (10) and (11), into one local matrix representation gives 33 [K([a])] = [P([a])][a] + [F([a])] = [A,] [B, ] p " m ' [B , ] p [N d ] " 0 [CJ + tD,] [n] + 0 = [ 0 ] (38) where [Ai ] , [Bi]p, [Ci] , and [Di] are defined in (16), (17), (24), and (23), respectively. In this case the [F([a])] term is not actually a function of [a], and it follows from (36) that [Y([a])] = 0. [P([a])] on the other hand has the drift term [Ci] which is a function of [())], thus d[P([a])][a] = • o o • •[*]• 0 d[CJ [n] (39) So d [Cj][n] must be evaluated in order to find [X([a])J as defined by (35). Following are the details of a full derivation of d [Ci][n] into the form, [Qi] d[(j>] : Through (24) and the chain rule d[C,][n] =-q/ n e V[N] T [N][n](djx e )V[N][(J)]dA - q / Q e V[N] T[N][n]u. eV[N] d[(|)] dA (40) Now du,e is evaluated: du,e = ( d[i e / d[(J)]) d[<j)] = ( du.e I dlEI) (dlEI / d[(j>]) d[(J>] = (du, e / dlEI) (d(V(J)TV(t)) 1 / 2 / d[(j>]) d[<J>] = (d|x e / dlEI) (d(V(()TV(l) ) 1 / 2 / d(V<t>)) ( d(V(|>) / d[(j>]) d[<|)] = ( du.e / dlEI) ((V(])TV(t) Tm V())T) ( d(V[N][(|>]) / d[<|)]) d[<|)] = ( d | i e / dlEI) (V(()T / IEI) (V[N] d[(j)] / d[<\>]) d[<|)] = ( d|xe / dlEI) IEI"1 (V[N][(|)])T V[N] d[$] Substituting this result into (40) gives 34 d [q][n] = - q / Q e V[N]T-[N][n] (( dLi e / dlEI) IEI"1 (V[N][(J)])T V[N] d[(J)]) V[N][(|)] dA -qLeV[N]T-[N][n]^ieV[N]d[(j)]dA . Rearranging terms leaves d[q][n] = [Qi]d[(H = " qLe V [ N ] T - V[N][(|>] [N1W ( d^ e 7 d l E I ) I E I " ! (V[N][<|>])T V[N] dA d[<|>] - qfQ eV[N] T-[N][n]Li eV[N]dAd[(t)] . Now, (39) can be rewritten as d[P([a])][a] = [X([a])]d[a] = • o o • W • o o- •dm 0 d[C,] [n] [Q,] 0 d[n] It follows from (37) that [K([a])] x = [P([a])] + [X([a])] + [Y([a])] = [A,] [B,] p [Q,] [C,] + [D,] (41) (38) and (41) are used to build the global Newton-Raphson formulation as defined by (31) and (32): "A[<Mi+1" • [A] A[n] i + 1 [Q]\ [C] i +[D] [A] [B] p 0 [C] j + [D] _1_ " [B] p[N d] • \ w1. I 0 I '[ct)]i+1' [n] i + ' = + A[n] i + 1_ w1. (42) Note that only [C] and [Q] must be computed for each iteration. On the other hand for an N-node problem , 2N equations must be solved simultaneously at each Newton-Raphson iteration as 35 opposed to having to solve two sets of N equations at each Gummel iteration. This places more demand on computing power both in terms of memory and speed. It should also be noted that any nonlinear term, possibly to account for generation-recombination in the continuity equations, would follow the same basic procedure as that of the nonlinear drift term for determination of [K([a])] x. This Newton-Raphson current conservation technique, (42), does converge to a solution given a "suitable" initial guess, as opposed to the implementation of Gummel's method to the same model which, for doping levels beyond a realitvely low point, diverges from the solution regardless of the initial guess (See Section 4.2). Unfortunately, using charge neutrality and zero potential everywhere as an initial guess most often leads to diverging results and is therefore not "suitable". On the other hand, for a good intial guess (42) converges very rapidly to the solution. In fact, Newton-Raphson methods converge quadratically near the solution, while direct iteration methods converge only linearly as the solution is approached [28]. Combining the robustness of the Gummel approach to the advection-diffusion model with the efficiency of the Newton-Raphson approach to the current conservation model yields a good system for solving the coupled Poisson-continuity equations. As it turns out, relatively few iterations of Gummel's method are required to achieve an appropriate initial guess for which the Newton-Raphson method can then be implemented to quickly converge to the solution (See Section 5.2). Applying this procedure to the strip test of the previous section requires a total of about seven iterations (four Gummel iterations and three Newton-Raphson iterations) as opposed to nine iterations required when applying only A D G . Figure 8(b) shows the carrier density distribution and potential field solutions to the strip problem for both models as well as analytic solutions based on the depletion approximation which is discussed in most semiconductor theory texts [26], [2]. Throughout the remainder of this paper, the Newton-Raphson current conservation approach and the combined Newton-Raphson current conservation and Gummel advection-diffusion approach are referred to as C C N R and A D C C , respectively. Applying the Newton-Raphson method to the advection-diffusion model was not considered here because there is no apparent advantage it could offer. On observation of (14), it is 36 clear that a Newton Raphson approach would be quite computationally intensive, much more so than either the Gummel approach to the advection-diffusion equations (ADG) or even CCNR, due to the many nonlinear terms that would have to be calculated at each iteration. In other words, even for a best case scenario in which it offered better convergence properties than CCNR, its computational inefficiency would be more than enough to diminish its usefulness. Also, the current conservation model has the nice property of current conservation which the advection-diffusion model does not. Qualitatively speaking, the current conservation model's stiffness matrix is analogous to K 2 in Section 4.1, and Figure 6 shows iteration behavior analogous to that of the C C G and CCNR approaches when applied to the Schottky barrier strip test. Also, the advection-diffusion model's stiffness matrix is analogous to either K j or K3 for which Gummel's method proves to be very robust. The ADCC scheme ensures convergence by first applying direct iteration (Gummel's method) to the or K3 case (advection-diffusion model) until a reasonable initial guess is achieved before implementing Newton-Raphson iterations on the K 2 case (current conservation model) for faster convergence. 37 5 Application Examples of application of the previously described A D G and ADCC methods are demonstrated here. In each example, the methods are used to compute the potenital field and carrier distribution throughout a particular MESFET device. From these results various electronic properties of the device can be derived such as currents, conductances, and capacitances. Example 1 demonstrates current and small signal conductance calculations based on the analyses in [22]. This example is also used to compare the performances and results of A D G and A D C C . Example 2 is a full analysis of a MESFET which describes, in detail, the derivation of a small signal circuit model from the potential field and carrier distribution computed with A D C C . Preceding the examples, in the following section, certain restrictions regarding the mesh coarseness of the F E M problem are specified that ensure, to a certain degree, convergence and accuracy of the computed solution. 5.1 Mesh Spacing For numerical stability and physically meaningful results, general statements can be made concerning restrictions on mesh coarseness. The Debye length is the characteristic decay length of charge in a semiconductor (See [26] and [2]) and is given by L D = ( k T 8 / q 2 N d ) 1 / 2 . Qualitatively speaking, if the lower order polynomial expansion of an element that has dimensions of the order h is to reasonably approximate a decaying charge distribution, it is safe to say that h should not be larger than the Debye length. This restriction is necessary to achieve physically meaningful results. The phrase "dimensions of the order h" refers to the approximate distance between any two nodes of an element. 38 For numerical stability, another limitation on the size factor h of an element is Reynold's number: R c = vh/D =s 2, where v is the velocity and D is the diffusivity of the carriers. In mathematics, this is also referred to as Peclet's number and is of importance in flow problems in general. Attempts have been made [27] to remove this restriction by modifying the weighting functions used in the standard F E M formulation; these are known as Petrov-Galerkin techniques. Put simply, the standard Galerkin F E M discretization (as demonstrated in this paper) can be shown to, by its very nature, add artificial negative diffusion to flow equations which, when combined with large drift constants, i.e. a large Peclet number, may lead to spurious node-to-node oscillations or "wiggles" in the computed solution. These wiggles are strictly numerical artifacts and, in the case of nonlinear flow problems, may even degrade the convergence and stability of an iterative solution scheme. The Petrov-Galerkin F E M discretization "solves" this problem essentially by adding positive artificial diffusion to cancel the negative artificial diffusion in the standard formulation. Unfortunately, such diffusive schemes lower the order of accuracy of the numerical solution from second order to first order. In fact, it is questionable to use such wiggle-suppression schemes at all (See Gresho and Lee [29].), and if wiggles do appear, simple local mesh refinement is probably the safest remedy. Fortunately, in the case of GaAs the drift velocity of the electron carriers cannot exceed approximately 2 x 10^ cm/s (barring velocity overshoot or any other quantum phenomena [26]), so Reynold's number does not dictate the mesh fineness for the submicron gate MESFET problem as much so as does the Debye length. 5.2 Simulations Following are two examples of F E M simulations for submicron gate MESFETs. Similar to the demonstrations in [22], the first example briefly analyzes a MESFET for the main purpose of 39 comparing the performances of the A D G and ADCC schemes. The second describes a complete analysis of another MESFET with a thinner active region of higher doping concentration in which conductances and capacitances of the device at a certain bias are predicted. In both examples, doping levels are relatively low in comparison to practical scenarios. This is done in order to keep the number of nodes at "reasonable" levels for implementation on a conventional workstation. The F E M algorithm is programmed entirely in Professional M A T L A B v4.2, and all F E M simulations are carried out on a Sun SPARC 10 workstation. Before proceeding to the examples, some details concerning current calculations, norms, and the electron mobility expression are specified here. As in Barnes and Lomax [23] following the idealogy of the current conservation model, currents [I] and charges [Q] travelling through and existing on all nodes are expressed conveniently as [I] = J r [N]TJ a-n dS = ([C] + [D])[n] (43) and [Q] = J r [N] T eV[N]-ndS[(l)] = s[A] [® + e [ B ] p [n] - q / Q [ N ] T N d dA (44) where local stiffness matrices [Ci], [Di], [Ai] , and [B]]p are defined in (23), (24), (16), and (17), respectively. These representations simply state that a node's current and charge are the accumulations of weighted normal component contributions of electric current density and electric flux density, respectively, over the portion of the boundary, T, spanned by that node's global shape function. It should be clear that these definitions satisfy the integral form of Maxwell's equations [30] for steady-state current and charge as node spacing becomes infinitesimal, i.e. / r Jn dS = I and J r D n d S = Q 40 For comparing node potentials resulting from the last two iterations, a scaled Euclidean norm is used: A<|> = (l/maxCM1))- (2 n ( V - K i A ) 2 ) m For comparing node carrier densities from the last two iterations, the following norm is used: An = ( 2 n ( l - n ^ V ) 2 ) 1 7 2 In contrast to the scaled Euclidean norm, this one places more importance on local nodal changes regardless of the values of the carrier densities at those nodes. Since carrier densities vary on exponential scales, this is a more suitable norm. In the above definitions, i is the iteration number, and n is the node number As in [22], the following expression is used for the nonlinear mobility term: Li e = (u,0 + v sIEI 3 / IE 0I 4) / (1 + IEI 4 / IE 0 I 4 ) , (45) where u.Q is the low field intensity electron mobility, v s is the electron saturation velocity, and IEI is the electric field intensity. The constants are chosen as u.Q = 5000 c m 2 s~' V ~ ' , v s = 0.85 x 10^ c m s ' ' , and IEQI = 4000 V cm"' . Figure 9 shows a plot of the electron drift velocity with respect to electric field intensity. 41 0 2 4 6 8 10 12 14 16 18 20 Electric Field Intensity (kV/cm) —> Figure 9. Electron drift velocity versus electric field intensity. 5.2.1 Example 1 Figure 10 shows the configuration of the cross section of the GaAs MESFET to be analyzed. The cross section is 1 x 2 um, and the gate length is 0.4 um. There is an active region with a uniform doping concentration of N d = 5 x 1 0 ^ cm"^, and a buffer zone with a uniform doping concentration of N d = 1 x 1 0 1 J cm" J and thickness 0.2 um. In order to reduce the number of necessary nodes, the Schottky barriers under the source and drain are neglected, and ohmic contacts are assumed in place of heavily doped wells that have the same effect. The mesh used for this problem was generated using a semi-automatic mesh generation algorithm based on Zienkiewicz and Phillips [31] and is shown in Figure 11. As can be seen, the mesh is finer near the contact edge singularities and underneath the gate since these are the regions of maximal gradients in carrier density and potential field. 42 0.3 \im 0.4 Lim G 0.3 Lim D E =5. N , = I 0 1 3 c m " 3 2 Figure 10. MESFET device for Example 1. 0 0.2 0.4 0.6 0.8 1 1.2 1.4 1.6 1.8 2 Figure 11. Mesh Grid for Example 1 - 959 nodes, 448 elements. Dirichlet boundary conditions exist along the contacts and just below the buffer zone. The boundary potentials at the source and drain (ohmic contacts) are set equal to the voltage on the source and drain contacts, V d and V s , while the boundary potential along the gate boundary is set to -0.6852 V using (28a). V s = 0 V is used as the reference voltage throughout this example. The boundary carrier densities at the source and drain contacts are set equal to the doping level forcing charge neutrality as required for ohmic contacts, while the carrier density n s just under the gate 43 contact is 1.62 x 10^ cm"^ satisfying (28b). The boundary under the buffer zone is given an artificial condition of zero potential and charge neutrality to approximate the effects of a deep substrate. Of course, the artificial charge neutrality condition results in "leaking" current through this artificial boundary, and current conservation is not preserved. This is done intentionally in order to keep the number of A D G iterations low. It turns out that inclusion of the term (15), which is required to specify zero net current across a boundary in the advection-diffusion model, increases the number of iterations necessary to reach a solution with A D G , in this case by a factor of four. Neumann boundary conditions are enforced on the potential and carrier density at the semiconductor-air borders between the contacts, and along the left and right edges of the device. These boundary values are zero suggesting no current, diffusion and drift, across the boundaries, i.e. perfect insulation. Applying these conditions is accomplished simply by neglecting the surface integral terms in the F E M formulations for both A D G and A D C C . Figures 11 and 12 show the potential field and electron carrier distribution calculated with the ADCC method for V s = V g = 0 V and V d = 3 V . The Schottky barrier spans most of the region under the gate leaving only a thin channel of a high concentration of electrons just underneath. The channel is "pinched o f f at the drain side by a depleted region of carriers. Figure 13. Electron carrier distribution. 45 Figure 14 shows a plot of the current density field, and it is clear that the concentration of current is greatest along the narrow channel underneath the barrier. Also, after closer inspection of 1 C Q Figure 13, a "bump" indicating an accumulation of carriers above the 5 x 10 cm" J level exists along this channel. This bump and the adjacent valley of carriers constitute what is known as a stationary or Gunn domain, and is a result of velocity saturation effects as predicted by one-dimensional FET analysis [26]. As the drain voltage is increased further in the saturation regime, the stationary domain absorbs the bulk of the extra potential difference along the channel. x 10" 12 10 \ \ \ v \ w \ \ \ s \ \ \ \ \ \ \ \ \ \ \ \ \ \ •.W-A w \ \ \ \ \ \ \ \ \ \ \ \ ' ' ' ' / / J I ' ' * ' S / i I ' ' ' s s / / 1 I \ \ \ \ \ \ S \ \ \ \ \ s \ \ N \ \ N E ^ 4 -2 > N S V -, V V X X v . \ *V v. . \ \ N N V • 10 15 x (cm) x 10 20 5 Figure 14. Current density field. Simulations were carried out for varying gate and drain voltages. A plot of I D vs. V _ , for Vg = 0, 0.5 V , is given in Figure 15. As expected, there is an active regime and a saturation regime of operation. There is also a region between these that shows negative resistance. Such regions are typically avoided in the design of amplfiers since they may cause the device to oscillate. Figure 16 shows a plot of 1^ versus Vg for V d = 2.5 V . 46 3001 r Drain current versus drain-to-source voltage. 200 Vgs.= 0 V + X ADCC ADG Figure 15. Lj vs. V d for Example 1. Drain current versus gate-to-source voltage (Vds = 2.5 V). Figure 16. I d vs. Vg ( V d = 2.5 V) for Example 1. From these two plots, a quick calculation of the voltage gain G of this device operating in saturation mode, or V d > 2 V , and at V„ = 0 can be made: G = 0 I d / dV ) / ( 3 l d / d V d ) - 116.3/7.2 - 16. Figure 15 also shows drain current versus drain-to-source voltage results for both the ADG and the A D C C methods for V d = 0,0.5 V . The A D G and ADCC curves closely match each other 47 in the saturation regime but differ considerably in the triode regime and between the triode and saturation regimes. It is assumed that the crude finite element approximation of the term containing the gradient of the mobility, V u , e , in the advection-diffusion equations may introduce significant error in representing the drift current term, particularly when the electric field intensity along the channel is in the range of the "hill" in Figure 9, i.e. where [x e is highly nonlinear. Since the current conservation model does not suffer from such approximations and has fewer terms to approximate than the advection-diffusion formulation, the ADCC plot is assumed to be a better approximation of the true current. For this problem, the A D G method required between thirty and forty iterations, while the A D C C method required only about eight iterations for change tolerances < 10"^. Of the ADCC iterations, the first two or three were A D G iterations and the remaining were CCNR iterations. It should be emphasized that these results as well as those for the strip test in Section 4.2 and shown in Figure7(b) come from two completely different F E M formulations, the current conservation model and the advection-diffusion model, solved by two completely different techniques, the Newton-Raphson method and Gummel's method, respectively. So, the fact that the A D G and A D C C results agree quite nicely, particularly in the saturation regime, confirms with some certainty that the results are indeed valid as approximate solutions to the MESFET problem. 5.2.2 Example 2 Next, a device similar to the one in the previous example is fully analyzed. A "full analysis" constitutes predicting the small signal characteristics of the MESFET for the purpose of constructing a small signal circuit model useful for implementation in circuit analysis applications. The circuit model used here is adopted from Pulfrey and Tarr [26], and a similar model is presented by Curtice [32]. The device used in this example is the same as the one described in Example 1 with a few modifications: The active layer is 0.4 um thick with doping density N d = 2 x 1 0 ^ cm"^, and the buffer zone is 0.6 um thick. Also, rather than forcing a neutral charge condition, a zero net current condition is forced on the lower artificial boundary. This ensures current conservation in the 48 system assuming the current definition of (43). The mesh used is shown in Figure 17. It is made up of 1453 nodes and 688 elements. Potential and carrier density field solutions were computed for a series of bias conditions using A D C C . The number of iterations required to reach the solutions starting from an initial guess of zero potential and neutral charge is about forty, of which only the last three or four are CCNR iterations. For a sufficiently good initial guess, CCNR iterations reached a solution within three or four iterations. x IP" 4 S o u r c e G a t e Drain x ( c m ) — > X1Q-4 Figure 17. Mesh grid for Example 2 - 1453 nodes, 688 elements. Figure 18 shows how the potential and carrier fields change as the drain voltage is varied. As the drain voltage is increased from zero voltage (Figure 18(a)-(b)), the depletion region expands on the drain side and begins to pinch-off the channel of carriers. In Figure 18(c)-(d), the channel is completely pinched off and there is relatively little change in the carrier density distribution as the drain voltage is further increased. As a result, the drain current increases very slowly for higher drain voltages, and the device is described as being in the saturation regime. It should be noted that, unlike in Example 1, there is no Gunn domain present in this device. Typically, very small scale MESFET devices do not have this trait, the reason for this being that carrier transit times are so small such that charge accumulations do not have enough time to form before the carriers leave the device. A small carrier wall exists along part of the artificial boundary in all cases. It is a result of satisfying the zero net current boundary condition across this boundary, the potential for which 50 en 's o > + T3 z CO -a > ll c o tu 3 CJ) 51 has a Dirichlet condition. The carrier wall supplies the necessary diffusion current that cancels out the drift current caused by the potential gradient on the boundary and is therefore solely a by-product of the artificial boundary condition with no physical meaning. A plot of the current density field is given in Figure 19, and as expected, the majority of the current exists mainly along a narrow channel from the drain to the source and between the Schottky barrier-induced depletion region and the buffer zone. x 10 Current flow for Vds = 4 and Vgs = -0.5 12h 10 E & 4 -2 • 1 ,'V- • I 1 \ \ - . 1 I I I 1 10 15 x (cm) Figure 19. Current density field. 20 x 10" Drain current versus drain voltage plots for different gate voltages are shown in Figure 20. As expected, both a triode and saturation regime are apparent from these plots. Unlike the device in Example 1, there is no region of negative differential resistance in any of the current plots. It is believed that there is a correlation between the existence of a Gunn domain and the appearance of a negative differential resistance in the plots. Under certain bias conditions and depending on its dimensions (particularly the gate length) and doping levels, a MESFET device may behave similarly to a transferred electron device (TED) operating in the limited space-charge accumulation (LSA) mode [33]. TEDs, also known as Gunn diodes, are typically also made of GaAs and are 52 used for their negative resistance characteristics in amplifier and oscillator applications. For a small-signal characterization, the bias conditions of interest are chosen to be V g S = -0.5 V and V d s = 4 V . A l l small-signal conductances and capacitances are sought in order to build a small-signal model of the MESFET. From the data used to create the plots in Figure 20, the small signal conductance g d s is roughly approximated as g d s = dId/dVds - A I d / A V d s = 2.06/1 mho/m = 2.06 mho/m . (46) where A I d = Id(V d s=4.5) - I d (V d s =3.5) and A V d s = 4.5 -3.5 = 1. Next, the drain current must be calculated for different gate-source potentials. Figure 21(a) shows a plot of the results. From this data the transconductance g m can be computed: g m = dld/d\gs ~ AI d /AV„ s = 16.17/1 mho/m = 16.17 mho/m . (47) Similarly, the conductance g „ s is given as (See Figure 21(b).) 53 g = ai„/aV - Al / A V g s = 15.97/1 mho/m = 15.97 mho/m . (48) Drain current versus gate-to-source voltage (Vds = 4 V). Gate current versus gate-to-source voltage (Vds = 4 V). Vgs(V) —> Vgs(V) —> (a) (b) Figure 21. (a) Drain and (b) gate current versus gate voltage. By linearizing the currents and voltages around the DC bias points, the method outlined in Appendix B can be applied for finding the small signal capacitances of the device. This is done simply by replacing [V] and [Q] with perturbations A[V] and A[Q], respectively, so that the resulting capacitance matrix satisfies A[Q] = [C]A[V], where A[Q] = [5Q d 5 Q g 5 Q S ] T and A[V] = [5V d 5 V g 5 V S ] T . Using (44), electrode charges are computed from the carrier and potential solutions for several electrode potentials in the vicinity of the bias voltages on the drain, gate, and source. Figure 22 shows plots of the charges on each electrode for varying electrode potentials. 54 Electrode charges versus drain-to-source voltage (Vgs = -0.5 V). + Qd x Qg o Qs 3.5 4 4.5 Vds (V) —-> (a) Electrode charges versus gate-to-source voltage (Vds = 4 V). -X Qg O Qs E O -2.5 -2 -1.5 -1 -0.5 Vgs(V) - - > (b) Electrode charges versus source voltage (Vds = 4 V, Vgs = -0.5 V). + Qd X Qg 0 Qs -10 ' ' 1 1 -1 -0.5 0 0.5 Vs (V) —-> (C) Figure 2 2 . Electrode charges for varying potentials on (a) the drain, (b) the gate, and (c) the source. 55 From the data in the plots, the capacitance matrix is constructed: [C] 100.04 26.66 -179.22 -22.7 288.9 -280.52 -0.26 -3.28 6.6 pF/m As discussed in Appendix B , the diagonal terms give the capacitances to ground. Since this is an active device, this capacitance matrix is not symmetric. So, the next step is to decide which off-diagonal terms are to be used to represent the mutual capacitances between the drain, gate, and source terminals. Assuming that the device is in common source mode, the gate acts as the input, and it is desired to represent the behavior of the charges on the gate as they redistribute themselves in response to a change in the applied input voltage. This behavior can be expressed by the capacitance elements C21 and C 2 3 [34]. Also, since the manner in which the charge redistributes itself at the output (the drain terminal) in response to a change in the drain voltage is sought, the capacitance element C i 3 is used as the mutual capacitance between the drain and the source, which in this case is an additional capacitance to ground. The resulting mutual capacitances and self-capacitances for this MESFET device in common source configuration are given as follows: C G S = 280.5 pF/m C G D = 22.7 pF/m C D S = 179.2 pF/m C G G = 288.9 pF/m C D D = 100.0 pF/m (49) Compiling the results of (46)-(49) yields the small signal circuit model shown in Figure 23. An exponential term is introduced in the current source expression to take into account the time it takes for an electron to travel from the source to the drain. This delay time is taken to be approximately 10 ps, which is the time required for an electron to travel 1 um at a velocity of 1 x 10^ cm/s. x = transit time (~ 10 ps) R G and R D = parasitic resistances g m = 13 mho/m g g s = 16.2 mho/m gjg = 2.1 mho/m C G G = 288.9 pF/m C D D = 100.0 pF/m C D S = 179.2 pF/m C G S = 280.5 pF/m s C G D = 22.7 pF/m Figure 23. Small signal model of MESFET in Example 2. It should be noted that the rough approximations of the partial derivatives in the above analysis were used in the context of a simple example with the intention of demonstrating direct application of the sampled data points in the plots. For applications in circuit analysis software, such as SPICE, each set of sampled data could be replaced by a polynomial curve fit. It follows that all conductances and capacitances could then be represented by continuous polynomial expressions as functions of the terminal voltages. 57 6 Extensions to FEM Model In monolithic microwave integrated ciruits (MMICs), MESFETs are implanted into an intrinsic GaAs substrate, and strip-like conductors feed the source, gate, and drain of the devices. Typically, source and drain strips are on the order of tens of microns in width, extending far past the feeds of the MESFET's active region (See Figure 24). These strip conductors introduce considerable capacitances of their own that are not completely accounted for in the previously described F E M analysis. These capacitances are parasitic degrading the performance of the device in terms of speed and power consumption [32], [35]-[37]. Following are descriptions and recommendations of how Green's integral equations and boundary integral equations (BIE) solved with the method of moments (MoM) can be used in conjunction with the F E M model of the MESFET to give a more accurate model of the device. The substrate and ground plane shown in Figure 24 extend infinitely along the x-axis. GaAs Substrate x Figure 24. MESFET with strip conductor feeds. 6.1 Method of Moments and Integral Equations The method of moments (MoM) is a very efficient numerical technique originally 58 introduced by Harrington [5] for finding approximate solutions to integral equations of the form: where r is a position in space, a(r) is known in some domain Q, b(r') is unknown and to be determined, and c(r, r') is a known function, typically a Green's function for a particular linear differential equation satisfying a set of boundary conditions. The M o M consists of expanding the unknown b(r') into a sum of N known basis functions with N unknown coefficients, setting up a set of N equations for N known values of a(r), and solving for the unknown coefficients by a matrix inversion. 6.2 Parasitic Capacitances and Green's Integral Equation The strip electrodes of a MESFET can be thought of as having free charges sitting on their surfaces when bias potentials are applied to them. As in [6], the goal is to find these charges in order to determine the parasitic capacitances of the strip electrodes which are then simply added together with the MESFET capacitances resulting from the F E M simulation in the circuit model. The big approximation and assumption that should be noted here are that the GaAs substrate is essentially charge-free (a perfect insulator) and that the strips and that the charge-free air-substrate-ground configuration can be decoupled from the active region of the MESFET. This reduces the problem to a linear Poisson equation which governs the behavior of the static field resulting from a set of free charges restricted to some finite domain (the strip surfaces). A general solution of the potential field for the linear Poisson equation is the Green's integral equation . [5]: (50) 59 where (|)(r) is the potential field, p(r') is the charge distribution on some domain Q, and G(r, r') is the Green's function of the system. In this case, the strip surfaces constitute the domain of the problem, and the Green's function [38] satisfies the boundary conditions that exist at the substrate interfaces and the condition of zero potential at infinity. Since the potentials on the electrodes and the Green's function are known, the charge distribution can be found using the M o M . As in , New York: Macmillan, 1968. [6], the strip conductors are approximated to be infinitesimally thin. In contrast to , New York: Macmillan, 1968. [6], the problem is solved in two dimensions rather than three, and complex images and a multipipe model of the strips are used to give a much more compact and efficient M o M solution to the problem rather than using conventional image theory and pulse basis functions. Brief overviews of these techniques follow. 6.2.1 Complex Images Complex images [39] are used for the purpose of determining an approximate closed-form Green's function for a multilayered substrate. For the case of a point charge above a ground plane, conventional image theory states that the potential field above the ground plane is identical to the field for the case in which the ground plane is replaced by a mirror image of the original charge with a sign opposite that of the original charge. This conveniently allows the Green's function of the original system to be represented as the sum of two free-space Green's functions which have simple, known analytic expressions. When dealing with a point charge embedded in a multilayered (more than two layers) structure, conventional image theory (, New York: Macmillan, 1968. 60 [6]) leads to representation of the Green's function by a singly (three layers) or multiply (more than three layers) infinite series, which often converges very slowly. The use of complex images allows for a much more compact (typically ~5 terms) approximate Green's function. A complex image has both a complex position and a complex amplitude. Since static fields are real, the complex images of a point source always occur in complex conjugate pairs. Complex images are computed numerically by approximating the Green's function in the spectral domain with complex exponentials using Prony's method. The spectral domain Green's function is computed using the transmission line method outlined in [12]. 6.2.2 Multipipe Model The multipipe model [7] is used for the determination of an appropriate expansion of p(r'). Each strip is represented by a series of pipes each with constant unknown charge. Gauss-Chebyshev quadrature is used to approximate the Green's integral equation, (50), by a finite summation which essentially results in the replacement of the unknown continuous charge distribution by a finite series of point charges with unknown values. Then, applying the results of Wheeler's work in which he derives a single pipe model of a strip in free space assigns a radius to each of the individual point charge locations, and hence pipes are formed . The result is a series of pipes with different radii and unknown charges that takes the place of the original strip. The charges on each pipe are related to the charge distribution on the strip by a simple relationship. With forced potentials on the strips and the appropriate Green's function, the charges are found using the M o M . The moment matrix resulting from such an approach is very compact for high computational efficiency. 6.2.3 Results Extending the problem of Example 2, the source, gate, and drain of the MESFET are fed by strips numbered 1,2, and 3, respectively. Strip 1 and Strip 3 are each 10 um wide, and Strip 2 is 0.4 um wide. In applying the M o M , thirteen pipes are used to model each of Strip 1 and Strip 61 3, while six are used for Strip 2 yielding a 32x32 moment matrix. Five complex image terms are sufficient to accurately approximate the Green's function of the air-substrate-ground plane system. Using the procedure outlined in Appendix B , the capacitances of the strips are calculated to be C l l = C 3 3 = 2 3 0 - ° P F / m C 2 2 = 148.7 pF/m *-T2 = ^23 = 72-9 pF/m C 1 3 = 98.2 pF/m Superposition of these parasitic capacitances with those computed in Example 2 gives the modified circuit model shown in Figure 25. The extra capacitance has the overall effect of slowing x = transit time (~ 10 ps) R Q and Rrj = parasitic resistances g m = 13 mho/m ggS =16.2 mho/m g,js = 2.1 mho/m C ' G G = 437.6 pF/m C ' D D = 330.0 pF/m C ' D S = 277.4pF/m C ' G S = 353.4 pF/m § C ' G D = 95.6 pF/m Figure 25. Modified small signal model of MESFET in Example 2. the device down. Figure 26 compares magnitude versus frequency plots of the transfer functions for both models. The analytic expression for the transfer function of the circuit model is 62 V - J O J T g d s + MCDD + cDS) -torn w 2 C G D ( C G G + C G S ) - ja»CG Dg gs ggs + J « ( C G D + C G G + C G S ) As can be seen in Figure 26, inclusion of the parasitic capacitances in the modified circuit model causes the 3 dB cutoff frequency to decrease by more than one half, from 7 GHz to 3 Ghz. This illustrates the level of importance of parasitic capacitances in the performance of a MESFET. ~i 1 1 1 1 1 i r~ 3 GHz Modified Model - Original Model 7 GHz 0 2 4 6 8 10 12 14 16 18 20 Frequency (GHz) -—> Figure 26. Gain versus frequency plots for circuit models. 6.3 Artificial Boundary Condition Adjustments In both Example 1 and Example 2, the potential is set to zero on the artificial boundary in order to emulate the effects of the ground plane of the substrate which is actually much farther away. The thickness of the substrate is typically on the order of 100 \im making it impractical to include it in the solution domain of the F E M problem due to the many extra nodes that would be required in such a situation. In addition to this, the zero diffusion and drift current conditions on the left and right boundaries are not entirely realistic since some current does, in reality, travel outside of the active region of the device and into and out of the portions of the electrode strips that 63 extend beyond the active region as depicted in Figure 24. Following, possible solutions to these deficiencies using integral equation techniques are proposed. 6.3.1 Green's Integral Equation A simple approach is to modify the F E M solution domain by extending the buffer zone around the entire active region. Zero net current boundary conditions are imposed on the boundary around the buffer zone. Using the Green's integral equation and M o M as explained earlier, the charges on the strip conductors are calculated. From these charges the field along the artificial boundary is computed and is used as a potential boundary condition in the F E M model. Applying this procedure to the MESFET in Example 2, the buffer zone is made to surround the active region and is 1 urn thick. The mesh is shown in Figure 27. The potential and carrier density fields computed with F E M and the new artificial boundary conditions for V d s = 4 V and V g S = 0 V are shown in Figure 29. Drain current versus drain-source voltage for V g S = 0 is plotted in Figure 28. There is considerable difference between this current plot and the one in Figure 20 even though the device being modelled is the same. Figure 27. Modified mesh grid for Example 2 - 1677 nodes, 800 elements. 64 Drain current versus drain-to-source voltage (Vgs = 0 V). Vds (V) —> Figure 28. Drain current versus drain-source voltage. These results are questionable and still under investigation due to the following reasons. Recall (See Section 2.) that ohmic contact conditions were applied at the boundaries between the source and drain conductors, and the active region of the MESFET in place of the heavily doped wells that exist in reality. This simplification greatly reduced the number of necessary nodes in the model (See Section 5.1). Unfortunately in this case, this simplification introduces a large step discontinuity in the carrier concentration condition along the source and drain portions of the boundary. Specifically, the portion of the boundary between the active region and the source and drain conductors has a carrier density condition of 2x10^ cm"^, while the portion of the boundary between the buffer region and the source and drain conductors has a carrier density condition of 1.621x10^ cm~3 due to the Schottky barrier. This discontinuity results in unrealistic, nonphysical current behavior in the vicinity of the step. Of course, if the ohmic contact approximations would be replaced by the actual heavily doped contact wells, this abrupt step discontinuity does not occur, but the number of nodes required as a result of this would be costly. This scenario has not been investigated here. Another potential problem is that the artificial carrier walls that appear along most of the buffer zone boundary are quite high (See Figure 29.), and it is questionable that they have no effect on the true carrier and potential distributions in and around the active region. 65 o ; = 4 V or 3.14 x 1 0 1 6 c m " 3 = -0.7211 V o r 16210 cm" 3 Figure 29. Artificial boundary adjustment results. 66 6.3.2 Boundary Integral Equation Instead of using the Green's integral equation to find the potential along the artificial boundary, the above approach could be further improved by incorporating an M o M matrix of a boundary integral equation (BIE) directly into the original F E M formulation thereby coupling the interior of the F E M domain to the exterior. The BIE expresses the relationship between the potential at a point on the boundary and the potential on the entire boundary. As outlined in [9], the BIE that satisfies the Laplace equation is where G(r,r') is the Green's function satisfying the exterior boundary conditions, (j)(r) is the potential along the boundary T, and a(r') = (V(()(r')-n), the potential flux normal to F. Using the M o M to discretize this equation , a matrix expression of a at a finite number of points in terms of the potential at those same points is derived: where [(j)] and [a] are nodal values of (j) and a, respectively, along T. An approximate expression for a is obtained by adding some shape functions [M] to (52), so a « [M] [a]. Next, although [9] uses a variational formulation to complete the derivation, it is shown here how the familiar Galerkin's method can also be used to arrive at a similar final result. Substituting the approximation for a and (52) into the surface term in (10) (global version) results in a new expression for this term: J r[N]TeV[N][(t)]-ndS = - J r [ N ] T e a dS = - ( f r [ N ] T £ [M]dS [ B ] " 1 ^ ] ) [<|>] . (53) (|>(r)/e = J r(()(r')(VG(r,r')-n)dr' - / r G ( r , r ' ) a(r') dr' (51) [A][<|>] = [B][a] [a] = [BYl[A] [(M (52) 67 The minus sign is introduced since the normal vector of the interior F E M domain points in the opposite direction to the normal vector of the exterior BIE domain. The term enclosed in the brackets is a matrix, and it can be thought of as the stiffness matrix for a "boundary element" comprised of and coupling all the nodes on the outer boundary of the F E M problem. This stiffness matrix differs from the one resulting from the variational approach in that it is not symmetric. The boundary element is treated as any other finite element, and the components of its stiffness matrix are added to the global stiffness matrix in the same manner. This hybrid method has proven effective in open electromagnetic scattering problems [40], and its application to the MESFET problem is currently being studied. 68 7 Conclusion This paper presented a two-dimensional finite element analysis of the stationary semiconductor device equations. Exploiting the best features of two different F E M formulations, the current conservation model and the advection-diffusion model, along with two different numerical solutions, Gummel's method and the Newton Raphson method, a robust and efficient solution approach was developed. This technique was put to the test in two sample problems in which potential and carrier density distributions were computed inside GaAs MESFETs. In order to demonstrate the usefulness of these field distributions for the purpose of circuit analysis, various MESFET parameters such as electrode currents, voltage gain, capacitances, and conductances were calculated and incorporated into a small-signal circuit model of one of these devices. Implementation considerations, such as stability and accuracy restrictions, of the F E M solution were also addressed. Finally, various extensions to the F E M approach involving the application of the method of moments (MoM) were briefly discussed and partially demonstrated. These extensions consist of Green's integral and boundary integral equation solutions that are implemented in such a way as to compensate for the limitations of a finite domain and to improve physical accuracy of the F E M solution. The ultimate goal of such extensions is to couple the finite domain to the rest of the circuit of which it is a sub-region. 7.1 Accomplishments The ability to model the steady-state fields inside semiconductor materials is an important step in properly accounting for active regions in MMICs. Through the development of a new formulation, it has been shown that the finite element method can be used reliably in the solution of the classical semiconductor equations at least in the case of majority carrier devices, which are most prominent in MMIC designs. Additionally, it has also been suggested that the choice of a finite element approach over a finite difference scheme has its advantages in light of the "big picture" of 69 MMIC modeling. The similarity in the mathematical structures of the finite element method and the method of moments together with the inherent flexibility of boundary handling in the finite element method offers a direct means for the coupling of the two methods. This is necessary if the mutual field coupling among all the components of an MMIC is to be considered. Failure to account for these component interactions may lead to inaccurate performance predictions in the overall design of an M M I C . 7.2 Future Work Application of the F E M formulations presented in this paper is not restricted to majority-carrier devices and can easily be applied in the simulation of other semiconductor devices. Also, these F E M formulations, as they stand, are not restricted to just two dimensions and, theoretically, can be applied to three dimensional problems as well, without any modification of the expressions or notation shown here. This would be possible by simply using three-dimensional finite elements in place of the two-dimensional elements implemented here. These two points open the door to a substantial amount of further investigation of and testing for the effectiveness of this finite element analysis in dealing with a myriad of semiconductor problems. Also, the application of boundary integral equations to account for exterior effects on the finite element domain has yet to be accomplished. Hybrid boundary integral finite element methods have successfully been applied to linear open electromagnetic scattering problems and open static field problems. Use of a boundary integral equation discretized with the method of moments leads to a convenient boundary element which is incorporated into the finite element formulation as any other element. For linear problems, this leads to one matrix equation that can be solved simultaneously; iterative techniques are not required. Herein lies the advantage of this hybrid method, as opposed to other exisiting iterative hybrid schemes, for dealing with nonlinear problems (nonlinear interior, linear exterior), since any other iterations besides the already necessary iterations due to the nonlinear aspect of the problem, e.g. Gummel or Newton-Raphson iterations, would prove to be quite costly in terms of speed. The difficulties with the boundary 70 integral approach are rooted in the challenge of accurately evaluating integrals over singularities (self-terms) and chosing appropriate basis functions both of which are necessary for an accurate moment method formulation. There is much room for experimentation in this area. 71 Bibliography [1] D. M . Pozar, Microwave Engineering, Addison-Wesley Publishing Company Inc., 1990. [2] P. A . Markowich, The Stationary Semiconductor Device Equations, New York: Springer-Verlag, 1986. [3] J. Svacina, Analysis of Multilayer Microstrip Lines by a Conformal Mapping Method, IEEE Transactions on Microwave Theory and Techniques, Vo l . 40, No. 4, pp. 769-772, April 1992. [4] C. Wei, R. F. Harrington, J. R. Mautz, and T. K . Sarkar, Multiconductor Transmission Lines in Multilayered Dielectric Media, IEEE Transactions on Microwave Theory and Techniques, Vol . 32, No. 4, pp. 439-449, April 1984. [5] R. F. Harrington, Field Computaion by Moment Methods, New York: Macmillan, 1968. [6] N . G. Alexopolous, J. A . Maupin, and P. T. Greiling, Determination of the Electrode Capacitance Matrix for GaAs Fet's, IEEE Transactions on Microwave Theory and Techniques, Vo l . MTT-28, pp. 459-466, May 1980. [7] G. E. Howard, J . J . Yang, and Y . L . Chow, A Multipipe Model of General Strip Transmission Lines for Rapid Convergence of Integral Equation Singularities , IEEE Transactions on Microwave Theory and Techniques, Vol . 40, No. 4, pp. 628-636, April 1992. 72 [8] K . S. Oh, D. Kuznetsov, and J. E. Schutt-Aine, Capacitance Computations in a Multilayered Dielectric Medium Using Closed-Form Spatial Green's Functions, IEEE Transactions on Microwave Theory and Techniques, Vol . 42, No. 8, pp. 1443-1453, August 1994. [9] O. C. Zienkiewicz,77ie Finite Element Method, London: McGraw-Hill, 1977. [10] G. R. Buchanan, Schaum's Outline of Theory and Problems of Finite Element Analysis, New York: McGraw-Hill, 1995. [11] D. S. Burnett, Finite Element Analysis, Reading, Massachusetts: Addison-Wesley, 1987. [12] R. Crampagne, M . Ahmadpanah, and J. Guiraud, A Simple Method for Determining the Green's Function for a Large Class of MIC Lines Having Multilayered Dielectric Strucures, IEEE Transactions on Microwave Theory and Techniques, Vol . MTT-26, No. 2, pp. 82-87, February 1978. [13] Z. Pantic and R. Mittra, Quasi-TEM Analysis of Microwave Transmission Lines by the Finite-Element Method, IEEE Transactions on Microwave Theory and Techniques, Vol . MTT-34,No. 11, pp. 1096-1103, November 1986. [14] A . Khebir, A . B. Kouki, and R. Mittra, Higher Order Asymptotic Boundary Condition for the Finite Element Modeling of Two-Dimensional Transmission Line Structures, IEEE Transactions on Microwave Theory and Techniques, Vol . 38, No. 10, pp. 1433-1437, October 1990. 73 [15] H . K . Gummel, A Self-consistent Iterative Scheme for One-dimensional Steady State Transistor Calculations, IEEE Transactions on Electron Devices, Vol . Ed-11, pp. 455-465, 1964. [16] A . de Mari, An Accurate Numerical One-Dimensional Solution of the PN-Junction Under Arbitrary Transient Conditions, Solid State Electronics, Vol . 11, pp. 1021-1053, 1968. [17] M . Adler, A Method for Achieving and Choosing Variable Density Grids in Finite Difference Formulations and the Importance of Degeneracy and Band Gap Narrowing in Device Modeling, Numerical Analysis of Semiconductor Devices, Dublin: Boole Press, 1979. [18] A . F. Franz, G. A . Franz, S. Selberherr, C. Ringhofer, and P. Markowich, Finite Boxes -A Generalization of the Finite Difference Method Suitable for Semiconductor Device Simulation, IEEE Transactions, Vol . ED-30, pp. 1070-1082, 1983. [19] M . S. Mock, Basic Theory of Stationary Numerical Models, Process and Device Modeling, North-Holland: Elsevier Science Publishers B . V . , pp. 159-195,1986. [20] R. E . Bank., W. M . Coughran, Jr., W. Fichtner, D. J. Rose, and R. K . Smith, Computational Aspects of Semiconductor Device Simulation, Process and Device Modeling, North-Holland: Elsevier Science Publishers B . V . , pp. 229-265,1986. [21] J .J . Barnes and R . J . Lomax, Two-Dimensional Finite Element Simulation of Semiconductor Devices, Electron. Letters, Vol . 10, pp. 341-343, August 8, 1974. 74 [22] J. J. Barnes, R. J. Lomax, and G . I . Haddad, Finite-Element Simulation of GaAs MESFET's with Lateral Doping Profiles and Submicron Gates, IEEE Transactions on Electron Devices, Vol . Ed-23, pp. 1042-1048, September 1976. [23] J .J . Barnes and R . J . Lomax, Finite-Element Methods in Semiconductor Device Simulation, IEEE Transactions on Electron Devices, Vol . Ed-24, pp. 1082-1089, August 1977. [24] P. E. Cottrell and E. M . Buturla, Steady State Analysis of Field Effect Transistors via the Finite Element Method, International Electron Devices Meeting Technical Digest, Washington, DC, pp. 51-54, December 1975. [25] G . L . Tan, X . L . Yuan, Q. M . Zhang, W. H . K u , and A . J. Shey, Two-Dimensional Semiconductor Device Analysis Based on New Finite-Element Discretization Employing the S-G Scheme, IEEE Transactions on Computer-Aided Design, Vo l . 8. No. 5, pp. 468-478, May 1989. [26] D. L . Pulfrey and N . G. Tarr, Introduction to Microelectronic Devices, Englewood Cliffs, New Jersey: Prentice Hall, 1989. [27] A . N . Brooks and T. J. R. Hughes, Streamline Upwind!Petrov-Galerkin Formulations for Convection Dominated Flows with Particular Emphasis on the Incompressible Navier-Stokes Equations, Computer Methods in Applied Mechanics and Engineering , Vol . 32, pp. 199-259, 1982. [28] C Froberg, Introduction to Numerical Analysis, Reading, Massachusetts: Addison-Wesley, 1965. 75 [29] P. M . Gresho and R. L . Lee, Don't Suppress the Wiggles - They're Telling You Something! , Computers and Fluids, Vol . 9, pp. 223-253, Great Britain: Pergamon Press Ltd., 1981. [30] S. V . Marshall and G. G. Skitek, Electromagnetic Concepts and Applications, Englewood Cliffs, New Jersey: Prentice Hall, 1990. [31] O. C. Zienkiewicz and D. V . Phillips, An Automatic Mesh Generation Scheme for Plane and Curved Surfaces by 'Isoparametric' Co-ordinates, International Journal for Numerical Methods in Engineering , Vol . 3, pp. 519-528, 1971. [32] W. R. Curtice, A MESFET Model for Use in the Design of GaAs Integrated Circuits, IEEE Transactions on Microwave Theory and Techniques, Vol . MTT-28, pp. 448-456, May 1980. [33] S. Y . Liao, Microwave Circuits and Devices, Prentice Hall Inc., 1980. [34] T. Takada, K . Yokoyama, M . Ida, and T. Sudo, A MESFET Variable-Capacitance Model for GaAs Integrated Circuit Simulation, IEEE Transactions on Microwave Theory and Techniques, Vol . MTT-30, pp. 719-723,1982. [35] R. L . Van Tuyl, C. A . Liechti, and C. A . Stolte, Gallium Arsenide Digital Integrated Circuits, Air Force Avionics Lab., A F W A L , WPAFB, Technical Report AFAL-TR-26-264, April 1977. [36] J. Maupin, P. Greiling, and N . Alexopoulos, Speed Power Tradeoff in GaAs FET Integrated Circuits, paper presented at 1st Specialty Conference Gigabit Logic for Microwave Systems, Orlando, Florida, May 1979. 76 [37] R. L . Van Tuyl, C. A . Liechti, R. E . Lee, and E. Gowen, GaAs MESFET Logic with 4-GHz Clock Rate, IEEE Journal on Solid State Circuits, Vol . SC-12, pp. 485-496, October 1977. [38] W. E. Boyce and R. C. DiPrima, Elementary Differential Equations and Boundary Value Problems, John Wiley & Sons Inc., 1986. [39] Y . L . Chow, J. J. Yang, and G. E. Howard, Complex Images for Electrosatatic Field Computation in Multilayered Media , IEEE Transactions on Microwave Theory and Techniques, Vol . 39. No. 7, pp. 1120-1125, July 1991. [40] P .P . Silvester and G. Pelosi, Finite Elements for Wave Electromagnetics Methods and Techniques, New York: IEEE Press, 1994. [41] W. R. Smythe, Static and Dynamic Electricity, McGraw-Hill, 1968. 77 Appendix A Derivation of Shape Functions The following derivation of shape functions is an extension of Problem 3.2 in [10]. A function f in a two-dimensional space can be approximated by a polynomial of the form f(x,y) = C ^ C x . y ) + C 2 g 2 (x ,y) + . . . + C M g M (x ,y ) (A. l) where gj(x,y) is some arbitrary known function of space, and Cj is some arbitrary constant. Given values of this function at M specific points in space gives a set of M equations and M unknowns: fCxj.yj) = C 1 g 1 ( x 1 , y 1 ) + C 2 g 2 (x 1 , y 1 ) + . . . + C ^ ^ x ^ ) f ( x 2 , y 2 ) = C 1 g 1 ( x 2 , y 2 ) + C 2 g 2 ( x 2 , y 2 ) + . . . + C M g M ( x 2 , y 2 ) f ( x M > y M ) = c i g i ( x M - y M ) + c 282 ( x M>yM) + • • • + c M g M (x M , y M ) or m = rate] where fj = f(xj,yj) and Xjj = gj(xj,yj) . So, the solution for the coefficients is given by [C] = [X] _ 1 [f]. Now, (A.l) can be rewritten as f(x,y) - [g][C] where [g] is a row vector with gj = g^(x,y). Substituting for [C] in this expression gives 78 f(x,y) ~ [ g l t X ] - 1 ^ =[N][f] =N 1 (x ,y)f (x 1 ,y 1 )+ N ^ x . y ) ^ ^ ) + . . . + N M (x ,y) f (x M ,y M ) (A.2) where Nj (x,y) is the shape function for point i . So, (A.2) expresses the original function in terms of it's values at a finite number of points, M , and the shape functions at these points. Note that this procedure is not limited to two dimensions and is just as easily applicable to higher dimensions. As a simple example, an expression for the shape functions of a linear triangular element is derived. Here , the three points, or nodes, at which the function to be approximated is known are located at the three vertices of the triangle, (x^,y^), (X2,y2)> and ^ 3 ^ 3 ) . Next, [g] is chosen to be [1 x y]T. This is the same as saying that gj(x,y) = 1, g2(x,y) = x, and g 3(x,y) = y in ( A . l ) . The matrix [X] is constructed: [X] 1 J 1 It follows that [N] = [1 x y] 1 x, 1 x 2 1 X , (A .3) Solving for an analytic expression for Nj (x,y), requires some considerable amount of algebra. Analytic solutions become much more difficult to obtain for higher order elements. Nevertheless, in practice there is no need to have analytic expressions for the shape functions, and (A.3) is ideally suited for numerical applications. Once [X]" ' has been computed, it is a simple matter to calculate the value of the shape functions at any point in space simply by plugging the coordinates 7 9 into [g] and solving (A3). This method for computing shape functions can be applied to a space of any number of dimensions. On the other hand, choosing an appropriate set of functions [g] and sampling points is quite arbitrary. For triangular elements of any order, Pascal's triangle can be used to find the proper [g]. 80 Appendix B Capacitance Calculations for a Multi-Port System A circuit model composed of N nodes interconnected by capacitors representing the capacitive coupling for a system of N conductors is sought. It is assumed that a means is available for calculating charges, Q n , on the conductors given any arbitrary set of potentials, V n , forced on the conductors. The charge on conductor 1 is defined as Q l = C 1 1 V 1 + C 1 2 V 2 + • • + C m V N Similarly, for the charge on conductor 2 Q 2 = C 2 1 V ! + C 2 2 V 2 + . . . + C 2 N V N and so forth for the charges on the remaining N-2 conductors. Putting these relationships between the conductor voltages and charges into matrix form gives [Q] = [C][V] where [Q] is an N x l column vector with components Q n , [V] is an N x l column vector with components V n , and [C] is an NxN matrix with components C n m . [C] is referred to as the capacitance matrix [41]. By definition, the N components of the first column of [C] are simply the values of the charges on the N conductors for the case of a potential of value one on the first conductor and zero on the rest. The N components of the second column of [C] are the values of the charges on the N conductors for the case of a potential of value one on the second conductor and zero on the rest, and so on. The magnitude of each individual component of [C], I C n m l , is the value of the circuit capacitor between node n and m for n * m, mutual capacitance, or the value of the capacitor from n to ground for n = m, self-capacitance, in the equivalent circuit model. For example, two microstrip conductors are represented by the equivalent circuit model as shown in Figure B1. IC icni 4= 12' 2 I 9-t I C 2 2 I Figure B1. Circuit model for two strip conductors. It should be noted that for a passive system, C n m = C m n , due to reciprocity
- Library Home /
- Search Collections /
- Open Collections /
- Browse Collections /
- UBC Theses and Dissertations /
- A two-dimensional finite element analysis of the stationary...
Open Collections
UBC Theses and Dissertations
Featured Collection
UBC Theses and Dissertations
A two-dimensional finite element analysis of the stationary semiconductor device equations Chavez, Patrick Pablo 1997
pdf
Page Metadata
Item Metadata
Title | A two-dimensional finite element analysis of the stationary semiconductor device equations |
Creator |
Chavez, Patrick Pablo |
Date Issued | 1997 |
Description | The ability to model the steady-state field inside active structures, such as a transistor, is an important aspect of monolithic microwave integrated circuit (MMIC) design. This paper focuses on such an active zone of semiconductor material, and presents a finite element analysis of the classical semiconductor equations. The semiconductor equations are very nonlinear and govern the potential and carrier density distributions in semiconductor materials. A previously developed finite element method (FEM) formulation of these equations, referred to as the current conservation model, is re-introduced and re-derived with compact matrix notation. It is shown how this formulation can be solved with the Newton-Raphson iterative scheme. Then, a newly developed FEM formulation, referred to as the advection-diffusion model, of the continuity equations is described in detail. It is shown by example how this formulation solved with Gummel's iterative technique is very numerically robust. These two different solution methods of the steady-state system of coupled Poisson and continuity equations are combined into a final solution algorithm that exploits their strengths. As a specific example, GaAs MESFETs are the focus of implementation, and the resulting potential field and carrier density distributions are used, to calculate various MESFET parameters such as electrode currents, voltage gain, capacitances, and conductances. Finally, various extensions to the FEM approach involving the application of the method of moments (MoM) are briefly discussed and partially demonstrated. These extensions are intended to compensate for the assumptions and simplifications, mainly with respect to the artificial boundary conditions, used in the original stand-alone FEM formulations. |
Extent | 5902343 bytes |
Genre |
Thesis/Dissertation |
Type |
Text |
FileFormat | application/pdf |
Language | eng |
Date Available | 2009-03-24 |
Provider | Vancouver : University of British Columbia Library |
Rights | For non-commercial purposes only, such as research, private study and education. Additional conditions apply, see Terms of Use https://open.library.ubc.ca/terms_of_use. |
DOI | 10.14288/1.0065230 |
URI | http://hdl.handle.net/2429/6387 |
Degree |
Master of Applied Science - MASc |
Program |
Electrical and Computer Engineering |
Affiliation |
Applied Science, Faculty of Electrical and Computer Engineering, Department of |
Degree Grantor | University of British Columbia |
GraduationDate | 1997-11 |
Campus |
UBCV |
Scholarly Level | Graduate |
AggregatedSourceRepository | DSpace |
Download
- Media
- 831-ubc_1997-0383.pdf [ 5.63MB ]
- Metadata
- JSON: 831-1.0065230.json
- JSON-LD: 831-1.0065230-ld.json
- RDF/XML (Pretty): 831-1.0065230-rdf.xml
- RDF/JSON: 831-1.0065230-rdf.json
- Turtle: 831-1.0065230-turtle.txt
- N-Triples: 831-1.0065230-rdf-ntriples.txt
- Original Record: 831-1.0065230-source.json
- Full Text
- 831-1.0065230-fulltext.txt
- Citation
- 831-1.0065230.ris
Full Text
Cite
Citation Scheme:
Usage Statistics
Share
Embed
Customize your widget with the following options, then copy and paste the code below into the HTML
of your page to embed this item in your website.
<div id="ubcOpenCollectionsWidgetDisplay">
<script id="ubcOpenCollectionsWidget"
src="{[{embed.src}]}"
data-item="{[{embed.item}]}"
data-collection="{[{embed.collection}]}"
data-metadata="{[{embed.showMetadata}]}"
data-width="{[{embed.width}]}"
async >
</script>
</div>
Our image viewer uses the IIIF 2.0 standard.
To load this item in other compatible viewers, use this url:
http://iiif.library.ubc.ca/presentation/dsp.831.1-0065230/manifest