Citation |

- Permanent Link:
- http://digital.auraria.edu/AA00003824/00001
## Material Information- Title:
- A PN semiconductor device simulator using MATLAB
- Creator:
- Pace, Shawn Randall
- Publication Date:
- 2011
- Language:
- English
- Physical Description:
- ix, 114 leaves : illustrations ; 28 cm
## Subjects- Subjects / Keywords:
- MATLAB ( lcsh )
MATLAB ( fast ) Semiconductors -- Junctions ( lcsh ) Numerical analysis -- Computer programs ( lcsh ) Numerical analysis -- Computer programs ( fast ) Semiconductors -- Junctions ( fast ) - Genre:
- bibliography ( marcgt )
theses ( marcgt ) non-fiction ( marcgt )
## Notes- Bibliography:
- Includes bibliographical references (leaf 114).
- General Note:
- Department of Electrical Engineering
- Statement of Responsibility:
- by Shawn Randall Pace.
## Record Information- Source Institution:
- |University of Colorado Denver
- Holding Location:
- |Auraria Library
- Rights Management:
- All applicable rights reserved by the source institution and holding location.
- Resource Identifier:
- 781635746 ( OCLC )
ocn781635746 - Classification:
- LD1193.E54 2011M R34 ( lcc )
## Auraria Membership |

Downloads |

## This item has the following downloads: |

Full Text |

A PN SEMICONDUCTOR DEVICE SIMULATOR USING MATLAB
by Shawn Randall Pace B.S.E.E., University of North Florida, 2007 A thesis submitted to the University of Colorado Denver in partial fulfillment for the degree of Master of Science Electrical Engineering 2011 This thesis for the Master of Science Electrical Engineering degree by Shawn Randall Pace has been approved by Hamid Fardi Robert Grabbe /!//?-/ II Date Pace, Shawn Randall (M.S., Electrical Engineering) A PN Semiconductor Device Simulator Using MATLAB Thesis directed by Professor Hamid Fardi ABSTRACT The purpose of this project is to develop a functional PN semiconductor device simulator that is modular in nature in order to allow for flexibility during programming and to allow for future modification with relative ease. The device modeling program is developed using MATLAB In addition, the programs main goal is to provide a simulator tool that can supplement device modeling, semiconductor device physics, and numerical analysis to construct basic semiconductor devices which can be studied using standard numerical analysis techniques. MATLABs capability and inherent nature of handling matrices and matrix operations makes this approach an excellent technique to develop numerical analysis algorithms. The program solution will be used to examine device parameters such as carrier statistics, device potential, and internal electric fields. The device solution is compared to analytical approximations in order to further strengthen the understanding between theory and exact numerical solutions and how those solutions are obtained. This abstract accurately represents the content of the candidates thesis. I recommend its publication. Signed Hamid Fardi ACKNOWLEDGEMENT The main credit and acknowledgement for the development of this project goes to Dr. Hamid Fardi, without whose knowledge and support would have made the creation of this project nearly impossible. His knowledge of numerical analysis related to semiconductor physics and device physics operation as a whole was immensely appreciated and invaluable. Furthermore, appreciation should be given to the University of Colorado Denver, whose Electrical Engineering department and faculty provided the necessary training and knowledge over the course of my work there to make this project possible. It is my sincere hope that this tool can be used in a valuable training manner and supplemental learning aid. TABLE OF CONTENTS Figures.....................................................................vii Tables.......................................................................ix Chapter 1. Introduction...........................................................1 2. Semiconductor Device Physics...........................................4 2.1 Semiconductor Material Parameters......................................4 2.2 Device Equations.......................................................5 2.3 Carrier Densities.....................................................12 2.4 Mobility Models.......................................................13 2.5 Recombination and Generation Models...................................15 3. Numerical Analysis....................................................18 3.1 Taylor Series Expansion...............................................18 3.2 Numerical Differentiation.............................................21 3.3 Numerical Integration.................................................23 3.4 Tri-diagonal Matrices.................................................25 3.5 PN Junction DC Analysis...............................................28 3.6 Numerical Solution of PN Junction.....................................51 4. SDS Results Investigated..............................................56 4.1 SDS Output Compared to Analytical Approximations......................56 v 5. Semiconductor Device Simulator....................................67 6. Summary............................................................73 Appendix A. Integral Form of the Current Density Equation......................74 B. Derivation Techniques of the Matrix-vector Coefficients............81 C. Semiconductor Device Simulator Implemented........................85 References...............................................................114 vi LIST OF FIGURES Figure 2.1 INFINITESIMAL SLICE FOR DETERMINATION OF INCREASE IN CARRIER DENSITY IN X..................6 2.2 ELECTRON CONCENTRATION DIFFUSION AFFECT...........9 3.1 INTEGRATION USING THE TRAPEZOIDAL RULE............23 3.2 TRI-DIAGONAL MATRIX...............................25 3.3 BASIC P-N JUNCTION DIODE..........................29 3.4 DIVISION POINTS FOR DC ANALYSIS...................32 3.5 MATRIX FORM OF THE MATRIX-VECTOR EQUATION.........40 3.6 COEFFICIENT A OF MATRIX-VECTOR EQUATION...........41 3.7 COEFFICIENT B OF MATRIX-VECTOR EQUATION...........41 3.8 COEFFICIENT C OF MATRIX-VECTOR EQUATION...........42 3.9 COEFFICIENT F OF MATRIX-VECTOR EQUATION...........42 3.10 UPDATED BIAS POTENTIAL BASED ON PREVIOUS POTENTIAL.........................................54 4.1 ABRUPT STEP JUNCTION PROPERTIES...................59 4.2 EQUILIBRIUM DEVICE POTENTIAL......................60 4.3 EQUILIBRIUM ELECTRIC FIELD........................61 4.4 DEVICE POTENTIAL WITH UNEQUAL DOPING DENSITIES....62 vii 4.5 SIMWINDOWS COMPARISON TO FIGURE 4.4............62 4.6 SDS OUTPUT SHIFTED TO ZERO.....................63 4.7 SIMWINDOWS ELECTRIC FIELD......................63 4.8 SDS ELECTRIC FIELD.............................64 4.9 SDS DEVICE POTENTIAL WITH EXTERNAL APPLIED BIAS...................................65 4.10 SIMWINDOWS POTENTIAL WITH EXTERNAL APPLIED BIAS...................................66 5.1 SDS BASIC FLOWCHART............................72 A. 1 INTEGRAL FORM OF CURRENT DENSITY EQUATIONS.....79 viii LIST OF TABLES Table 2.1 Material parameters for silicon...........................................4 2.2 Mobility parameters......................................................14 2.3 Ionization generation parameters.........................................17 3.1 Partial derivatives required for matrix coefficients.....................43 3.2 Numerical constants......................................................54 3.3 Physical constants.......................................................55 IX 1. Introduction Although there are several powerful industry and student level device physics simulators already available, the hope is that with this new tool and supplemental report, the theory of numerical analysis applied to device physics can be studied at a more fundamental level. Traditional simulators allow for device simulation but do not offer any insight into how that solution was obtained. The addition of the report that supplements the Semiconductor Device Simulator will develop the basic tools necessary to understand the operation of the program and allow future modifications as necessary. This project was developed in the programming environment provided by MATLAB which is developed and supported by Mathworks. Almost all academic institutions use MATLAB heavily in the learning process and it is readily available to most students, hence the choice of this programming environment. Additionally, MATLAB is an extremely powerful numerical analysis tool that makes the solution of a program like this one rather straightforward. MATLABs capability and inherent nature of handling matrices and matrix operations makes this an excellent tool to develop numerical analysis algorithms. The current semiconductor simulation tools available do not lend to updates or modifications in a simple or straightforward manner. In fact most tools are not open to be changed in any manner unless the source code is available, which is rare. With the open nature of MATLAB the solution process can be studied throughout and all variables and parameters are available for examination. The program could have the capability to be standalone which allows for the use of the tool without the requirement of the MATLAB environment.* The symbol used throughout this manual is to indicate a future development not yet built into the current version of the program. Almost all sections in this manual include a future development section to indicate where the code can be easily updated to include new functions or models as necessary. 1 The remaining discussion in this section gives a brief overview of the contents included in the following sections. Chapter 2 gives a brief introduction to the devices physics required to build a numerical analysis program to simulate semiconductor devices in general. Device physics equations (continuity, Poissons, current, etc.) are discussed as are the models, such as mobility, SRH (Shockley-Read-Hall) recombination, and generation. Boundary conditions are also discussed. Lastly, silicon material parameters are given and briefly discussed. Currently no other semiconductor material is used in this program, but the addition of other materials is rather straightforward*. Chapter 3 is broken into several sections. The first section will give a brief introduction to numerical analysis techniques such as differentiation using the difference method, integration using the trapezoidal technique, Taylor series expansion, and basic tri-diagonal matrix solutions. The next section will discuss the numerical parameters required to support the program, such as mesh spacing, iteration limits, convergence tolerances, etc. The final section will give an in-depth view of setting up and solving the numerical equations for a basic PN junction device. The section will cover converting the basic device physics equations in Chapter 2 into numerical equations, setting up the tri-diagonal block matrix, performing the recursive solution algorithm and determining convergence criteria. Chapter 4 will provide simulation results and compare them to basic theory (depletion width, built-in potential, etc.) as well as with other Simulation tools, such as SimWindows. Chapter 5 will discuss the program as written and the various functions built to support it. Each function will be discussed and explained. It should be noted however, that the MATLAB code itself is heavily commented as a supplement to this section. Chapter 6 gives a brief summary of the report. 2 Appendix A Integral Form of the Current Density Equation Appendix B Derivation Techniques of the Matrix-vector Coefficients Appendix C Semiconductor Device Simulator Implemented A final note on the Introduction is in order. Currently this program supports basic 1-D PN junction devices with uniform mesh spacing. Although the program is hardcoded for a Step Junction, the junction profile is easily changed just as the numerical parameters (convergence tolerance, relaxation factor, device length and doping) are without any further modification to the program. Modifications for 2-D solutions, refined mesh spacing, Shottky diodes, and BJTs could be further implemented*. The same can be said for different numerical solution techniques if desired. 3 2. Semiconductor Device Physics In this chapter, a brief discussion of simple device physics will be covered along with the differential equations governing them as well as any necessary models required to establish a numerical solution of the differential equations. Semiconductor physical parameters will be given in Section 2.1, device equations discussed in Section 2.2, carrier densities in Section 2.3, mobility models in Section 2.4 and recombination- generation models in Section 2.5. 2.1. Semiconductor Material Parameters A tabulated list of the semiconductor physical parameters necessary to reach a numerical solution and some extra parameters used for future expansions are given below. Currently only parameters for Silicon material are used. The corresponding program variable is also shown. Table 2.1 Material parameters for silicon Nomenclature Program Variable Default Value Units Relative Permittivity er 11.7 - Silicon Permittivity es 1.036 x 10-12 F/cm Energy Gap (300K) EG300 1.08 eV Alpha EG_alpha 4.73 x 10-4 - Beta EG_beta 636 - Conduction Band Density (300K) NC300 2.8 x 1019 cm-3 Valence Band Density (300K) NV300 1.04 x 1019 cm-3 Electron Mobility Function* 1350 cm2/V-sec Hole Mobility Function* 495 cm2A,-sec Intrinsic Carrier Concentration ni (Function)* 1.45 x 1010 cm-3 SRH Electron Lifetime tau_nO 0.2 x 10-7 sec SRH Hole Lifetime tau_pO 0.2 x 10-7 sec SRH Electron Concentration Nsrh_n 5 x 1015 cm-3 4 Table 2.1 (Cont.) Nomenclature Program Variable Default Value Units SRH Hole Concentration Nsrh_p 5 x 1015 cm-3 Electron Affinity** - 4.17 eV Donor Energy Level** - 0.044 eV Acceptor Energy Level** - 0.045 eV Saturation Velocity** Function** - cm/sec Electron Auger Coef.** - 2.8 x 10-31 cm6/sec Hole Auger Coef.** - 9.9 x 10-32 cm6/sec * Designates parameters whose values are outputs of a function but are currently set to a constant. This is also in keeping with the theme of future expansion mentioned in the introduction chapter. These terms can be further modeled with little effort. ** Parameters are not currently used. 2.2. Device Equations Semiconductor device phenomenon is described and governed by Poissons equation (2.1). 0V = --(r+P-n) (2.1) where, T = N(x) Nd Na and the continuity equations for electrons and holes (2.2a and 2.2b). Â§l = L.ZL+G-u dt q dx (2.2a) dp__ J_ n dt q dx + G U (2.2b) 5 r is the effective doping concentration defined for the semiconductor, whether it be a PN Junction, Bipolar Junction Transistor, MOSFET, and so on. Equation (2.1) is the differential equation for the potential distribution in an arbitrarily doped semiconductor. This equation, depending upon the particulars, is difficult to solve analytically and in some cases can be impossible. Approximations must be made in most analytical solutions, alternatively numerical methods can be used which is the focus of this report. [1] Because much of semiconductor analysis is concerned with the spatial variation of carriers and potentials from one region of a device to another, Poissons equation is almost always encountered in the solution procedure. It is one of the most useful principles in semiconductor device analysis. Equations (2.2a) and (2.2b) are known as the continuity equations for electrons and holes respectively. These equations are used to describe the free carrier densities and their properties within the device by examining their interactions within a small infinitesimal volume that we can extend across the device as a whole. To derive a one-dimensional continuity equation for electrons, we consider an infinitesimal slice of thickness dx located at x as shown in Figure 2.1. _Â£ ( 1 G R Jn(x) _ Jn(x + dx)^ x + dx Figure 2.1 INFINITESIMAL SLICE FOR DETERMINATION OF INCREASE IN CARRIER DENSITY IN X 6 The number of electrons in the slice can increase because of the net flow into the volume and from net carrier generation in the device. The overall rate of electron increase equals the algebraic sum of (1) The number of electrons flowing into the slice, minus (2) The number flowing out, plus (3) The rate at which electrons are generated, minus (4) The rate at which they recombine. [ 1 ] The first two components (1) and (2) are found by dividing the currents at each side of the slice by the charge on an electron, and for now, the last two components (3) and (4) are simply denoted by G and R, respectively. The rate of change in the number of electron in the slice is then given by equation (2.3). A is the cross-section area of the slice, q is the electron charge, and G and R is the generation and recombination densities respectively. If we use Taylor series expansion on the second term then, dJ Jn(x + dx) = J(x) -I--dx + .... and substituting into equation (2.3) obtains the dx continuity equation for electrons as shown in equation (2.2b). A similar procedure will result in equation (2.2a) for the continuity equation for holes. To yield equations that can be used analytically or numerically we need to introduce the relationship between current density J and the carrier densities n, and p. The current transport equations or current densities are defined by a combination of the drift and diffusion affects on the free carriers within the device. The drift portion of the current expression is simply due to the presence of an electric field, if there is one. This electric field creates an overall drift velocity vd, and is proportional to the electric field and the mobility of the carrier as shown in equation (2.4). dx \ -q '/(*) Jn(x + dx)' s - (2.3) vd =-me (2.4) 7 The mobility of the free carrier is a parameter that describes the ease of motion of an electron in response to an electric field and is affected by several things. The mobility model and how it is used in the application of this program is discussed shortly. The current density flowing in the direction of the electric field can be found by summing the product of the charge on each electron times its velocity over all electrons n per unit volume, or that is The same procedure can be extended to the hole free carriers and then combined to give the total current density due to the applied electric field in equation (2.5). The other affect on the current density through a device is due to the diffusion of free carriers in the device if a special variation of carrier energies or densities exist. To investigate this affect we start by examining a portion of a semiconductor device that is of n-type and with an electron density that only varies in one dimension. The temperature across the region is average so that only the density varies with x and we then consider the number of electrons passing the plane at x=0 per unit time per unit area (Figure 2.2). n ^ n,drift = = fiq^E (2.5) 8 X Figure 2.2 ELECTRON CONCENTRATION DIFFUSION AFFECT Because the electrons are at finite temperature, they only have random thermal motion along one line even with no electric field applied. On average the number of electrons crossing the plane at x=0 from the left side start at approximately x=-X after a collision. The lambda term is defined as the mean-free path for an electron, or that is the average distance an electron will travel before encountering another collision and is given by X = vthTcn. Therefore the average rate per unit area of electrons crossing the plane at x=0 depends upon the density of electrons that started at x=- X and is -^N{-X)vth. The factor of Vi is due to the fact that half of the electrons move further to the left of x=- X after a collision and half to the right. Similarly, the rate at which electrons pass the plane x=0 from the right is ~^N{X)vth which gives a net rate of particle flow from the left is. 9 F = vth [n(-X) n(X)] which can be approximated at x=^ by the Taylor series expansion to give F = v,. dn . dn n(0) X (0) + X dx dx , dn and finally because each electron carries charge q, the particle flow corresponds to a current per unit area in equation (2.6). = -dF = qylv, dn dx (2.6) A more useful expression of equation (2.6) can be developed by using the theorem for equipartition of energy in this one-dimensional case. [1] This allows us to write ~mlv?h ^kT and X vlhTcn, and rearrange equation (2.6) into equation (2.7). hjiff = dDn ^ Wkere n kT Dn=ju q (2.7) The quantity Dn is known as the diffusion constant or Einstein relationship. It relates the two important coefficients that characterize free carrier transport by drift and by diffusion in a solid. Now with the drift and diffusion terms derived, the total electron current density can be expressed by equation (2.8) Jn = Jn.drift + ^n.diff = dMnnE + ^n (2.8) Following the same procedure for holes leads to equation (2.9) and combing the expresses the current density for both carriers as in equation (2.10) 10 (2.9) ^ P J P,drift Jp,diff J = Jn+Jp (2.10) Before we continue to the next section, the equations discussed so far are summarized in equations (2.1 la-e) below which are often referred to as the equations of state [2], and a brief discussion on non-thermal equilibrium conditions is given. The current SDS program does not employ these techniques as solution mechanisms. Equations of State: , dp dip Jp =-qDn qUn P H pdx Hhp dx (2.11a) dn dip Jn = 3" = 1.3-/.+C-(; dt q dx (2.1 Id) ^ =-e(r + -,) (2.1 le) Not all devices operate under thermal equilibrium conditions and so techniques are developed which allow the analysis of semiconductors that are not at thermal equilibrium to occur. We call this quasi-Fermi analysis. Through this analysis the use of quasi-Fermi potential levels are established which relate the current to the gradient of the carrier at that given quasi-Fermi level. The expressions for the quasi-Fermi carrier densities and current densities are shown below in equations (2.12-2.15). 11 n = m( exp (2.12) ( E E ") f = n, exp l kT ) V qWfi-VfnY kT p = n, exp f E E \ c- cfp kT nj exp . dEfn dÂ¥fn Jn= Mnn r~ = -qM" ax ax dE d\j/ fP dx (2.13) (2.14) (2.15) 2.3. Carrier Densities The carrier models used in this program follow the Boltzmann statistics approximation of the Fermi-Dirac distribution model. The carrier densities are approximated in equations (2.16) and (2.17). n~Nc exp f F E ^ r-c kT = "ie eXP ' qty-h)' kT (2.16) p ~ Nv exp (Ev~EfP) = nie exp 'qtyp-V')' { kT j kT K / (2.17) nie is the effective intrinsic carrier concentration and has dependencies on temperature and total doping concentration. This project neglects band gap narrowing due to heavy doping levels but could be implemented as discussed in [3]*. The intrinsic carrier concentration then is only modeled as a function of temperature as shown below in equation (2.18). 12 (2.18) f nie(T) = y/Nc(T)Nv(T)cx p v Eg(T)' 2kT , Nc and Nv are known as the density of states and along with the band-gap energy, E also have dependencies on temperature as shown in equations (2.19-2.21). Eg(T) = Eg(300) + a 3002 300 + p T2 T + p (2.19) NAT) = T 300 AL(300) (2.20) NV(T) = Nv (300) (2.21) Nc (300) and Nv (300) along with a and p, are predefined quantities given in section 2.1. 2.4. Mobility Models The mobility of an electron or hole to move throughout the silicon crystal lattice is representation of the responsiveness of a free carrier to move through the conduction and valence bands in the presence of an electric field. This phenomenon is governed by several complex interactions at the quantum mechanical level. [4] Temperature plays a major role because at higher temperatures the lattice will vibrate more causing the overall free carrier mobility to lower. Dopant levels also have a significant role since the dopant sites create local distortions in the lattice allowing the free carriers to scatter more. 13 The models built into SDS allow for effective doping dependency and electric field dependency as shown in equation (2.22). There is also a factor included to take into account for high-injection levels in high power devices. These models can be turned on or off to account for any dependency that is desired. The constants and fitting parameters are tabulated below. M = /^max Mn 1 + \N(x)\ V "ref J + f 1 1 + f V where, / = \a 2.04 N ref (2.22) Table 2.2 Mobility parameters Nomenclature Program Variable Default Value Units Electron Mobility Concentration Reference* Nrefn 8.5 x 1016 cm-3 Electron Mobility Alpha* alphan 0.72 - Electron Mobility Beta* betan 2 - Electron Critical Electric Field* Ecn 8.0 x 103 V/cm Electron Maximum Mobility* mu_maxn 1330 cm2/V-sec Electron Minimum Mobility* mu_minn 65 cm2/V-sec Hole Mobility Concentration Reference* Nrefp 6.3 x 1016 cm-3 Hole Mobility Alpha* alphap 0.72 - Hole Mobility Beta* betap 1 - Hole Critical Electric Field* Ecp 1.95 x 104 V/cm 14 Table 2.2 (Cont.) Nomenclature Program Variable Default Value Units Hole Maximum Mobility* mu_maxp 495 cm2/V-sec Hole Minimum Mobility* mu_minp 47.7 cm2/V-sec * Constants to calculate field dependant mobility terms, when used. 2.5. Recombination and Generation Models The SDS program models both generation and recombination for use in determining current flow in the specified device when solving both the continuity equations and Poissons equation. The most fundamental type of recombination defined as Shockley-Read Hall recombination. The recombination rates are defined in equation (2.23). u=un=up _______np ~ nf_______ Tp{" + n{)) + Tn(p + p0) (2.23) where, n(t = n(. exp kT and Po = n, exp kT \ The electron and hole lifetime are defined by r and rp, respectively and are defined below in equation (2.24) where o, is the capture cross-section and Nt, the density of trapping centers. T. = and Tn = (2.24) 15 Another method to calculate the lifetime as used in this program is from [3] and is shown in equation (2.25) below. T. = "n 0 N j _|_ total and Tn = > o N N | _j_ total (2.25) SRH.n N SRH.p The constants are tabulated below. Another type of recombination is Auger recombination, characterized by electron-hole pairs disappearing without the aid of trapping centers as in SRH recombination, and is mainly used when there are large high-injection rates. Auger recombination is defined by equation (2.27). This model and the SRH model are added together when used in conjunction to provide the overall recombination rate as shown in equation (2.26). To use the Auger model, the partial derivative terms used in Chapter 3 -Numerical Analysis will need significant modification since Auger recombination is a function of both carrier densities. U=Usm+UM0 (2.26) where, Umjg =r(n2p-p2n) (2.27) The generation model used in this program models the basic generation affect due to impact ionization as defined in equation (2.28). It is useful in characterizing reverse applied bias breakdown. In low electric field applied bias cases the generation rates are typically negligible. [4] The constants and fitting parameters are tabulated below. G=G, =Gr=Ua\j,\ + ar\jl,\) (2.28) where, f / \ / an = o exP - 1 f\ and ap = api) exp \ V 1 i y J V 3 / 16 Table 2.3 Ionization generation parameters Nomenclature Program Variable Default Value Units Electron Alpha* alpha_nO 2.25 x 107 Hole Alpha* alpha_pO 3.08 x 106 Electron Energy* E_nO 3.26 x 106 V/cm Hole Energy* E_pO 1.75 x 106 V/cm Fitting Parameter* m 1 - 17 3. Numerical Analysis In this chapter, we will proceed in setting up the numerical equations in order to solve the physical device equations established in the previous chapter. A procedure for setting up the equations in numerical form, linearizing those equations, and solving the subsequent block tri-diagonal matrix equation will be given in a step by step manner. However, first a brief introduction into some of the basic numerical analysis techniques required to solve the device equations will be given below. All of these techniques are used in transforming the differential device equations in to linearized numerical equations suitable for a numerical analysis solution. 3.1. Taylor Series Expansion The Taylor series is named after the English mathematician Brook Taylor (1685- 1731) and the Maclaurin series is named in honor of the Scottish mathematician Colin Maclaurin (1698-1746) despite the fact that the Maclaurin series is just a special type of Taylor series. But the idea of representing particular functions as sums of power series goes back to Newton, and the general Taylor series was known to Scottish mathematician James Gregory in 1668 and to the Swiss mathematician John Bernoulli in the 1960s. Taylor was apparently unaware of the work of Gregory and Bernoulli when he published his discoveries on series in 1715 in his book Methodus incrementorum directa et inversa. Maclaurin series are named after Colin Maclaurin because he popularized them in his calculus textbook Treatise of Fluxions published in 1742. [5] 18 As we have seen in Chapter 2 Semiconductor Device Physics, the basic equations describing the physics of semiconductors are differential in nature and nonlinear. In order to perform numerical analysis on these equations, they must be linearized and in order to do this, the Taylor series expansion is performed about a point of interest. Additionally, for our equations of interest, the linearization must be done in each variable of interest (n, p, and v|/) and this can be seen later in section 3.5. The general form of the Taylor series is show in equation (3.1) below and is referred to as the Taylor series expansion of the function f at a (or about a, or centered at a). [5] rw = Z /(a), v, -(x-a) =f(a) + (x-a) n=0 n\ f'(a) + 2! ix.af+nsi 3! {x-af +... (3.1) For the special case of a =0, the Taylor series becomes equation (3.2) and is given the special name Maclaurin series. In our analysis, we will be using the standard Taylor series expansion as the expansion point is not always zero. = = f(0) + ^-(x) + ^-(x)2 +^p-(x)3 +... (3.2) n=0 m 1! 2! 3! Since we are dealing with a system of nonlinear equations in which we are interested in multiple variables (n, p, and v|/), each equation that is linearized must be done so around each variable and so therefore the partial derivatives of each variable is introduced. Extending the above discussion on series expansion to a multiple variables is rather straightforward. We will also introduce the concept of perturbation state (or delta state). Let us re-examine equation (3.1) and transform it into the form show in equation (3.3) below. We assume that a is a point such that f(a) = 0 and this point is defined as the equilibrium point of the system f (x) = f(x) since f (x) = 0, when x = a. Since f(a) = 0, the nonlinear differential equation can be approximated near the equilibrium point by equation (3.4). 19 (3.3) jr f(x) = f(a) H----1 x=a (x a) + higher order terms dx f\x) = ^-\x__a{x-a) (3.4) When the expansion point, in our case a, is sufficiently close to x, the higher order terms in equation (3.3) can be neglected as they are very close to zero. To complete the linearization, we define the perturbation state, or delta state (8x = x-a) and using the fact that 8f(x) = f (x) we obtain the linearized model show in equation (3.5). This can be easily applied to multi-variable systems and is show in equation (3.6). This is the form employed in section 3.5 to linearize the system. In summary, this technique allows us to linearize a system of nonlinear equations about a point of expansion. In this analysis, the concept is a more generic approach. As will be seen later, we linearize the current equations and generation and recombination equations as they involve nonlinear terms of an exponential nature. Since the three equations of interest are the two current continuity equations and Poissons equation, and since the current continuity equations are differential expressions of the current equations, these must also be linearized. This allows the linearized current and recombination/generation equations to be substituted in directly. This will be more closely examined in section 3.5. Sf(x) = ^-\*=Sx dx (3.5) (3.6) 20 3.2. Numerical Differentiation It is often necessary, especially in the procedures following, to estimate the derivative of a function and is one of the fundamental requirements in performing finite difference method techniques. The procedure used to perform this is called the difference method and is applicable to first or second order derivatives. The procedure can be used to extend to higher order derivatives but is not discussed here. Additionally, more elegant difference techniques can be used in order to lower the inherent numerical error associated with calculating the derivative of a function but is also not discussed here. Numerical differentiation is performed using a series expansion and in this instance a Taylors series expansion is also used. If we assume that fix) has as many continuous derivatives as necessary, then from Taylors formula, equation (3.7) is found, which is the same form as equation (3.1) but with the expansion about x, instead of a, which was arbitrary for the above discussion. [6] f(x + h) = f(x) + ^f (x) + -^f'(x) + -^f"'(x) + ... (3.7) One way to think of equation (3.7) is that we are expanding the function fix) about x, within small increments of h. So if we replace the arbitrary a in equation (3.1) with x+h, equation (3.7) is a straightforward application. If we solve equation (3.7) for f (x), the equation can be written as h 2! 3! and neglecting higher order terms gives equation (3.8), / (x) = + ^f&l' for small value of h. h (3.8) 21 Equation (3.8) is called the forward divided difference approximation and graphically is the representation of the slop of the line passing through points (x, f(x)) and (x+h, f(x+h)). There are two other main approximations for a first order derivative and they are directly obtained from Taylors series expansion as well. Assume we replace h with -h in equation (3.7) and repeat the same process we obtain equation (3.9) below. Equation (3.9) is referred to as the backward divided difference approximation. / (x) = ^, for small value of h. (3.9) h Finally, if we replace h with -h in equation (3.7) and then subtract the resulting equation from the old one, we obtain /'(*) = f(x + h)~ f(x-h) h hA 2 h - /"'(*)- /<5) (*)-... 3! 5! and neglecting higher order terms gives, / (x) = + ^for small value of h. (3.10) 2 h Equation (3.10) is known as the central difference approximation. Each numerical technique to establish the derivative of a function lends to round-off or truncation errors as we are approximating the value. Keeping higher order terms of the expansion will lower this error but is not necessary for this application. 22 3.3. Numerical Integration Just as with derivatives, calculating the integrals of functions numerically is also of interest for this application. The main use for this technique, recalling the device continuity equations from Chapter 2, is to calculate the current density across the device using the recombination and generation terms that have previously been calculated. Since mesh spacing is also defined, use of a simple trapezoidal integration rule is rather straightforward. Let us suppose we would like to evaluate the definite integral of the form show in equation (3.11) and which is graphically represented in Figure 3.1. (3.11) x0=a X! x2 x3 *n-2 *n-1 Xn-b X Figure 3.1 INTEGRATION USING THE TRAPEZOIDAL RULE 23 To calculate the integral (area under the curve), we begin by splitting up the interval [a, b] into n subintervals, each with uniform mesh spacing as shown in Figure 3.1. This gives equation (3.12). h = and x^a + ih for i = 0,1,2,..., n (3.12) n Now by examining each subinterval, the area under the curve can be approximated by the sum of the areas of the trapezoids. The area of the /th trapezoid is therefore A = | L/C*,-! ) + /(-*, )] with the sum or total area Tn being obtained by extending over the entire area. Tn = A + A2 + ...+A-i + A, = |[/Uo ) + fix, )] + ^[/(*I ) + fix 2 )] + ... + |[/CX.2 ) + fiXn_| )] + ![/(*_, ) + /(*)] Finally, by collecting like terms we arrive at equation (3.13). Tn =^[fix0) + fixn)\+h!Â£f(xi) (3.13) ^ i=i Equation (3.13) is referred to as the composite trapezoid rule as the integral is estimated by a composite of the areas of n trapezoids. [6] 24 3.4. Tri-diagonal Matrices This section will briefly discuss matrices with respect to diagonal dominate matrices and especially tri-diagonal matrices. A tri-diagonal matrix, as shown in Figure 3.2, is simply a matrix which has non-zero entities on the main diagonal, super-diagonal (upper), and sub-diagonal (lower) and frequently arises in the solution of numerical differential equations, such as in this project. a,l a\2 0 0 0 2. a 22 a2i 0 0 0 a)2 a3i a34 0 0 0 43 a4A a45 0 0 0 aS4 a55 Figure 3.2 TRI-DIAGONAL MATRIX The main content of this section will discuss a recursive method for solving a block tri-diagonal matrix equation as this is what is employed in the solution of this project. The algorithm is borrowed from [4] directly but is a standard algorithm used frequently. There are many methods to solve these types of matrices, from the very basic naive Gaussian elimination, to more complex solutions, like LU decomposition, Jacobi method, etc. As we will see in Section 3.5 below, for a 1-D solution, a block tri-diagonal matrix solution must be obtained. The recursive algorithm for this procedure is as follows. Assume the matrix equation that needs to be solved is of the form in equation (3.14). A(N)Sy(N -1) + B(N)Sy(N) + C(N)Sy(N +1) = F(N), 2 where matrices A, B, and C are 3x3 dimension and vector F is of 3x1 dimension. Since A, B, and C are 3x3 matrices we have a block tri-diagonal matrix as opposed to a simple tri-diagonal matrix. A direct solution would be to construct this matrix per equation (3.14) and use any of the aforementioned solution algorithms but this is computationally expensive. It would be much more efficient to only have to calculate the matrix blocks and solve for the delta state (or error value) and not have to manipulate every element in the matrix, including the majority of empty or zero elements. The use of the recursive algorithm allows for this computation savings. To continue, equation (3.14) is transformed into the following equation (3.15) involving unknown vectors at two points, N and N+l, instead of three. B (N)Sy(N) + C (N)Sy(N + l) = F (N) (3.15) Equation (3.15) can be solved for 5y(N) explicitly using basic matrix math, which yields equation (3.16) below. fy(N) = B'(N)F'(N)- B'(N)C' (N)Sy(N + \) (3.16) If we then take equation (3.16) and reduce the division point N by one, we can substitute this result back into equation (3.14), rearranging the terms to arrive at equation (3.17). The following steps illustrate this procedure. Sy(N -1) = B (N -l)F (N -\)-B (N -1)C (N -\)S\'(N) (division point N reduced by one) Substituting back into equation (3.14) yields, A(N)[b\N -1 )F (N -1 )-B (N-1 )C (N l)Jv(A0] +B(N)Sy(N) + C(N)Sy(N + 1) = F(N) Rearranging terms yields, [B(N)-A(N)B (N-\)C (N-l)]Sy(N) + C(N)Sy(N + \) W* ^ / = F(N)-A(N)B(N-\)F(N-l) 26 If we now compare equation (3.17) with equation (3.15) we obtain the recursive relations given below in equation (3.18), but for one less division point calculation. B(N) = B(N)-A(N)B(N-\)C(N-\) C'(N) = C(N) and F\N) = F(N)-A(N)B(N-\)F(N-\), 3 compare equation (3.18) with equation (3.14) for N=2 specifically and obtain equation (3.19). B(2) = B(2), C (2) = C(2), F (2) = F(2) (3.19) We know that Sy(N-l) for N=2 is Sy(l), which is equal to zero. Sy(l) corresponds to the boundary condition error, but the boundary condition is fixed and so there is no error associated with that term and therefore reduces to zero. Finally, starting with equation (3.19) at N=2, then using equation (3.18) for the rest of the division points (N=3, 4,..., L-l) we can use equation (3.16) to calculate the error values, 5y(N), across the whole device (N=L-1, L-2,...,2) starting with 8y(L)=0, since this is also a fixed boundary condition. At this point, matrices A, B, C, and F are arbitrary with the dimensions defined above. As we will see in the next section, these matrices will be defined and so calculation of the above algorithm becomes very straightforward. Once A, B, C, and F are calculated across the device using the recursive relations given, the error potential is solved for directly. The definition of these matrices comes from performing the finite element methods (difference method) on the device equations and then linearizing those equations. 27 3.5. P-N Junction DC Analysis In this section, a procedure for performing DC analysis of a 1-D, P-N junction device will be developed. The basic P-N junction diode is examined as further development techniques for other devices follow the same basic procedures and solutions used in this section. For a 1-D Bipolar Junction Transistor, for example, only slight modification of this procedure is required and the resulting simulator program can be modified to use either solution depending upon the device specified. Another example of modification is time domain analysis, which requires extra division spacing but ultimately requires little modification to the base program with the conceptual idea that these further implementations can be readily incorporated at a later date.* Equations (3.20a-e) and Figure 3.1 show the basic setup for a 1-D P-N junction diode. The overall doping profile as of yet is undefined but could be a step junction, linear junction, and so on. , ^ dp dip (3.20a) dn dip J = aDn QjUnn , OX ox (3.20b) dp 1 n P = p+G Un, dt q dx p p (3.20c) 1 +r l, a,"",'a* " (3.20d) = +p-n), where T = ND-NA dx Â£ (3.20e) For DC analysis the time derivatives in equations (3.20c and 3.20d) will become zero as the device is considered to be in steady state. 28 (N=1) (N=L) Figure 3.3 BASIC P-N JUNCTION DIODE Before we can proceed with numerical analysis of the device show in Figure 2, a few points of interest should be noted. Many of the models used in the solution procedure are of the basic type, however can be expanded easily to incorporate future models. For example, the intrinsic carrier concentration model is currently only temperature dependant and the mobility model is fixed across the device. Both could be modified to incorporate band gap narrowing for heavily doped devices or electric field dependency respectively.* More discussion on each model used in this section along with its capabilities is further discussed in Chapter 5 Semiconductor Device Simulator and also in Chapter 2. The recombination model employed in this solution is of the basic Shockley-Read Hall type. It is the most fundamental type of recombination method and so therefore is most applicable to basic PN junction devices. Equations 18C and d imply that the recombination and generation terms are the net terms and can constitute different models, such as SRH and Auger recombination. In this case, the net recombination rate from both models is used. Furthermore, the recombination centers only occupy a single energy level directly in the middle of the forbidden energy band. This implies that the equilibrium hole concentration is equal to the equilibrium electron concentration which is in turn equal to the intrinsic carrier concentration, no=po=ni. For higher power devices, Auger recombination, as discussed in Chapter 2, can also be considered but is not explicitly used in this solution procedure. This allows us to use the basic SRH recombination equation from Chapter 2, repeated as shown here in equation (21). 29 (3.21) u = uP=v pn nf ^(P + ni) + TP(n + ni) If we only want to evaluate the device in low electric field applications technically the generation terms are not necessary and are not supported in some programs, such as in [3]. However, this program includes the generation terms and as with recombination we assume that the net hole generation terms are equal to the net electron generation terms. These terms are useful in supporting devices with extremely high or reverse applied biases. We can examine the process of avalanche breakdown by incorporating these terms, however when low applied biases are used, these terms are essentially zero. See more discussion in Chapter 5 on the generation function. Repeating the equations from Chapter 2 for the one dimensional case gives us equation (3.22) below. Before we proceed to develop the numerical equations from equations (3.20) necessary to solve the device let us discuss the boundary conditions that are also necessary. Since we are performing a finite element analysis, an initial value or boundary value must be given. We will develop a two-point (Dirichlet type) boundary-value problem. Finite difference techniques can then be employed on the device physics equations and Taylors series expansion to linearize the equations. This leaves a tri-diagonal block matrix that can be solved iteratively and finalized with Newtons iteration principle to arrive at a converged solution. We need to define the appropriate boundary conditions. With respect to Figure 3.3, p and n regions are located at x=0 and x=w, respectively, where w is the device length. For our analysis, we assume that the n-region is fixed at zero potential, or that is the external contact is grounded and the p-region is at any applied voltage. The thermal equilibrium analysis that follows also starts with the p-region fixed at ground. where (3.22) dx 30 The boundary conditions need to satisfy a zero space charge region at each contact. Usually the semiconductor forming an ohmic contact with the metal has sufficiently high doping concentration. If there were a nonzero space charge region in such a highly doped semiconductor then there would be a correspondingly high electric field which would lead to breakdown conditions at the metal-semiconductor contact, so therefore the zero space charge condition must be satisfied. Future work can include the Shottky boundary conditions necessary to satisfy Shottky devices.* HO) + pi0) n{0) = 0, T(w) + p(w) n(w) = 0, p(0)n{0) = nf, p(w)n(w) = nf, ^(0) = V--j-ln 'pi 0)" , wiw) = V~-In niw) 9 >h 9 n, where, 9 = kT (3.23a) (3.23b) (3.23c) Equation (3.23a) satisfies the zero space charge condition specified above, while equation (3.23b) satisfies the thermal equilibrium condition necessary at each contact. Equations (3.23a and 3.23b) are used to establish the two-point boundary conditions that are self-consistent among the two equations and the result is shown in equation (3.24). P(0) = - n(w) = - n 3 1 + 1 + 2n, V 2 J'(0)j T(w) 2 1 + - f 2tli 1 1 + -il/2 ' >, n( 0) = Pi 0) ,-11/21 >, p(w) = n(w) (3.24) Conditions for the potential boundary conditions given in equation (3.23c) are deduced by the use of the quasi-Fermi potentials briefly discussed in Chapter 2. These are re-written here in equation (3.25) for reference. w = 0n------In p 9 P_ ii i >1 3 n 3. u 3. where, 9 = kT (3.25) 31 If the definition of the zero potential point on the n-region and applied bias on the p- region is extended to equation (3.25) it follows that(pp(w) = (pn(w) = 0 and (pp{0) = (pn(0) = V Substituting these into equation (3.25) yields equation (3.23c) directly. With equation (3.23c) and equation (3.24) we now have a two-point boundary condition with respect to the three unknown variable, n, p, and v|/ for use in the solution of the device. We can now proceed with the numerical analysis of the device equations but in order to do that, we need to transform the differential equations into difference equations with the variables defined at a finite number of division points (mesh points) as shown in Figure 3. Main division points are defined by N, with N=1 and N=w corresponding to the device endpoints x=0 and x=w respectively. Auxiliary points are defined such that they are exactly half the way between main division points. This solution uses uniform mesh spacing and is sufficient for most classes of devices, however when the extent of the depletion region becomes much smaller than the neutral bulk regions it is sometimes more beneficial to have course mesh spacing in the neutral regions and finer mesh spacing in the depletion region. The SDS program can be easily modified to incorporate non-uniform mesh spacing*. X=0 1 M-1 -O- N-1 |*- h(M-1)- h(N)- M O- N - h(M) h(N+1)] M+1 ------O- N+1 N+2 I h(M+1)! X=w L Figure 3.4 DIVISION POINTS FOR DC ANALYSIS In the following analysis, variables n, p, and vp are estimated on the main points, while the derivatives of these quantities are estimated on auxiliary points. The device equations in (3.20a-e), which need to be transformed into difference equations, are repeated here in equation (3.26a-e) neglecting the time derivative. [4] dp to-W'P dy/ dx (3.26a) 32 (3.26b) Bn By/ Jn =(lDn^-(iMn , ox ax 1 BJ Q = ~~q'~d^ + Gp~Up' =-^L+c,-t'.. q dx ay_ q dx2 = (F+p-n), where V = ND-NA (3.26c) (3.26d) (3.26e) The equations are now transformed into the difference equations. A brief discussion on equations (3.27a and 3.27b) will follow shortly (including the coefficients for the lambda terms and how they are obtained, but equations (3.27c-e) are obtained directly by using the difference methods employing the techniques discussed in Section 3.2 - Numerical Differentiation. JPW) JAM) q h(M) q h(M) [Apt(M)p(N) + Ap2(M)p(N + \)\ \XnX{M)n{N) + An2{M)n{N + \)\ (3.27a) (3.27b) where, i tits ,A,,V'(N)-ys(N + 1) j ,,,, ,...y/(N)-y/(N + \) 4.1 (M) = Mp(MV ---- AJM) = VAM)Y ---- l-ep i i ,*,MN)-y'W + v Ap2(M) Mp(M) ^nl(M) /J(M) /3(M) = 0E(M)h(M) = 6[y/(N)-y/{N +1)] q Jp{M)-Jp{M-\) h(N) -G(N) + U(N) = 0 (3.27c) q J(M)-Jn(M-1) h (N) + G(N)-U(N) = 0 (3.27d) 33 (3.27e) ytf{N -1) + y2if/{N) + y2i//(N +1) = ~[r(N) + p(N) n(N)] As you can see in the continuity equations (3.26c and 3.26d), the derivative of the current density with respect to the spatial variable, x, has been transformed into the difference approximation in equations (3.27c and 3.27d). Also Poissons equation in (3.26e) was transformed directly to equation (3.27e) also by using the difference method. This time, the second order central difference method was used since Poissons equation is a second order differential equation. The transformation of equation (3.26e) into difference form is shown below in equation (3.28). 1 h (N) y/{N + \)-y/(N) h(M) y/{N)-\fr{N-1) h(M -1) -^[r(A0 + p(A0-(A0] e (3.28) Throughout the solution process, the mesh spacing is continually defined at each division point, whether it is a main division point or auxiliary point. Since we are using uniform mesh spacing and the auxiliary spacing is equal to the main spacing, equation (3.28) could be further transformed into the standard looking central difference equation for a second order derivative. Mesh spacing is left as a function of the division point so that future expansion of non-uniform mesh spacing is rather straightforward*. Lastly, if we examine equation (3.28) with equation (3.27e) the gamma coefficients are obtained directly. K(AO = 1 h(M -\)h (N) y2(N) = h(N) h(M -\) + h(M) y2(N) = 1 h(M)h'(N) At this point in the solution process we have all of the equations transformed into numerical difference equations. The current density equations can be plugged into the current continuity equations and combined with Poissons equation; the three equations can be then solved simultaneously. However, the current density equations, as well as the generation and recombination terms, are not linear and so the application of Newtons iteration principle cannot yet be used. The next step in the solution process is to linearize these equations about the fundamental variables of interest, n, p, and vj/. 34 However, before we proceed to linearize the equations let us return to equations (3.27a and 3.27b), where if we compare these to equations (3.26a and 3.26b) it is obvious they were not directly obtained from the difference technique. Before continuing, please proceed to Appendix A, where a complete derivation and explanation of equations (3.27a and 3.27b) is given. Because the continuity equations contain expressions for the current densities which themselves contain derivatives of the fundamental variables, n, p, and vp, the derivatives for the current densities are defined at main points as well. Therefore, it is necessary to define the spacing between the auxiliary points shown in equation (29). Equation (3.30) shows the hole current density equation linearized through the use of Taylors expansion theorem since it involves nonlinear functions. Since each equation is a function of one or more fundamental variables, the use of partial derivatives is applicable. h'(N)=-[h(M -\) + h(M)] (3.29) Jp(M) = Jp(M) + SJp(M) (3.30) Solving for SJp(M) gives 35 The hole current density is used as an example and the rest of the linearization process for the electron current density, recombination and generation equations is the exact same. Once the continuity equations are linearized, since they contain derivatives of the current densities, the delta state in equation (3.29) can be substituted in along with the rest shown below. dJAM) 8Jn{M)~ Sp(N) + + dp(N) dy/{N) aJ AM) p -Sp(N +1) 8y/(N) + dp(N +1) dl//(N +1) Sys(N +1) dJ(M -1) SJp(M-l) -Â£ Sp{N -1) + ' dp(N -1) dJ(M-\) + dJp(M-1) dy/(N 1) 8y/(N -1) + dp(N) dJAM -1) Sp(N) dy/(N) 8i//(N) dn(N) dn(N + h dn(N +1) KW) di//(N + 1) Sy/(N + \) dn(N 1) ^ K(m-\) dn(N) Sn(N) oy/(N 1) oi/s(N) dG(N) s /hr dG(N) /Xl lx 3G(/V) 8G{N) ~ ^ 1) + ^ ^w(^V-l) + - + dp(N-l) dG(N) Bn(N-l) dy/{N-1) 8y/{N -1) - dG(N) dG(N) _ Sp(N)-\--------Sn(N) H---------5y/(N) + #/(/V) dp(N) y ' dn(N) dG(N) dp(N+ \)~ry' ' a(/V + l) dG(AO af/(A0 diy(N) s ,Kt aG(/v) _ ^ aG(/v) ... ,, 8p{N +1) + ^ Sn(N +1) + Sy/(N +1) dy/(N +1) dp(N) -8p(N) + dn(N) Sn(N) 36 Replacement of Jp(M) by Jp(M) + dJp(M) and Jp(M -1) by Jp (M -1) + dJ (M -1) and so on in equation (25c and d) and rearranging the terms yields equation (3.31). <7 SJp(M)-SJp(M -1) -SG(N) + SU(N) h (N) 1 -1) = .-. P--------G(N) + U(N} q h(N) \_ q SJn(M)-SJJM-l) h'(N) -SG(N) + SU(N) 1 Jn(M)-Jn(M-1) q h{N) -G(N) + U(N) (3.31) We now have the continuity equations (3.31) linearized in terms of incremental variables. The incremental variables (delta state) are given above and these can now be substituted into equation (3.31) which yields the two continuity equations involving only fundamental variable increments as unknowns. Due to the amount of physical space required to show the substitution and rearrangement process, the two equations are given directly in equation (3.32) and (3.33). It is a simple algebraic process but care must be taken that all signs are kept in their proper order and place. Continuity equation for holes... dJp(M- 1)+ dG{N) qh(N) dp(N -1) dp(N -1) 8p(N -1)... 1 ~dJp(M) dJ{p(M -1) dG(N) dU(N) qh(N) dp(N) dp(N) dp(N) dp(N) + a/J(M) dG(N) qh'(N) dp(N +1) dp(N + 1) Sp(N)... Sp{N +1) lG^N)- Sn{N -1)... an(N -1) dG(N) dU\N) dn(N) dn(N) Sn(N)... 37 8y/{N -1)... dG(N) dn(N +1) Sn(N + 1) 1 dJp(M -1) dG(N) qh (N) dy/{N -1) dy/(N-\) 1 ~dJp(M) dG(N) qh'(N) dy/(N) dl//(N) di//(N) Sy/(N)... 1 qh (N) dy/(N +1) 3G(AQ dy/{N +1) dy/{N +1)... = Kw )-Kw- 1)1 + G\N)-u(N) qh (N) L J (3.32) Continuity equation for electrons... 3G()(/V) dp(N -1) dG\N) + dp(N +1) 1 Sp(N- 1) + Sp(N + 1)... dG(N) dU(N) dp(N) dp(N) Sp(N)... dJn(M-\) | 3G(AQ + + qh (N) dn(N -1) dn(N-\) dJn(M) ayn(M-i) dn(N) Â£n(N-l)... I (AO dn(N) dJ(M) ^ aG0(/V) + qh (N) dn(N + 1) 3(Af + l) 3n(N) dn(N) Sn(N + \)... + dJJM -1)+ 3G(A0 qh (N) dy/(N-\) dl//(N-\) + < + 1 1 dy/(N) dl//(N) dJn(M) | 3G(AQ + Sy/(N 1)... 3G0(yv)| a^(^v) J ---------T----- Sl//(N + \)... qh (N) dy/(N + \) 3^(A^ + 1)J -i)]-G0(yv> + G(yv) qh(N)L J Sy/(N)... (3.33) 38 Poissons equation does not involve nonlinear functions and does not technically need to be linearized in incremental variables but for ease of calculation and to match the continuity equations we have done so in equation (3.34). This allows us to setup up a standard block tri-diagonal matrix. The linearization process is not repeated but is the same for the continuity equations. For example, p(N) is replaced by p(N) + dp(N) and so on and then the terms rearranged and grouped. Sp(N)- Sn(N) + y (N)Si//(N -1) + y2 (N)Si//(N) + y3 (N)Sy/(N +1) Â£ Â£ = -^-[r(AO + p(AO-n0(AO]... (3.34) - y (N)Siy{) (N-1) y2 (N)Sy/{) (/V -1) y3(N)Si// (N +1) These three equations (3.32-3.34) are meant to be solved simultaneously and in order to do so we can transform the three equations into a matrix-vector representation that allows us to employ the recursive algorithm discussed earlier. It will also allow us to develop the coefficients of each matrix. Equation (3.35) is the matrix-vector representation of the two continuity equations and Poissons equation combined. A(N)Sy(N -1) + B(N)Sy(N) + C(N)Sy(N + 1) = F(N), for 2 < N < L -1 (3.35) To evaluate the matrix coefficients A, B, C, and F let us examine equation (3.35) in matrix form for three division points, N-l, N, and N+l as shown in Figure 3.5 and then expand each equation in the matrix to form equations (36a-c). 39 8p(N-1) 8n(N-l) 8y/(N -1) ~au(N) al2(N) axi(N) bu(N) bn(N) bxi(N) cxx(N) cx2(N) cX3(N) 8p(N) 7i,W a2i(N) a22(N) a22(N) b2x(N) b22(N) b2i(N) c2X(N) c22(N) c2i(N) 8n(N) = fnW aix(N) ai2(N) a(N) bM(N) bn(N) b}i(N) c3l(V) cn(N) cn(N) 8y/(N) 8p(N + \) 8n(N + \) 8y/(N +1) Figure 3.5 MATRIX FORM OF THE MATRIX-VECTOR EQUATION an (N)5p{N -1) + an (N)Sn(N -1) + al3(N -1 )8yr(N -1) +buSp(N) + bn(N)Sn(N) + bu(N)Si//(N) (3.36a) +cuSp(N + l) + cn8n(N + l) + cl38tp(N + \) = fu(N) a2x (N)Sp(N -1) + a22 (N)Sn(N -1) + a23 (N -1 )Sy/{N -1) +b2XSp(N) + b22(N)Sn(N) + b23(N)Si//(N) (3.36b) +c21 Sp( N +1) + c22Sn{ N +1) + c238y/{ N +1) = /l2 (N) a3x (N)Sp(N -1) + an(N)Sn(N -1) + a33(N -1 )Sy/{N -1) +b3xSp(N) + b32(N)Sn(N) + b(N)Sy/(N) (3.36c) +c3XSp( N +1) + cn8n{ N +1) + ci3Sy/(N +1) = fx3 (N) Now if we compare equations (3.36a-c) with equations (3.32-3.34) respectively, the matrix coefficients for A, B, C, and F are obtained directly. The matrix coefficients are tabulated below for easy reference in Figures 3.6-3.9 40 1 dJp(M -1) dG(N) qh (N) dp(N- -1) WN-\)' 1 dJap(M -1) dG(N) qh (V) dy/(N- -1) dy/{N-1) 1 -1) , dG(N) qh (N) dn(N - -1) dn(N -1) 1 -1) , 5G(A) qh(N) dy/(N- -1) a^(v-i) 0G(AQ Bn(N -1) dG(N) Bp(N-1) 0, ~ 0, ^33 Y\ (AO h(M -\)h(N) Figure 3.6 COEFFICIENT A OF MATRIX-VECTOR EQUATION bu = qh(N) dJp(M) dJp(M-\) dp(N) dp(N) bn bn = dG"{N) dU(N) dn(N) dti(N) 1 \dJp(M) dJp(M -1) dl//(N) dt//(N) qh (N) 3G(AQ dU(N) dp(N) dp(N) dG(N) dy/(N) 9G(AQ dU(N) dp(N) dp(N) qh (AO [ dl//(N) dl//(N) bM=-i b?,2=bn = r2(N) Â£ Â£ dG(N) dn(N) dG(N) dlf/{N) dU(N) dn(N) + -T 1 h(M -1) h(N) Figure 3.7 COEFFICIENT B OF MATRIX-VECTOR EQUATION 41 CI1 CI3 = c22 1 dJp(M) dG(N) dG(N) qh (V) + dp(N +1) C'2 dn(N +1) 1 dG(N) dG(A0 qh (AO dp(N + \) dy/{N +1) Cl' ~ dp(N +1) 1 dJn(M) dG(N) qh '(AO dn(N +1) dn(N -1-1) 1 a/(M) , dG(N) qh '(AO a^(v + o a^(A/ + u c21 o, o, Cj2 y^hi) 1 h(M)h (N) Figure 3.8 COEFFICIENT C OF MATRIX-VECTOR EQUATION / =-{J0p(M)-Jp(M-\)] + Gq(N)-U0(N) qh (N) L J fi2=-^-\jn(M)-Jn(M-\)]-G()(N) + U0(N) qh (N) L J f]J=-^[T(N)+ p0(N)-n(N)]... - y (N)y(N -1) y2{N)y/\N) y3(N)ys(N +1) Figure 3.9 COEFFICIENT F OF MATRIX-VECTOR EQUATION Now that we have the coefficients for the matrix-vector equation determined they can be calculated across the whole device in order to solve equation (3.35). However, as we have seen the recursive algorithm is more efficient and requires less computations then solving equation (3.35) directly. Refer to section 3.4 for the recursive algorithm and appropriate relationships. Before the coefficients for equation (3.35) can be evaluated the partial derivatives embedded in the coefficients must be determined first. The partial derivatives are determined directly from the current density, recombination, and generation equations 42 (3.27), (3.21), and (3.22) respectively. The derivation techniques for several of the partial derivatives are given in Appendix B, you may refer to that now for more information. The required partial derivatives are given in Table 3.1 below for reference. Their calculations are given below in summarized form. No equations numbers are associated with the calculations; they must be referenced to the table numbers below. It may seem excessive but the number of partial derivatives is not as extensive as the list appears. Most of the derivatives can be coded into one function as we will see in Chapter 5. The first table below simply lists all the partial derivatives required to support the matrix coefficients A, B, C, and F. Table 3.1 Partial derivatives required for matrix coefficients 1) (lfl) dJp(M-1) (1 b) dJp(M) (1 c) cUp(M) (\d) dp(N-1) MN) dp(N) dp(N +1) dJp(M-1) (2a) cUp(M- 1) (2b) dJp(M) (2c) dJp(M) (2d) dyr(N-l) dlf/(N) dy/(N) dl//(N + 1) 1) Bn(N-1) (3 a) dJn(M-1) dn( N) (3 b) Mn(M) dn(N) (3c) dJ(M) dn(N + 1) (3 d) K(M-X) (4 a) 37(M-1) (4b) cUn(M) (4c) dJ(M) (4 d) dy/(N -1) dl//(N) dys(N) di//(N + 1) dU(N) dn(N) (5b) dG(N) dp(N- 1) (6a) dG(N) MK) (6b) dG(N) dp(N + 1) (6c) dU(N) M N) (5a) dG(N) dn(N- 1) (la) dG(N) dn(N) (lb) dG(N) dn(N + 1) (7c) dG(N) dy/(N -1) m dG(N) dl//(N) m dG(N) dt//(N + 1) (8c) Referring to the current density equation from Appendix A at division point M-1 and M repeated below in equation (3.37), the partial derivatives for (la-d) are obtained directly. 43 (3.37) M" - = " 'W* +-l)p(Af)] 1) 0/7(A^ 1) dJp(M-1) a^(^v) topm = dp(N) dJp (M) _ 3p(W + l) ~ 7 77^pi (M I) /?(M -1) MM -1) Xp2(M-l), ^i ,(M), h(M) p' h(M) Xp2(M), (1) (lc) (W) The partial derivative of the current density with respect to the potential, vp is more complicated and requires use of the product rule and chain rule for derivatives since the lambda terms are functions of the potential themselves. The results are tabulated below. The lambda terms calculated as a part of the current density equations are repeated here for reference in equation (3.38). The partial derivatives for (2a-b) are tabulated below. /I, (M) = jUp(M) + \) y/{N)-y/{N +1) (M-\) = jUp(M-\) (M-\) = jup(M-\) 1 -em) W(N \)-y/(N) y/(N -\)-if/(N) l e 0(M) fi(M) = 9\y/{N)-y/{N +1)] P(M -1) = e[y/{N -1) ^(AO] where (3.38) 44 dyKN-1) h(M-\){ dyKN-1) dyKN-1) (2a) where, d\\(M ~ 1) //P(M-1) tep2(M-\)= Vp(M-1) a^(yv-i) _(i-Vw-'>) P(M- 0(M- (l-e*"-") dJp(M-l)_ dJp(M-l) dyKN) dyKN-\) where, 9^,(M-1)= 8^,(M-1) a^(yv) ~ a^(yv-i) a^cM-i)^ a^2(M-i) a^(yv) a^(yv-i) a/;(M) dyKN) KM) P(N) dApi(M) dyKN) +p(N+l) ayM) a^Ao (2c) where, Mplm_ MP(M) j P{M)e~PiM) a^(^v) _ (\-e'P(M)) (l-e~PiM)) a^2(M)_ JUP(M) , /KM*** dy/(N) 45 afp(M) afp(M) dyKN+l) ~ di/KN) where, dApl(M) a/t,(M) a^(yV + l)_ dy/{N) dXpl{M) = dJp(M) dy/(N + l) dy/(N) The procedure for evaluating the electron related partials derivatives in 3a-d and 4a-d are exactly the same. The partial derivatives for these terms are tabulated below. 1) dn(N-l) dJn(M-1) dn(N) = dn(N) h(M -1) (M-l), q h(M -1) q h{M) ayn(M) dn(N +1) q h(M) A2(M), (3a) (3 b) (3c) (3a) dyKN-1) h{M-1) l a^yv-l) di/XN-1) J where, dAn](M-\) = jun(M-\) dXp2(M-1) a^(yv-i) mp(m-1) a^-i) aAn2(M-i)^//(M-i) a/ip,(M-i) a^(A^-l) jup(M-\)' dys(N-l) dJn(M-1) dy(N) q h(M -1) n(N -1) ain,(M-1) dy/(N) + n(N) dXn2(M-\) m 46 where, dAn[{M-\) = dAni(M -1) jUn(M -1) dAp2{M-1) dy/(N) dy/{N-1) jUp(M-\) di//(N-\) dAn2(M-\) = dAn2(M-1)^ jUn(M-1) a^CM-l) diy(N) dy/(N-\) Mn(M -1) a^(/V-l) tym= q dyKN) h(M) (A0 a^AO +n(A+l) a^2(M) a^Ao (4c) where, dAjM) ^yn(M)dAp2(M) diy(N) jUp(M) dy/{N) dAn2(M) ^jUn(M)dApl(M) di//(N) jUp(M) a^(AO ay(M) a^(/v + i) 4 h(M) n(N) a/in,(M) a^ + D + n(A + l) dAn2(M) dy/{N +1) (4J) where, a/l,(M) ^ BAni(M)= yn(M)dApl{M) dy/{N + \) dys(N) UP{M) dy/{N) dAn2(M) dA2(M)_ jun(M) dApi(M) dy/{N +1) dy/(N) jUp(M) di//(N) The amount of partial derivatives for the current density equations may seem excessive, but as will be explained in Chapter 5, only one function for the lambda coefficients and partial derivative lambda expressions are required to calculate all of the above which allows the calculation for each of the matrix coefficients to be more efficient. The partial derivatives for the recombination and generation terms are also given below. Their derivations are similar and briefly discussed in Appendix B as well. 47 Recombination terms: dU(N) _ n(N)-fnU(N) dp(N) Tn(p(N) + ni) + Tp(n(N) + ni) (5a) p(N) + rpU(N) dU(N) ____ dn(N) Tn(p(N) + ni) + T (n(N) + ni) (5b) Generation terms: x,M.X).a,N) (6a) dp(N-\) h(M -1) 7,1 p dG(N) MH) = + m,. mu ar(N) (6b) dG(N) m. dp(N + \) h(M) = ^77-Ap2(M)-an(N) (6c) where, h(M) h(M -1) mn = ;---, mh =----;---- fl 2 h(N) 2h(N) ap(N) = aPo exP and E(N) = maE(M -1) + mhE(M) The choice of the sign for the partial derivative terms depends upon the sign of the current density at the given mesh point, that is J (N) >0 or Jp (N) < 0 determines + or respectively when, Jp(N) maJp(M -1) + mbJp(M) 48 m An[{M-l)-an{N) {la) ZG{N) _ . dn{N-\) h(M -1) dG(N) dn(N) m h(M -1) K2{m-\) + mu h{M) a{N) {lb) dG{N) dn{N +1) m. -+- b h{M) Xnl{M)-an{N) (7c) h{M) h{M -1) ------- ffl ~ --------- 2ti{N) h 2h{N) an(N) = anO eXP En0 |Â£(/V)| and E{N) = maE{M -l) + mbE{M) The choice of the sign for the partial derivative terms again depends upon the sign of the current density at the given mesh point, that is Jn {N) >0 or Jn {N) < 0 determines + or respectively when, Jn{N) = maJn{M -\) + mhJn{M) The choice of the sign or + for the remaining partial derivatives below depends upon whether Jn{N),Jp{N),E{N)>0 or Jn{N),J p{N), E{N) < 0 respectively. For example, G^n is positive is Â£(A0>0 but negative if Â£(iV) < 0, however G 3I is negative is E{N) > 0 but positive if E{N) < 0. The term that governs the sign follows each GÂ¥ij below. 49 dG(N) dy/(N -1) ~ GÂ¥\\ + GÂ¥\2 + GÂ¥U + GÂ¥\A (8a) <*G(N) - 3G(7V) r r r r 3^ + 1)Gi+G+C'm+G'')4 (8c) where, _ _L fy ( AT\ P Gy/\\ ~ ~ n aP\N) 2 h(M-1) E(Nf q ma(N) dJ (M -1) G,=-2-2----^-------, Jp{N) E(N) >12 q a^(N-l) m. vn 0 G'":' +A(M-l)a" JmIN) E(N) Â¥14 q dy/{N-\) ' Gr3i=+T^7T a,W)-Â£j h(M) E(N) Gvn ~ + _mhap(N) dJp(M) q dy/(N) G^=*Trt:a.m^ h(M) K(AO E(NY q JAN) J{ff 34 dlf/( N) E(N) JAN) E(N) 50 ^21 ~ ^Vll ^V_31 GÂ¥22=ap{N) m, 3./_Â£m dl//(N) dJp(M -1) q di//(N-1) ^V23 ^Vl3 (//33 GvU=an(N) q Mn(M) dl//(N) ma dJn(M -1) q 0^(A-1) J M*) With all of the coefficients for the matrix-vector form of the continuity equations and Poissons equation tabulated and the partial derivatives for those coefficients calculated, we can use them to develop an algorithm which calculates the thermal equilibrium or applied DC bias potential and carrier densities throughout the device. 3.6. Numerical Solution This section will discuss the numerical solution process used to solve the matrix- vector equation (3.35), including development of the trial values, or initial guess, and finally Newtons iteration principle to update the trial values and check for convergence of the solution. A basic solution process is also discussed which was used to develop the program given in Appendix C. This solution algorithm involves solving the two continuity equations and Poissons equation as previously discussed. The process of transforming these equations into finite difference equations allows the use of the implicit method as a solution procedure which is simply a multidimensional Newton iteration problem. Equations (3.32-3.34) have been developed such that full coupling exists between the three fundamental variables, n, p, and v|/ or that is each fundamental variable is allowed to change during each iteration. In this manner the equations can be solved simultaneously in a straight forward fashion. Other methods exist to solve the numerical equation developed including decoupled methods (Gummel method [3]) but are not discussed here. They have their various differences and benefits for specialized situations but for a generic solution algorithm for most possible device conditions (low injection vs. high injection for example) the couple Newton iteration 51 procedure is preferred. It will increase the overall number of computations per iteration but decrease the overall computation time for convergence. The solution to equation (3.35), repeated below for reference and show in matrix form, can involve a rather significant amount of matrix elements. An efficient solution procedure is required or the computation time, even for simple devices, can become excessive and limits the use of such simulator tools to any practical sense. Use of the coupled Newton iteration principle with the recursive solution procedure will suit the needs of this program adequately. A(N)Sy(N-\) + B(N)Sy(N) + C(N)Sy(N + l) = F(N), for 2 ~p(N)~ Sp(N) y(N) = n(N) , Sy(N) = Sn(N) y(N) Sy/{N) B( 2) C(2) 0 0 Sy( 2) F{2) A( 3) B( 3) C(3) 0 Sy( 3) F( 3) 0 A(4) B( 4) c(4) : 0 A{L 2) B(L- 2) C(L 2) Sy(4) F(4) 0 0 Sy(L- 2) F(L-2) 0 0 A(L-l) B(L-l) Sy(L-l) F(L-l) Newtons iteration principle requires that we have an initial guess or starting value to calculate the error state for the first iteration. The calculated error is then added to the initial value and the process repeated until convergence is reached. All devices, whether simple 1 -D P-N junctions or elaborate 2-D bipolar devices, should have the thermal equilibrium condition solved for first, or that is no external applied bias. Not only does it allow for faster computation speeds on the applied bias conditions but allows the user to obtain a good feel for whether the devices equilibrium conditions are being met as intended. 52 For the thermal equilibrium condition, trial potentials are simply provided based on the charge neutral assumption that the boundary conditions were established under. n2 Carrier densities are given by pT = -T, nT = p, where T is the effective doping concentration throughout the device. The zero-bias trial values for the potential are given by equation (3.39) below. -In e ____ r p region n region (3.39) Now that the trial potentials are defined we can conceptualize a solution procedure, which is given below. A flow chart representation is provided in [4] as well as the authors version in Chapter 5, which represents this programs solution. (1) Define physical and numerical parameters (2) Define iteration limit and set current iteration number to zero (3) Use the matrix-vector coefficients, A, B, C, and F to calculate the error for each iteration using the recursive algorithm discussed in Section 3.4. This is equivalent to solving equation (33) using direct Gaussian elimination for example. (4) Improve the initial guess or trial potentials by adding the calculated error. (5) Check for convergence against defined convergence tolerance (6) End program if solution has converged otherwise update current iteration limit by 1 and return to step (3) repeating the process until solution has converged or iteration limit is reached. Assuming that a converged solution for the equilibrium condition has been reached, extending the solution to an external biased condition only requires updating the equilibrium condition potential with a simple proportional relationship and then repeating the process defined above. Technically updating the potential is not necessary but will lower the total iteration number for the DC condition. This relationship is given in equation (3.40) and illustrated in Figure 3.10 below. 53 i//t(N) = (3.40) v2-v, y/,(L)-y/x{L) V -V ^(N)+- 2 1 y/,{L)-y/^{L) V\{L) Figure 3.10- UPDATED POTENTIAL BASED ON PREVIOUS POTENTIAL [4] A tabulated list of physical constants and user defined constants as well as a tabulated list of numerical constants required to develop the program in Appendix C is given below. Table 3.2 Numerical constants Nomenclature Program Variable Default Value Units Device Length length 10 um Convergence Tolerance etol 1 x 10-6 - # Grid Points L 1001 - Iteration Number it 0 - Iteration Limit itl 15 - 54 Table 3.3 Physical constants Nomenclature Program Variable Default Value Units Dielectric Permittivity Constant eo 8.854 x 10-14 F/cm Electronic Charge q 1.602 x 10-19 C (A-sec) Boltzmanns Constant k 8.62 x 10-5 eV/K Thermal Voltage Vt 0.02586 V Boltzmann Factor (Thermal Voltage) l/Vt 38.6 1/V Temperature T 300 K Donor Concentration Nd 1 x 1015 cm-3 Acceptor Concentration Na 1 x 1015 cm-3 Left Contact VI 0 V Right Contact Vr 0 V 55 4. SDS Results Investigated In this chapter, we will review the output of the SDS program and compare it to analytical results as well as to other simulators, notable SimWindows. [7] SimWindows is also a 1-D semiconductor device simulator capable of simulating the same devices as the SDS program. The first section will compare the SDS results with respect to standard analytical models, for example the built-in potential with respect to using the depletion approximation as opposed to solving the numerical solution which makes no such assumptions. The purpose of this is to verify the working operation of the program versus general semiconductor operation that is learned at the classroom level. 4.1. SDS Output Compared to Analytical Assumptions Most basic PN Junction analysis is performed using the analytical device expressions as defined in Chapter 2. To make these differential equations mathematically tractable, assumptions must be made in order to arrive at a solution. One of the main assumptions used in this approximation is the Depletion Approximation. [ 1 ] The Depletion Approximation is a way to simplify Poissons equation such that the potential throughout the device can be calculated without the need of the carrier concentrations, n and p. The approximation is performed by assuming that the transition region between the n and p bulk regions, or that is depletion region, is totally depleted of mobile carriers and that the mobile majority carrier densities abruptly become equal to corresponding dopant concentrations at the edges of the depletion region. The charge density is therefore zero everywhere except in the depletion region, where it equals the ionized dopant concentration. Additionally, we generally assume that the dopant atoms are complete (100%) ionized and free to contribute to current conduction. This reduces Poissons equation to what is shown in equation (4.1) below. 56 (4.1) d2(p dx2 If we further only evaluate Poissons equation with respect to the n-region of the PN junction, equation (4.1) reduces to equation (4.2) below. This reflects no acceptor atoms in the n-region and can easily be integrated from an arbitrary point in the n- region depletion region to the edge of the region resulting in equation (4.3). Similar results are obtained in the p-region in equation (4.4). d2(/> dE (4.2) qN, E(x) =----(*-*) 0 QNa E(x) =-------(jc + jc ) x (4.4) If we further integrate equations (4.3) and (4.4) one more time, we arrive at expressions for the potential in each respective n and p-region shown in equations (4.5) and (4.6). n
-xp |

Full Text |

PAGE 1 A PN SEMICONDUCTOR DEVICE SIMULATOR USING MA TLAB by Shawn Randall Pace B . S.E.E., University of North Florida , 2007 A thesis submitted to the Univer ity of Colorado Denver in partial fulfillment for the degree of Master of Science Electrical Engineering 2011 f PAGE 2 This thesis for the Master of Science Electrical Engineering degree by Shawn Randall Pace has been approved by Hamid Fardi I 1 I 19-( I I Date PAGE 3 Pace, Shawn Randall (M.S., Electrical Engineering) A PN Semiconductor Device Simulator Using MATLAB Thesis directed by Professor Hamid Fardi ABSTRACT The purpose of this project is to develop a functional PN semiconductor device simulator that is modular in nature in order to allow for flexibility during programming and to allow for future modification with relative ease. The device modeling program is developed using MATLAB . In addition, the program's main goal is to provide a simulator tool that can supplement device modeling, semiconductor device physics, and numerical analysis to construct basic semiconductor devices which can be studied using standard numerical analysis techniques. MA TLAB ' s capability and inherent nature of handling matrices and matrix operations makes this approach an excellent technjque to develop numerical analysis algorithms. The program solution will be used to examine device parameters such a s carrier statistics , device potential , and internal electric fields. The device olution is compared to analytical approximations in order to further strengthen the understanding between theory and exact numerical solutions and how those solutions are obtained. This abstract accurately represents the content of the candidate's thesis. I recommend its publication . Signed Hamjd Fardi PAGE 4 ACKNOWLEDGEMENT The main credit and acknowledgement for the development of this project goes to Dr. Hamid Fardi, without whose knowledge and support would have made the creation of this project nearly impossible. His knowledge of numerical analysis related to semiconductor physics and device physics operation as a whole was immensely appreciated and invaluable . Furthermore, appreciation should be given to the University of Colorado Denver , whose Electrical Engineering department and faculty provided the necessary training and knowledge over the course of my work there to make thjs project possible. It is my sincere hope that this tool can be used in a valuable training manner and supplemental learning aid. PAGE 5 TABLE OF CONTENTS Figures . . .. . .. . .. . .. . .................. . .. ...... .. . ................... . .. .................. .. . .. . . vii Tables ................... . . . . .. . .. . .. ................ .. . .. . .. .. . . .. . .. ............................ .ix Chapter I . Introduction . . .. . .. .. . . . ... . ...................... .. . .......... ................ . . .. .. ...... 1 2 . Semiconductor Device Physics ....... .................. . .. .. . ....... .. . .............. .4 2 . 1 Semiconductor Material Parameters . ....................................... .. . ...... .4 2.2 Device Equation s ...... . .. ................. . ........ . .. . .. .. . .. . .. . .. . .. . .. .............. 5 2.3 Carrier Densities ..... .. . .. .. . .. .. .. . . . . .. . .. ................................. ....... .. .12 2.4 Mobility Models . . .. . .. ....................... .. . .. . .. .. . . ... . . . . .. ....... . .. . .. ....... 13 2 . 5 Recombination and Generation Model . ......... .. . . ..... . ........................ 15 3. Numerical Analysi s . .. . . ............... .. . .. . .. ........................ .. . . .. . . .. .. . . .18 3.1 Taylor Series Expansion . . .. ...... . .. ..... ... . . .. . ....... . ....... . .. . .. ............... 18 3.2 Numerical Differentiation .................... .. . . . . . ......... .. . .. .. ........ . ...... .. 21 3.3 Numerical Integration ................................................................. 23 3.4 Tri-diagonal Matrices . ............................................................. .... 25 3.5 PN Junction DC Analysis ................ . .. . .. . .. . .. ........... . . .. . .. . .. . .. . .. ..... 28 3.6 Numerical Solution of PN Junction ................................................. 51 4. SDS Re s ult s Investigated . . ................. . . .. .. . . .. ...... ........... . .. .. . . . .. .... 56 4.1 SDS Output Compared to Analytical Approximations ....... .. . ...... ........... 56 v PAGE 6 5 . Semiconductor Device Simulator . . .. . .................................... .. ...... . . 67 6 . Summary . . .. . .. . .. . ..... . ............. .. .. .. . . . . .......... . . .. ...................... . ... 73 Appendix A. Integral Form of the Current Den s ity Equation ........ .. . . . ... .................... 74 B. Derivation Technique s of the Matrix-vector Coefficient s . . .. . . . ................. 81 C. Semi c ondu c tor Device Simulator Implemented ........... .. . . .. .. . .. . ............ 85 Reference s .... . . ........ . . . ........................... .. .................................... ... . 114 vi PAGE 7 LIST OF FIGURES Figure 2 . 1 INFINITESIMAL SLICE FOR DETERMINATION OF INCREASE IN CARRIER DENSIT){ IN )( ..................... .. . .. . .. ........... 6 2 . 2 ELECTRON CONCENTRATION DIFFUSION AFFECT ....... .. . . . .. ...... . . 9 3.1 INTEGRATION USING THE TRAPEZOIDAL RULE .. ...................... 2 3 3.2 TRI-DIAGONAL MATRJ)( .......... . . . . . ....... . . . . . . . ............. . ...... . ...... . ... .............. . 25 3.3 BASIC P-N JUNCTION DIODE ......... .. .. . . .. .. . . ...... . . . ........ .. . . .. ...... . 29 3.4 DIVISION POINTS FOR DC ANAL ){SIS . .. ............ ........ .. . . ........ ..... 32 3 . 5 MATRJ)( FORM OF THE MATRJ)(-VECTOR EQUATION ................ . .40 3 . 6 COEFFICIENT A OF MATRJ)( VECTOR EQUATION . . .. ........ . . .. . .. .. . .41 3.7 COEFFICIENT B OF MATRJ)( VECTOR EQUATION ... ........ . ........ . . .41 3.8 COEFFICIENT C OF MATRJ)(-VECTOR EQUATION ........... . .......... .42 3 . 9 COEFFICIENT F OF MATRJ)(-VECTOR EQUATION ....... . . ... .......... .4 2 3 . 10 UPDATED BIAS POTENTIAL BASED ON PREVIOUS POTENTIAL . ...... . . . ..... ................... ... . . . . . . . ......... . . . . ............. . ...... . . ........ . . . . . . . .... 54 4 . 1 ABRUPT STEP JUNCTION PROPERTIES ..... . . . .. ......... .................. 59 4 . 2 EQUILIBRIUM DEVICE POTENTIAL. ............... . . .. .. . ....... .. .. ....... . 60 4 . 3 EQUILIBRIUM ELECTRIC FIELD ................. .. . .. . ............ ......... . . 61 4.4 DEVICE POTENTIAL WITH UNEQUAL DOPING DENSITIES . ......... 62 vii PAGE 8 4.5 SIMWINDOWS COMPARISON TO FIGURE 4.4 ............................. 62 4.6 SDS OUTPUT SHIFTED TO ZER0 .............................................. 63 4.7 SIMWINDOWS ELECTRIC FIELD .............................................. 63 4.8 SDS ELECTRIC FIELD ....... ...................................................... 64 4.9 SDS DEVICE POTENTIAL WITH EXTERNAL APPLIED BIAS . . .. . .......... . . .. ... ........................ ..................... .... 65 4.10 SIMWINDOWS POTENTIAL WITH EXTERNAL APPLIED BIAS .......... . .. .. . . .. ......... . .. .. . . .. .. . ............................ ... 66 5.1 SDS BASIC FLOWCHART. ...................................................... . . 72 A.l INTEGRAL FORM OF CURRENT DENSITY EQUATIONS ................ 79 viii PAGE 9 LIST OFT ABLES Table 2.1 Material parameters for silicon .. .. .. . . .. . . ...................... . . .. ........... .. . ... .4 2.2 Mobility parameters . . .. . ....................... .. . ....... ............................. 14 2.3 Ionization generation parameters .......................... .......................... 17 3.1 Partial derivatives required for matrix coefficients .............................. .43 3.2 Numerical constants ............... .. . .. .. ..... . ...................................... 54 3.3 Physical constants ....................................................... .. ............ 55 ix PAGE 10 1. Introduction Although there are several powerful industry and student level device physics simulators already available, the hope is that with this new tool and supplemental report , the theory of numerical analysis applied to device physics can be studied at a more fundamental level. Traditional simulators allow for device simulation but do not offer any insight into how that solution was obtained . The addition of the report that supplements the Semiconductor Device Simulator will develop the basic tools necessary to understand the operation of the program and allow future modifications as necessary. This project was developed in the programming environment provided by MA TLAB which is developed and supported by Math works . Almost all academic institutions use MATLAB heavily in the learning process and it is readily available to most students, hence the choice of this programming environment. Additionally , MA TLAB is an extremely powerful numerical analysis tool that makes the solution of a program like this one rather straightforward. MATLAB 's capability and inherent nature of handling matrices and matrix operations makes this an excellent tool to develop numerical analysis algorithms. The current semiconductor simulation tools available do not lend to updates or modifications in a simple or straightforward manner. In fact most tools are not open to be changed in any manner unless the source code is available , which is rare. With the open nature of MA TLAB the solution process can be st udied throughout and all variables and parameters are available for examination. The program could have the capability to be standalone which allows for the use of the tool without the requirement of the MA TLAB environment.* The symbo l " *" used throughout this manual is to indicate a future development not yet built into the current version of the program. Almost all sections in this manual include a future development section to indicate where the code can be easily updated to include new functions or models as necessary. 1 PAGE 11 The remaining discussion in this section gives a brief overview of the contents included in the following sections. Chapter 2 gives a brief introduction to the devices physics required to build a numerical analysis program to simulate semiconductor devices in general. Device physics equations (continuity, Poisson's, current, etc.) are discussed as are the models, such as mobility, SRH (Shockley-Read-Hall) recombination , and generation. Boundary conditions are also discussed. Lastly , ilicon material parameters are given and briefly discussed. Currently no other semiconductor material is used in this program , but the addition of other materials is rather straightforward* . Chapter 3 is broken into several sections. The first section will give a brief introduction to numerical analysis techniques such as differentiation using the difference method, integration using the trapezoidal technique, Taylor series expansion, and basic tri-diagonal matrix so lutions. The next section will discuss the numerical parameters required to support the program, such as mesh spacing, iteration limits , convergence tolerances , etc. The final section will give an in-depth view of setting up and solving the numerical equations for a basic PN junction device. The section will cover converting the basic device physics equations in Chapter 2 into numerical equations, setting up the tri-diagonal block matrix, performing the recursive solution algorithm and determining convergence criteria. Chapter 4 will provide simulation results and compare them to basic theory (depletion width, built-in potential, etc.) as well as with other Simulation tools, such as Sim Windows. Chapter 5 will discuss the program as written and the various functions built to support it. Each function will be discussed and explained. It should be noted however , that the MATLAB code itself is he avily commented as a supplement to this section . Chapter 6 gives a brief summary of the report. 2 PAGE 12 Appendix A Integral Form of the Current Density Equation Appendix B-Derivation Techniques of the Matrix-vector Coefficients Appendix C Semiconductor Device Simulator Implemented A final note on the Introduction is in order. Currently this program supports basic 1-D PN junction devices with uniform mesh spacing. Although the program is hardcoded for a Step Junction, the junction profile is easily changed just as the numerical parameters (convergence tolerance, relaxation factor , device length and doping) are without any further modification to the program. Modifications for 2-D solutions, refined mesh spacing, Shottky diodes, and BJT' s could be further implemented* . The same can be said for different numerical solution techniques if desired. 3 PAGE 13 2. Semiconductor Device Physics In thjs chapter, a brief discussion of simple device physics will be covered along with the differential equations governing them as well as any necessary models required to establish a numerical solution of the differential equations. Semiconductor physical parameters will be given in Section 2.1, device equations discussed in Section 2 . 2, carrier densitie s in Section 2.3, mobility model s in Section 2.4 and recombination generation model s in Section 2.5 . 2.1. Semiconductor Material Parameters A tabulated list of the semiconductor physical parameters necessary to reach a numerical solution and some extra parameters used for future expansions are given below. Currently only parameters for Silicon material are used. The corresponding program variable is also shown. Table 2.1 -Material parameters for silicon Nomenclature Program Variable Default Value Units Relative Permittivity er 11.7 Silicon Permittivity es 1.036 X 10-12 F/cm Energy Gap (300K) EG300 1.08 eV Alpha EG_alpha 4.73 X 10-4 -Beta EG_beta 636 -Conduction Band NC300 2.8 X 1019 cm-3 Density (300K) Valence Band NV300 1.04 X 1019 cm-3 Den sity (300K) Electron Mobility Function* 1350 cm2/V-sec Hole Mobility Function* 495 cm2/V-sec Intrinsic Carrier Concentration ni (F unction) * I .45 x 1010 cm-3 SRH Electron Life time tau_nO 0.2 X 10-7 sec SRH Hole Lifetime tau_pO 0.2x 10-7 sec SRH Electron Concentr atio n N sr h _ n 5 X 1015 cm-3 4 PAGE 14 Table 2.1 (Con't.) Nomenclature Program Variable Default Value Units SRH Hole Concentration Nsrh_p 5 X 1015 cm-3 Electron Affinity** -4.17 eV Donor Energy Lev e l** 0 . 044 eV Acceptor Energy Level** 0 . 045 eV Saturation Velocity** Function** -em/sec Electron Auger Coef . ** -2 . 8 X 10-31 cm6/sec Hole Auger Coef. ** -9 .9x 10-32 cm6/sec * Designates parameters whose values are outputs of a function but are currently set to a constant. This is also in keeping with the theme of future expansion mentioned in the introduction chapter. These terms can be further modeled with little effort. **Parameters are not currently used. 2.2. Device Equations Semjconductor device phenomenon is described and governed by Poisson ' s equation (2.1). iiljf q -=--(r+ p-n) dx2 (2.1) where , r = N(x) = N d N a and the continuity equations for electrons and holes (2.2a and 2.2b) . (2.2a) dp 1 dl p -=---+G-U dt q dx (2.2b) 5 PAGE 15 r is the effective doping concentration defined for the semiconductor, whether it be a PN Junction, Bipolar Junction Transistor, MOSFET, and so on . Equation (2 .1) is the differential equation for the potential distribution in an arbitrarily doped s emiconductor. Thi s equation, depending upon the particulars, is difficult to solve analytically and in s ome cases can be impossible. Approximations must be made in most analytical solutions, alternatively numerical methods can be used which is the focus of this report . [1] much of semiconductor analysis is concerned with the s patial variation of carrier s and potentials from one region of a device to another , Poisson's' equation i almost alway s encountered in the solution procedure . It is one of the most useful principles in semiconductor device analysis . Equations (2.2a) and (2.2b) are known as the continuity equations for electrons and holes respectively. These equations are used to describe the free carrier densities and their properties within the device by examining their interactions within a small infinitesimal volume that we can extend across the device as a whole. To derive a one-dimensional continuity equation for electrons, we consider an infinitesimal slice of thickness dx located at x as shown in Figure 2.1. E " J , ( x ) G R J ,(x+dx) 0 0 E., X x+dx Figure 2.1 -INFINITESIMAL SLICE FOR DETERMINATION OF INCREASE IN CARRIER DENSITY IN X 6 PAGE 16 The number of electrons in the slice can increase because of the net flow into the volume and from net carrier generation in the device. The overall rate of electron increase equals the algebraic sum of (1) The number of electrons flowing into the lice, minus (2) The number flowing out, plus (3) The rate at which electrons are generated, minus ( 4) The rate at which they recombine. [ 1] The first two components (1) and (2) are found by dividing the currents at each side of the slice by the charge on an electron, and for now, the last two components (3) and (4) are simply denoted by G and R, respectively. The rate of change in the number of electron in the slice is then given by equation (2.3). an Adx=(J"(x)-J"(x+dx)JA+(G-R)Adx ax -q -q (2.3) A is the cross-section area of the slice, q is the electron charge, and G and R is the generation and recombination densities respectively. If we use Taylor series expansion on the second term then, J "(x+dx) = J ,(x)+ aJ, dx+ .... ax and substituting into equation (2.3) obtains the continuity equation for electrons as shown in equation (2.2b) . A similar procedure will result in equation (2 .2a) for the continuity equation for holes. To yield equations that can be used analytically or numerically we need to introduce the relationship between current density J and the carrier densities n, and p. The current transport equations or current densities are defined by a combination of the drift and diffusion affects on the free carriers within the device. The drift portion of the current expression is simply due to the presence of an electric field, if there is one . This electric field creates an overall drift velocity v d , and is proportional to the electric field and the mobility of the carrier as shown in equation (2.4). (2.4) 7 PAGE 17 The mobility of the free carrier is a parameter that describes the ease of motion of an electron in response to an electric field and is affected by s everal things. The mobility model and how it i s used in the application of thi program is discussed shortly. The current density flowing in the direction of the electric field can be found by summing the product of the charge on each electron times its velocity over all electrons n per unit volume , or that i s II Jn.drift =I -qvi = -nqvd = nqj..i " E i=l The same procedure can be extended to the hole free carriers and then combined to give the total current density due to the applied electric field in equation (2.5). J drift = Jn,drift + J p.drift = ( qf-Ln n + qf..Lp P) (2.5) The other affect on the current density through a device i due to the diffusion of free carriers in the device if a special variation of carrier energies or densities exist. To investigate this affect we start by examining a portion of a semiconductor device that is of n-type and with an electron density that only varies in one dimension. The temperature across the region is average so that only the density varies with x and we then consider the number of electrons pas ing the plane at x=O per unit time per unit area (Figure 2.2). 8 PAGE 18 n -It 0 X Figure 2.2 ELECTRON CONCENTRATION DIFFUSION AFFECT Because the electron s are at finite temperature, they only have random thermal motion along one line even with no electric field applied. On average the number of electrons crossing the plane at x=O from the left side start at approximately x=-A after a collision. The lambda term is defined as the mean-free path for an electron, or that is the average distance an electron will travel before encountering another collision and is given by A= vrJ,r c n . Therefore the average rate per unit area of electrons crossing the plane at x=O depends upon the density of electrons that started at x=A and is N( -A.)v ,". The factor of '12 is due to the fact that half of the electrons move further to the left of x=A after a collision and half to the right. Similarly , the rate at which electrons pa s s the plane x=O from the right is N(A.)v," which gives a net rate of particle flow from the left i s . 9 PAGE 19 F = v,11 [ n( -A)-n(A)] which can be approximated at x=A by the Taylor series 2 expansion to give F v,. {[ n(O)-: A]-[ n(O)+: -v.,,A: and finally because each electron carries charge q, the particle flow corresponds to a current per unit area in equation (2 . 6). dn 1 = -qF = q/Lv,1 , d x n . diff (2.6) A more u s eful expression of equation (2.6) can be developed by using the theorem for equipartition of energy in this one dimensional case. [1] Thi s allows us to write m : = kT and A= v,11 r c,, and rearrange equation (2.6) into equation (2.7). 1 qD dn where, n . dijJ II dx (2.7) kT D , =-f.1, q The quantity D , i s known a s the diffu s ion constant or Einstein relation s hip. It relates the two important coefficients that characterize free carrier transport by drift and by diffu s ion in a solid. Now with the drift and diffusion terms derived, the total electron current density can be expressed by equation (2.8) iJn 1 , = 1 n . drifr + 1 n .diff = qJ.l,nE + qD , d X (2.8) Following the same procedure for holes leads to equation (2. 9) and combing the expresses the current den s ity for both carrier as in equation (2.1 0 ) 10 PAGE 20 (2.9) (2.1 0) Before we continue to the next section, the equations discussed so far are summarized in equations (2.11a-e) below which are often referred to as the equations of state [2], and a brief discussion on non-thermal equilibrium conditions is given. The current SDS program does not employ these techniques as solution mechanisms. Equations of State: 1 D dn n () If/ II -q II dX q f.LII dX dp 1 ()] p -=---+G-U dt q dX dn = _!_ . dl11 + G _ U dt q dX (2.11a) (2.llb) (2.11c) (2.11d) (2.1le) Not all devices operate under thermal equilibrium conditions and so techniques are developed which allow the analysis of semiconductors that are not at thermal equilibrium to occur. We call this quasi-Fermi analysis. Through this analysis the use of quasi-Fermi potential levels are established which relate the current to the gradient of the carrier at that given quasi-Fermi level. The expressions for the quasi-Fermi carrier densities and current densities are shown below in equations (2.12-2.15). 11 PAGE 21 1 -n dE! , -dlf/ J , " Jl, dx qJl , n dx 2.3. Carrier Densities The carrier models used in this program follow the Boltzmann statistics approximation of the Fermi-Dirac distribution model. The carrier densities are approximated in equations (2.16) and (2.17). N ( E J,-Ec) (q(lf/-f/J, ) ; n::::: c exp =n; e exp kT kT (2.12) (2.13) (2.14) (2.15) (2.16) (2.17) n;e is the effective intrinsic carrier concentration and has dependencies on temperature and total doping concentration. This project neglects band gap narrowing due to heavy doping levels but could be implemented as discussed in [3]*. The intrinsic carrier concentration then is only modeled as a function of temperature as shown below in equation (2.18). 12 PAGE 22 n;, (T) = N c (T)Nv (T) exp( ---'Eg::...._(T_)) 2kT (2.18) N c and N v are known as the density of states and along with the band-gap energy, E g also have dependencies on temperature as shown in equations (2.19-2 . 21). E (T) = E (300) +a[ 300 2 -_E._] g g 300+ fJ T + fJ 7i N c (T) = (_!_) 2 N c (300) 300 7i N v (T) = (_!_) 2 N v (300) 300 (2.19) (2.20) (2.21) N c (300) and N v (300) along with a and p, are predefined quantities given in section 2.1. 2.4. Mobility Models The mobility of an electron or hole to move throughout the silicon crystal lattice is representation of the responsiveness of a free carrier to move through the conduction and valence bands in the presence of an electric field. This phenomenon is governed by several complex interactions at the quantum mechanical level. [4] Temperature plays a major role because at higher temperatures the lattice will vibrate more causing the overall free carrier mobility to lower. Dopant levels also have a significant role since the dopant sites create local distortions in the lattice allowing the free carriers to scatter more. 13 PAGE 23 The models built into SDS allow for effective doping dependency and electric field dependency as shown in equation (2.22). There is also a factor included to take into account for high-injection levels in high power devices . These models can be turned on or off to account for any dependency that is desired. The constants and fitting parameters are tabulated below . Jl= flmax flmin -( 1=N (-X )----"J Ja""--+ flmin 1+ --+ f N ref np 1 2 ( 1 / Ja 1 where, Table 2.2 -Mobility parameters Nomenclature Program Variable Default Value Electron Mobilit y Concentration Nrefn 8.5 X 1016 Referenc e * Electron Mobility Alpha * alphan 0 . 72 Electron Mobility Beta * be tan 2 Electron Critical Electric Field* Ecn 8.0 X 103 Electron Maximum Mobility * mu_ maxn 1330 Electron Minimum Mobility * mu_ minn 65 Hole Mobility Concentration Nrefp 6.3 X J0J6 Referenc e * Hole Mobility Alpha * alphap 0.72 Hole Mobility Beta * be tap I Hole Critical Electric Field* Ecp 1.95 X 104 14 (2.22) Units cm-3 -V/cm cm2N-sec cm2N-sec cm-3 --V/cm PAGE 24 Table 2.2 (Con't.) Nomenclature Program Variable Default Value Units Hole Maximum Mobility * mu_maxp 495 cm2N-sec Hole Minimum Mobility * mu_minp 47.7 cm2N-sec * Constants to calculate field dependant mobility terms, when used. 2.5. Recombination and Generation Models The SDS program models both generation and recombination for use in determining current flow in the specified device when solving both the continuity equations and Poisson's equation. The most fundamental type of recombination defined as Shockley-Read Hall recombination . The recombination rates are defined in equation (2.23). 2 np-n. u =U =U = I " P r P ( n + n0 ) + r" ( p + p0 ) where, (E.-Eft] n o = n ; exp I kT I and (2.23) The electron and hole lifetime are defined by r" and rP, respectively and are defined below in equation (2.24) where (J' , is the capture cross-section and N, , the density of trapping centers. I r"=--(J'" v,,N, and (2 . 24) 15 PAGE 25 Another method to calculate the lifetime as used in this program is from [3] and is shown in equation (2 . 25) below. '(110 '(II = _ ____:.:_::...___ 1+ Ntnwl N SRH.II t"po and r = ----'-P 1 + Nrnral N SRH.p (2.25) The constants are tabulated below. Another type of recombination is Auger recombination, characterized by electron-hole pairs disappearing without the aid of trapping centers as in SRH recombination , and is mainly used when there are large high-injection rates. Auger recombination is defined by equation (2 . 27) . This model and the SRH model are added together when used in conjunction to provide the overall recombination rate as shown in equation (2.26) . To use the Auger model, the partial derivative terms used in Chapter 3 -Numerical Analysis will need significant modification ince Auger recombination is a function of both carrier densities. u = U SRH +U AUG (2 . 26) where, (2.27) The generation model used in this program models the basic generation affect due to impact ionization as defined in equation (2 . 28) . It is useful in characterizing reverse applied bias breakdown . In low electric field applied bias cases the generation rates are typically negligible. [4] The constants and fitting parameters are tabulated below. (2.28) where, 16 PAGE 26 Table 2 . 3 Ioni zation generation parameters Nomenclature Program Variable Default Value Units Electron Alpha * alpha _nO 2.25 X 107 Hole Alpha * alpha _pO 3 .08x 106 Electron Energy * E_nO 3.26 X 106 V/cm Hole Energy * E_pO 1.75 X 106 V/cm Fitting Parameter * m I -17 PAGE 27 3. Numerical Analysis In this chapter, we will proceed in setting up the numerical equations in order to solve the physical device equations established in the previous chapter. A procedure for setting up the equations in numerical form, linearizing those equations, and solving the subsequent block tri-diagonal matrix equation will be given in a step by step manner. However, first a brief introduction into some of the basic numerical analysis techniques required to solve the device equations will be given below. All of these techniques are used in transforming the differential device equations in to linearized numerical equations suitable for a numerical analysis solution. 3.1. Taylor Series Expansion The Taylor series is named after the English mathematician Brook Taylor (16851731) and the Maclaurin series is named in honor of the Scottish mathematician Colin Maclaurin (1698-1746) despite the fact that the Maclaurin series is just a special type of Taylor series . But the idea of representing particular functions as sums of power series goes back to Newton, and the general Taylor series was known to Scottish mathematician James Gregory in 1668 and to the Swiss mathematician John Bernoulli in the 1960s. Taylor was apparently unaware of the work of Gregory and Bernoulli when he published his discoveries on series in 1715 in his book Methodus incrementorum directa et inversa. Maclaurin series are named after Colin Maclaurin because he popularized them in his calculus textbook Treatise of Fluxions published in 1742. [5] 18 PAGE 28 As we have seen in Chapter 2Semiconductor Device Physics, the basic equations describing the physics of semiconductors are differential in nature and nonlinear. In order to perform numerical analysis on these equations, they must be linearized and in order to do this, the Taylor series expansion is performed about a point of interest. Additionally, for our equations of interest, the linearization must be done in each variable of interest (n, p , and '!') and this can be seen later in section 3.5. The general form of the Taylor series is show in equation (3.1) below and is referred to a the Taylor series expansion of the function fat a (or about a, or centered at a). [5] f(x) =I j < " l (a) (xa)" = f(a) +/(a) (x-a) n = O n! 1! +/(a) (x-a)2 + j "(a) (x-a)3 + ... 2! 3! (3.1) For the special case of a =0, the Taylor series becomes equation (3.2) and is given the special name Maclaurin series. In our analysis, we will be using the standard Taylor series expansion as the expansion point is not always zero. Since we are dealing with a system of nonlinear equations in which we are interested in multiple variables (n, p, and 'Jf), each equation that is linearized must be done so around each variable and so therefore the partial derivatives of each variable is introduced. Extending the above discussion on series expansion to a multiple variables is rather straightforward. We will also introduce the concept of perturbation state (or delta state). Let us re-examine equation (3.1) and transform it into the form show in equation (3.3) below. We assume that 'a' is a point such that f(a) = 0 and this point is defined as the equilibrium point of the system f(x) = f(x) since f'(x) = 0, when x =a. Since f(a) = 0, the nonlinear differential equation can be approximated near the equilibrium point by equation (3.4) . 19 PAGE 29 f(x) = f(a) + dx lx=a (x-a)+ higher order terms {(x)= dxlx=a(x-a) (3.3) (3.4) When the expansion point, in our case 'a', is sufficiently close to x, the higher order terms in equation (3.3) can be neglected as they are very close to zero. To complete the linearization , we define the perturbation state , or delta state (ox = x-a) and using the fact that of'(x) = f'(x) we obtain the linearized model show in equation (3.5). (3.5) This can be easily applied to multi-variable systems and is show in equation (3.6) . This is the form employed in section 3.5 to linearize the system. (3. 6) In summary, this technique allows us to linearize a system of nonlinear equations about a point of expansion. In this analysis , the concept is a more generic approach. As will be seen later, we linearize the current equations and generation and recombination equations as they involve nonlinear terms of an exponential nature. Since the three equations of interest are the two current continuity equations and Poisson's equation, and since the current continuity equations are differential expressions of the current equations , these must also be linearized. This allows the linearized current and recombination/generation equations to be ubstituted in directly. This will be more closely examined in section 3.5. 20 PAGE 30 3.2. Numerical Differentiation It is often necessary , e pecially in the procedures following, to estimate the derivative of a function and is one of the fundamental requirements in performing finite difference method techniques. The procedure used to perform this is called the difference method and is applicable to first or second order derivative . The procedure can be used to extend to higher order derivatives but is not di cussed here. Additionally, more elegant difference techniques can be used in order to lower the inherent numerical error associated with calculating the derivative of a function but is also not discussed here. Numerical differentiation is performed using a series expansion and in this instance a Taylor's series expansion is also used . If we assume that f(x) has as many continuous derivatives as necessary, then from Taylor's formula, equation (3. 7) is found, which is the same form as equation (3.1) but with the expansion about 'x', instead of 'a', which was arbitrary for the above discussion. [6] h . h2 ' h3 f(x+h)=f(x)+-f (x)+-f (x)+-f (x)+ ... 1! 2! 3! (3.7) One way to think of equation (3.7) is that we are expanding the function f(x) about x, within small increments of h . So if we replace the arbitrary 'a' in equation (3.1) with 'x+h', equation (3.7) is a straightforward application. If we solve equation (3.7) for f'(x) , the equation can be written as {(x) = f(x+h)-f(x) _!!:._ {(x)-!t_ { "(x)+ ... h 2! 3! and neglecting higher order terms gives equation (3.8), { (x) = f(x + f(x), for small value of h. (3.8) 21 PAGE 31 Equation (3.8) is called the forward divided difference approximation and graphically is the representation of the s lop of the line passing through points (x, f(x)) and (x+h, f(x+h)) . There are two other main approximations for a first order derivative and they are directly obtained from Taylor's series expansion as well. Assume we replace ' h ' with '-h ' in equation (3.7) and repeat the same process we obtain equation (3.9) below. Equation (3. 9) is referred to as the backward divided difference approximation . /(x)=f(x)-{(x-h)' forsmallvalueofh. (3.9) Finally , if we replace 'h' with '-h' in equation (3.7) and then subtract the resulting equation from the old one , we obtain /(x) = f(x+h)-f(x-h) _!t_ / " (x)-!!:.__ t PAGE 32 3.3. Numerical Integration Just as with derivatives, calculating the integrals of functions numerically is also of interest for this application. The main use for this technique, recalling the device continuity equations from Chapter 2, is to calculate the current density across the device using the recombination and generation terms that have previously been calculated. Since mesh spacing is al o defined, use of a simple trapezoidal integration rule is rather straightforward. Let us s uppo se we would lik e to evaluate the definite integral of the form show in equation (3.11) and which is graphically represented in Figure 3.1. I= {f(x)dx (3.11) y f X n 2 Xn 1 Xn=b X Figure 3 .1-INTEGRATION USING THE TRAPEZOIDAL RULE 23 PAGE 33 To calculate the integral (area under the curve) , we begin by s plitting up the interval [a, b] into n subintervals , each with uniform mesh spacing as shown in Figure 3.1. Thjs gives equation (3. 12) . b-a h = --and x i = a+ ih for i = 0,1, 2, ... ,n n (3.12) Now by examining each subinterval, the area under the curve can be approximated by the sum of the area s of the trapezoids. The area of the ith trapezoid is therefore with the sum or total area Tn being obtruned by extending over the entrre area. T,, ... + A n I + A,, h h = -[J(x0 ) + f(xl )] +-[f(xl) + f(xz )] + 2 2 h h + 2[f(x,_z) + f(x,_ l)] + 2[f(X,_1 ) + f(x,)] Finally , by collecting like terms we arrive at equation (3.13). Equation (3.13) is referred to as the composite trapezoid rule as the integral is estimated by a composite of the areas of n trapezoids. [6] 24 (3. 13) PAGE 34 3.4 . Tri-diago n a l Matrices This section will briefly discuss matrices with respect to diagonal dominate matrices and especially tri-diagonal matrices. A tri-diagonal matrix, as shown in Figure 3.2, is simply a matrix which has non-zero entities on the main diagonal , super-diagonal (upper), and sub-diagonal (lower) and frequently arises in the solution of numerical differential equations, such as in this project. aii ai2 0 0 0 a 2 I a22 a23 0 0 A= 0 a32 a33 a34 0 0 0 a43 a44 a4s 0 0 0 as4 as s Figure 3.2-TRI-DIAGONAL MATRIX The main content of thi section will discuss a recursive method for solving a block tri-diagonal matrix equation as this is what is employed in the solution of thi project. The algorithm is borrowed from [4] directly but is a standard algorithm used frequently. There are many method to olve these types of matrices , from the very basic naive Gaussian elimination, to more complex solutions, like LU decomposition, Jacobi method, etc. As we will see in Section 3.5 below , for a 1-D solution, a block tri-diagonal matrix solution must be obtained . The recursive algorithm for this procedure is as follows . Assume the matrix equation that needs to be solved is of the form in equation (3.14). A(N)by(N -1) + B(N)by(N) + C(N)by(N + 1) = F(N), 2 N L-1 (3. 14) for Oy(N) = [4J(N) , &z(N) , b"lf/(N)r , 2 N L -1 25 PAGE 35 where matrices A, B, and C are 3x3 dimension and vector F is of 3x 1 dimension. Since A, B, and C are 3x3 matrices we have a block tri-diagonal matrix as opposed to a simple tri-diagonal matrix. A direct solution would be to construct this matrix per equation (3.14) and use any of the aforementioned solution algorithms but this is computationally expensive. It would be much more efficient to only have to calculate the matrix blocks and solve for the delta state (or error value) and not have to manipulate every element in the matrix, including the majority of empty or zero elements. The use of the recursive algorithm allows for this computation savings. To continue, equation (3.14) is transformed into the following equation (3.15) involving unknown vectors at two points, N and N+ 1, instead of three. B . (N)8y(N) + c (N)8y(N + 1) = F ' (N) (3.15) Equation (3.15) can be solved for 8y(N) explicitly using basic matrix math, which yields equation (3.16) below. 8y(N) =B. (N)F' (N)-B ' (N)C. (N)8y(N + 1) (3.16) If we then take equation (3 .16) and reduce the division point N by one, we can substitute this result back into equation (3.14), rearranging the terms to arrive at equation (3.17). The following steps illustrate this procedure. 8y(N -1) =B. (N -l)F' (N -1)-B . (N -l)C. (N -1)8y(N) (division point N reduced by one) Substituting back into equation (3.14) yields , A(N)[ B . (N -I) F . (N -1)-B . (N -l)C' (N -1)8y(N)] +B(N)b'y(N) + C(N)b'y(N + 1) = F(N) Rearranging terms yields, [ B(N)-A(N)B. (N -l)C. (N -1) ]8y(N) + C(N)b'y(N + 1) = F(N)-A(N)B. (N -i)F. (N -1) 26 (3.17) PAGE 36 If we now compare equation (3.17) with equation (3.15) we obtain the recursive relations given below in equation (3.18), but for one less division point calculation . B . (N) = B(N)-A(N)B. (N -l)C. (N -1) c (N) = C(N) and F ' (N) = F(N)-A(N)B. (N -1)F. (N -1), 3$ N $ L -1 (3. 18) To finalize the algorithm since we would like to solve for 8y(N), the error values, we compare equation (3.18) with equation (3.14) for N=2 specifically and obtain equation (3.19). B . (2) = 8(2), c (2) = C(2) , F . (2) = F(2) (3. 19) We know that 8y(N-1) for N=2 is 8y(l), which is equal to zero. 8y(l) corresponds to the boundary condition error , but the boundary condition is fixed and so there is no error associated with that term and therefore reduces to zero. Finally, starting with equation (3.19) at N=2, then using equation (3.18) for the rest of the division points (N=3, 4, ... , L-1) we can use equation (3.16) to calculate the error values , 8y(N) , across the whole device (N=L-1 , L-2 , . .. , 2) starting with 8y(L)=0 , since this is also a fixed boundary condition. At this point, matrice s A, B, C , and Fare arbitrary with the dimensions defined above. As we will see in the next section , these matrices will be defined and so calculation of the above algorithm becomes very straightforward. Once A , B, C, and Fare calculated across the device using the recursive relations given, the error potential is solved for directly. The definition of these matrices comes from performing the finite element methods (difference method) on the device equations and then linearizing those equations. 27 PAGE 37 3 .5. P-N Junction DC Analysis In this section, a procedure for performing DC analysis of a 1 -D, P-N junction device will be developed. The basic P-N junction diode is examined as further development techniques for other devices follow the same basic procedures and solutions used in this section. For a 1-D Bipolar Junction Transistor, for example, only slight modification of this procedure is required and the resulting simulator program can be modified to use either solution depending upon the device specified. Another example of modification is time domain analysis, which requires extra division spacing but ultimately requires little modification to the base program with the conceptual idea that these further implementations can be readily incorporated at a later date. * Equations (3.20a-e) and Figure 3.1 show the basic setup for a 1-D P-N junction diode. The overall doping profile a of yet is undefined but could be a step junction, linear junction, and so on . __ D ap _ al.fl J p -q p ax qJ..L,, p ax ' (3.20a) an al.fl l, = qD, ax-qJ..L, n ax ' (3.20b) ap 1 aJ, -=---+G -U at q ax p P ' (3.20c) an= _ _!... aJ, + G _ U at q ax II II' (3.20d) (3.20e) For DC analysis the time derivatives in equations (3.20c and 3.20d) will become zero as the device is considered to be in steady state. 28 PAGE 38 X=O (N=1) p n Figure 3 .3BASIC P-N JUNCTION DIODE X=w (N=L) Before we can proce e d with numerical analysis of the device s how in Figure 2 , a few point s of intere s t s hould be noted. Many of the model s u s ed in the s olution procedure are of the ba s ic type , however can be expanded ea s ily to incorporate future models. For example , the intrin sic carrier concentration model is currently only temperature dependant and the mobility model i s fixed acro ss the device. Both could be modified to incorporate band gap narrowing for heavily doped device s or electric field dependency respectiv e ly. * More di s cu s sion on ea c h model u s ed in thi s ection along with it s capabilitie s i s further di s cu ss ed in Chapter 5 -Semiconductor Device Simulator and als o in Chapter 2. The recombination model employed in thi s s olution i s of the bas ic Shockley-Read Hall type. It is the mo s t fundamental type of recombination method and s o therefore i s most applicable to ba sic PN junction devices. Equation s 18C and d imply that the recombination and generation term are the net terms and can con s titute different models, such a s SRH and Auger recombination . In this ca s e , the net recombination rate from both model s i s used. Furthermore , the recombination centers only occupy a single energy level directly in the middle of the forbidden energy band. This implies that the equilibrium hole concentration i s equal to the equilibrium electron concentration which i s in turn equal to the intrin s ic carrier concentration , n o =p o=nj. For higher power device s , Auger recombination , a s discus s ed in Chapter 2 , can also be considered but i s not explicitly u s ed in this solution procedure. This allows us to u s e the basic SRH recombination equation from Chapter 2 , repeated a s hown here in equation (21 ). 29 PAGE 39 2 U =U =U = pn-n; p II ( r, p + n ; ) + r" (n+ n;) (3.21) If we only want to evaluate the device in low electric field applications technically the generation terms are not necessary and are not supported in some programs, such as in [3]. However, this program includes the generation terms and as with recombination we assume that the net hole generation terms are equal to the net electron generation terms. These terms are useful in supporting device with extremely high or reverse applied biases. We can examine the process of avalanche breakdown by incorporating these terms, however when low applied biases are used, the e terms are essentially zero. See more discussion in Chapter 5 on the generation function. Repeating the equations from Chapter 2 for the one dimensional case gives us equation (3.22) below. G = G p = G, = _l_(aPJl PI+ q where ( E1 , 0 J a P = a P0 exp lEI , (E,o J and a , = a ,0 exp lEf , E=-d'lf ax (3.22) Before we proceed to develop the numerical equations from equations (3.20) necessary to solve the device let us discuss the boundary conditions that are also necessary. Since we are performing a finite element analysis, an initial value or boundary value must be given. We will develop a two-point (Dirichlet type) boundary-value problem. Finite difference techniques can then be employed on the device physics equations and Taylor's series expansion to linearize the equations. This leaves a tri-diagonal block matrix that can be solved iteratively and finalized with Newton's iteration principle to arrive at a converged olution. We need to define the appropriate boundary conditions. With respect to Figure 3.3, p and n regions are located at x=O and x=w, respectively, where w is the device length. For our analysis, we assume that then-region is fixed at zero potential, or that is the external contact is grounded and the p-region is at any applied voltage. The thermal equilibrium analysis that follows also starts with the p-region fixed at ground . 30 PAGE 40 The boundary conditions need to satisfy a zero space charge region at each contact. Usually the semiconductor forming an ohmic contact with the metal has sufficiently high doping concentration. If there were a nonzero space charge region in such a highly doped semiconductor then there would be a correspondingly high electric field which would lead to breakdown conditions at the metal-semiconductor contact, so therefore the zero space charge condition must be satisfied. Future work can include the Shottky boundary conditions necessary to satisfy Shottky devices. * r(O) + p(O)-n(O) = 0, r(w) + p(w)-n(w) = 0 , p(O)n(O) = n;, p(w)n(w) = n;2 , 1/f(O) = V _ _.!._In[ p(O)J. 'lf(w) = V _ _.!._ ln[n(w)J. B n ; B n; where, B=_!L kT (3.23a) (3.23b) (3.23c) Equation (3.23a) sat isfies the zero space charge condition specified above, while equation (3.23b) sat isfie s the thermal equilibrium condition necessary at each contact. Equations (3.23a and 3.23b) are used to establish the two-point boundary conditions that are self-consistent among the two equations and the result is shown in equation (3.24). p(O) = _ {I +[l +( rr} n(w)={l+[l+(:(:JJ"} 2 n(O) = __!!j_ p(O) n 2 p(w)=-' n(w) (3.24) Conditions for the potential boundary conditions given in equation (3.23c) are deduced by the use of the quasi-Fermi potential s briefly di sc ussed in Chapter 2. These are re-written here in equation (3.25) for reference. (3.25) 31 PAGE 41 If the definition of the zero potential point on the n-region and applied bias on the p region is extended to equation (3.25) it follows that PAGE 42 (3.26b) (3.26c) (3.26d) (3.26e) The equations are now transformed into the difference equations. A brief discussion on equations (3.27a and 3.27b) will follow shortly (including the coefficients for the lambda terms and how they are obtained, but equations (3.27c-e) are obtained directly by using the difference methods employing the techniques discussed in Section 3 .2-Numerical Differentiation. 1 ,(M) = _q_[A,1 (M)p(N) + A ,2 (M)p(N + 1)] h(M) 1,(M) = _q_(A,1 (M)n(N) + A,2 (M)n(N + 1)] h(M) I where, A (M) = (M) 'lf(N)-'1/(N + 1) pi f.l, . ] -P PAGE 43 YJif/(N -1) + Y21f/(N) + Y31f/(N + 1) = _!l[f'(N) + p(N)-n(N)] c (3.27e) As you can see in the continuity equations (3.26c and 3.26d), the derivative of the current density with respect to the spatial variable, x, ha been transformed into the difference approximation in equations (3.27c and 3.27d). Also Poisson's equation in (3.26e) was transformed directly to equation (3.27e) also by using the difference method. This time, the second order central difference method was used since Poisson's equation is a second order differential equation. The transformation of equation (3.26e) into difference form is shown below in equation (3. 28). lf/(N + 1) -ljf(N) lf/(N) -ljf(N -1)] =-q (r(N) + (N)-n(N)] h (N) h(M) h(M -1) & p (3.28) Throughout the solution process, the mesh spacing is continually defined at each division point, whether it is a main division point or auxiliary point. Since we are using uniform mesh spacing and the auxiliary spacing is equal to the main spacing, equation (3.28) could be further transformed into the standard looking central difference equation for a second order derivative . Mesh spacing is left as a function of the division point o that future expansion of non-uniform mesh spacing is rather straightforward * . Lastly, if we examine equation (3.28) with equation (3.27e) the gamma coefficients are obtained directly. At this point in the solution process we have all of the equations transformed into numerical difference equations . The current density equations can be plugged into the current continuity equations and combined with Poisson's equation; the three equations can be then s olved simultaneously. However, the current density equations, as well as the generation and recombination terms, are not linear and so the application of Newton ' s iteration principle cannot yet be used. The next step in the solution process is to linearize these equations about the fundamental variables of interest, n , p, and 'I' 34 PAGE 44 However, before we proceed to linearize the equation s let us return to equations (3.27a and 3.27b), where if we compare these to equations (3.26a and 3.26b) it is obvious they were not directly obtained from the difference technique. Before continuing, please proceed to Appendix A, where a complete derivation and explanation of equations (3.27a and 3.27b) is given. Because the continuity equations contain expressions for the current densities which themselves contain derivatives of the fundamental variables, n, p, and \If, the derivatives for the current densities are defined at main points as well. Therefore, it is necessary to define the spacing between the auxiliary points shown in equation (29) . (3. 29) Equation (3. 30) shows the hole current density equation linearized through the use of Taylor's expansion theorem since it involves nonlinear functions. Since each equation is a function of one or more fundamental variables, the use of partial derivatives is applicable . 1 p(M) = a1 o (M) a1 o (M) " 8p(N)+ " 8p(N +1) I ap(N) ap(N + 1) (3.30) a1(M) a1 (M) + " 8 (N) + " 8 (N + 1) alf/(N) If/ alf/(N+1) If/ Solving for 81 P (M) gives a1 o(M) a1 o(M) 81 (M)"" P 8 (N) + P 8 (N + 1) P ap(N) p ap(N + 1) p a1 o (M) (Jj o (M) + " 8 (N) + P 8 (N + 1) a lf/(N) If/ a lf/(N + 1) If/ 35 PAGE 45 The hole current density is used as an example and the rest of the linearization process for the electron current density, recombination and generation equations is the exact same. Once the continuity equations are linearized, since they contain derivatives of the current densities , the delta state in equation (3.29) can be substituted in along with the rest shown below. aJ0(M) aJ0(M) 8JP(M)z P 8p(N)+ " 8p(N+l) ap(N) ap(N + 1) aJ0(M) aJ0(M) + P 8 (N)+ P 8 (N+1) alf/(N) If/ alf/(N+l) If/ aJ0(M -1) aJ0(M -1) 8J"(M-l)z P 8p(N-1)+ P 8p(N) ap(N -1) ap(N) aJ 0 (M -1) aJ 0 (M -1) + ' ' 81f/(N -1) + " 81f/(N) alf/(N-1) alf/(N) 8J 8n(N+l) " an(N) an(N + 1) + 8 (N)+ 8 (N + 1 ) alf/(N) If/ alf/(N+l) If/ 8J (M -I) z aJ 2 (M -l) 8n(N -1) + (M -1 ) 8n(N) " . an(N -1) an(N) + (M -1) 8 1f/(N _1 ) + (M -1) 8 1f/(N) alf/(N -1) alf/(N) 8G(N) z aco (N) 8 p(N -1) + aGo (N) 8n(N -1) + aco (N) 81f/(N -1) ap(N -1) an(N -1) a lf/(N -1) + aco (N) 8 (N) + aco (N) 8n(N) + aco (N) 8 (N) ap(N) p an(N) a lf/(N) If/ + aco (N) 8 p(N + 1) + aco (N) 8n(N +I)+ aco (N) 81f/(N + 1) ap(N + 1) an(N + 1) alf/(N + 1) 8U(N) z au\N) 8 (N)+ auo(N) 8n(N) ap(N) p an(N) 36 PAGE 46 Replacement of J r(M) by and J "(M -1) by (M -1) + aJ " (M -I) and so on in equation (25c and d) and rearranging the terms yields equation (3.31 ). 1 81 (M)-81 (M -1) " . " -8G(N)+8U(N) q h (N) J0(M) J0(M 1) =-..!_ " " -G0(N)+U0(N) q h (N) (3.31) ..!... 8Jn(M)-_8Jn(M -1) -8G(N)+5U(N) q h (N) = _.!__ J 2(M)-_J 2(M -1) -Go(N)+Uo(N) q h (N) We now have the continuity equations (3.31) linearized in terms of incremental variables . The incremental variables (delta state) are given above and these can now be substituted into equation (3. 31) which yields the two continuity equations involving only fundamental variable increments as unknowns. Due to the amount of physical space required to show the substitution and rearrangement process, the two equations are given directly in equation (3.32) and (3.33). It is a simple algebraic process but care must be taken that all signs are kept in their proper order and place . Continuity equation for holes ... . + 8 p(N -I) ... [ l (M -1) aG0 (N) ] qh (N) ap(N -1) ap(N -1) +{ . l _ -1)]-aG0(N) + aU0(N)}5 (N) ... qh (N) ap(N) ap(N) ap(N) ap(N) p +[ . 1 . _ aG0(N) ] 5p(N+I)aG0(N) . 5n(N-I) ... qh (N) ap(N + 1) ap(N + 1) an(N -1) -[aG0 (N) _avo (N)]. 5n(N) ... an(N) an(N) 37 PAGE 47 aco (N) . on(N + 1)-[ . 1 . (M -1) + aco (N) ] . O'lf(N -1) ... an(N + 1) qh (N) a 'lf(N -1) a 'lf(N -1) +{ . 1 _ -1)]aG0(N)}o N ... qh (N) a'lf(N) a'lf(N) a'lf(N) 'If( ) (3.32) [ 1 aG0(N) ] + qh.(N).a'lf(N+1) -a'lf(N+l) O'If(N+ 1 ) ... =-.1 [1(M)-J0(M -l)]+G0(N)-U0(N) qh (N) p P Continuity equation for electrons ... aG0(N) Op(N-l)+[aG0(N) _ aU0(N)]op(N) ... ap(N -1) ap(N) ap(N) aG0(N) + op(N + 1) ... ap(N + 1) +[-. 1 . (M -1) + aG0 (N) ] . On(N _1 ) ... qh (N) an(N -1) an(N -1) +{ . 1 _ -1)]+ aG0(N) _ aU0(N)}On(N) ... qh (N) an(N) an(N) an(N) an(N) +[ . 1 . + aG0(N) ]On(N+l) ... qh (N) an(N + 1) an(N + 1) +[-. 1 aG0(N) ]O'If(N-1 ) ... qh (N) a 'lf(N -1) a 'lf(N -1) +{ . 1 (N) ... qh (N) a 'lf(N) a 'lf(N) a 'lf(N) '1/ +[ . 1 . + aG0(N) ]O'If(N+1 ) ... qh(N) a'lf(N+l) a'lf(N+l) (3.33) =-, 1 [1 (M)1 (M -1)]G0 (N) + U0 (N) qh (N) II II 38 PAGE 48 Poisson's equation does not involve nonlinear functions and does not technically need to be linearized in incremental variables but for ease of calculation and to match the continuity equations we have done so in equation (3.34). This allows us to setup up a standard block tri-diagonal matrix . The linearization process is not repeated but is the same for the continuity equations . For example, p(N) is replaced by p0 (N) + 2Jp(N) and so on and then the terms rearranged and grouped. q t5 p(N)!ln(N) + r. (N)blf/(N -1) + y2 (N)blf/(N) + y3 (N)blf/(N + 1) c c =-q[r(N)+ p0(N)-n(N)J... c r. (N)blf/0 (N -l)-y2 (N)blf/0 (N -1)-y3 (N)blf/0 (N + l) (3.34) These three equations (3.32-3.34) are meant to be solved simultaneously and in order to do so we can transform the three equations into a matrix-vector representation that allows us to employ the recursive algorithm discussed earlier. It will also allow us to develop the coefficients of each matrix. Equation (3.35) is the matrix-vector representation of the two continuity equations and Poisson's equation combined . A(N)t5 y(N -1) + B(N)t5 y(N) + C(N)t5y(N + 1) = F(N), (3.35) for 2 N L-1 To evaluate the matrix coefficients A , B , C, and F let us examine equation (3.35) in matrix form for three division points , N-1, N , and N+1 as shown in Figure 3 . 5 and then expand each equation in the matrix to form equations (36a-c). 39 PAGE 49 la1 1(N) a1 2(N) a13(N) b11(N) b12(N) h1 3(N) C11(N) C 1 2(N) C13(N)l a2 1 (N) a22 (N) a23 (N) b21 (N) b22 (N) b23 (N) c21 (N) c22 (N) c23 (N) a31 (N) (N) a33 (N) b31 (N) b32 ( N ) b33 (N) c31 (N) c32 (N) c3 3 (N) p(N-1) 8n(N -1) l5f//(N-I ) 8 p(N) l5n(N) l5f//(N) l5p(N+1) 8n(N +I) l5f//(N +I) Figure 3.5-MATRIX FORM OF THE MATRIX-VECTOR EQUATION a11 (N)J p(N -1) + a1 2 (N)Jn(N -1) + a1 3 (N -1)Jf//(N -1) +b1115 p(N) + b1 2 (N)Jn(N) + b1 3 (N)I5fi/(N) +c11J p(N + 1) + c1 215n(N + 1) + c13Jf//(N + 1) = 1 (N) a2 1 (N)J p(N -1) + a22 (N)Jn(N -1) + a23 (N -I)I5fi/(N -1) +b21 15 p( N) + b22 ( N)Jn( N) + b23 ( N)l5fl/( N) +c21J p(N + 1) + c22Jn(N + 1) + c23Jf//(N + 1) = (N) (N)J p(N -1) + (N)Jn(N -1) + (N -I)Jfi/(N -I) +b3 1J p(N) + b32 (N)Jn(N) + b33 (N)I5fi/(N) +c3 1J p(N + 1) + c3215n(N + 1) + c3 3Jf//(N + 1) = (N) (3.36a) (3.36b) (3.36c) Now if we compare equations (3.36a-c) with equations (3.32-3.34) respectively, the matrix coefficients for A, B, C, and Fare obtained directly . The matrix coefficients are tabulated below for easy reference in Figures 3.6-3.9 40 PAGE 50 1 -1) aG0(N) a" =-. . -' qh (N) ap(N -1) ap(N -1) aG0(N) a = -----'-----'---1 2 an(N -1) 1 -1) ()G0(N) al3 = ' qh(N) af//(N-1) af//(N-1) aG0(N) a = -----'-----'---2 1 op(N -1) 1 -1) aG0(N) a22 = + _ ____:___.:____ qh (N) an(N -1) an(N -1) 1 a1 c M o ac o c N) a23 = + ------'----'--qh (N) a f//(N -1) a f//(N -1) 1 a 3 1 = 0, a32 = 0, a33 = Y1 (N) = h(M -l)h (N) Figure 3.6-COEFFICIENT A OF MATRIX-VECTOR EQUATION b = 1 II qh .(N) ap(N) ap(N) ap(N) ap(N) b = _ aG0(N) + aU0(N) 1 2 an(N) an(N) b = 1 _ -1)]-aG0(N) 1 3 qh . (N) a f//(N) a f//(N) a f//(N) b = aG0 (N) _ auo (N) 2 1 ap(N) ap(N) b = 1 22 qh. (N) an(N) an(N) an(N) on(N) b = 1 _ -1)]+ aG0(N) 23 qh . (N) af//(N) af//(N) af//(N) Figure 3.7-COEFFICIENT B OF MATRIX-VECTOR EQUATION 41 PAGE 51 1 dG0(N) ell= . . ' qh (N) dp(N + 1) dp(N + 1) dG0(N) c = -_ ____:_____:__ 1 2 dn(N + 1) 1 (M) CJG0 (N) c = . -__ ;_____:__ 1 3 qh . (N) d lf/(N + 1) d lf/(N + 1)' dG0(N) c =--2 1 dp(N + 1) 1 (M) dG0 (N) c = . +---22 qh . (N) dn(N + 1) dn(N +I) I (M) ()G0 (N) C23 = + __ ;_____:__ qh (N) d lf/(N + 1) d lf/(N + 1) 1 c 3 1 = 0 , C32 = 0, c33 = Y J (N) = h(M)h (N) Figure 3.8-COEFFICIENT C OF MATRIXVECTOR EQUATION ft1 =-.1 -l)]+G0(N)-U0(N) qh (N) ft2 = .1 -I)]-G0(N)+U0(N) qh (N) ft3 =q[r(N)+ p0(N)-n(N)J. .. c -]'j (N)If/0 (N -1)-y2 (N)If/0 (N)-y3 (N)If/0 (N + l) Figure 3.9-COEFFICIENT F OF MATRIXVECTOR EQUATION Now that we have the coefficients for the matrix-vector equation determined they can be calculated aero s the whole device in order to solve equation (3.35) . However, as we have seen the recursive algorithm is more efficient and requires less computations then solving equation (3.35) directly. Refer to section 3.4 for the recursive algorithm and appropriate relationships. Before the coefficients for equation (3.35) can be evaluated the partial derivatives embedded in the coefficients must be determined first. The partial derivatives are determined directly from the current density, recombination, and generation equations 42 PAGE 52 (3.27), (3.21), and (3.22) respectively. The derivation techniques for several of the partial derivatives are given in Appendix B, you may refer to that now for more information. The required partial derivatives are given in Table 3 . 1 below for reference. Their calculations are given below in summarized form. No equations numbers are associated with the calculations; they must be referenced to the table numbers below. It may seem excessive but the number of partial derivatives is not as extensive as the list appears. Most of the derivatives can be coded into one function as we will see in Chapter 5. The first table below simply lists all the partial derivatives required to support the matrix coefficients A , B , C , and F. Table 3.1 -Partial derivatives required for matrix coefficients (M -1) (M -1) (M) (M) (Ia) (1b) (lc ) (1d) ap(N -1) ap(N) ap(N) ap(N + 1) (M -1) (M -1) (M) (M) (2a) (2b) (2c) (2d) alf/(N-1) alf/(N) alf/(N) alf/(N+1) (M -1) (M -1) (M) (3a) (3b) (3c) an(N -1) an(N) an(N) a12 (M -1) a12 (M -1) (M) (4a) (4b) (4c) a lf/(N1) a lf(N) a lf/(N) an(N + 1) alf/(N + 1) aU0(N) an(N) aG0(N) (Sb) ap(N -1) (6a) aGo (N) (6b) aGo (N) ap(N) ap(N + 1) aU0(N) ap(N) aG0(N) (Sa) an(N -1) aG0(N) alf/(N -1) (7a) aG o(N) (7b) an(N) aG0(N) (8a) (8b) a f//(N) aG0(N) an(N + 1) aG0(N) alf/(N + 1) (3d) (4d) (6c) (7c) (8c) Referring to the current density equation from Appendix A at division point M-land M repeated below in equation (3. 37), the partial derivatives for (la-d) are obtained directly. 43 PAGE 53 J P(M -1) = q [ A P1(M -l)p(N -l)+A 2(M -1)p(N)J h(M -1) P J P(M) =-q-[ A1 , 1(M)p(N)+ A p2(M)p(N + 1)] h(M) -1) = q A (M -1) dp(N -1) h(M -1) p i ' (la) (M -1) = q A (M -1) dp(N) h(M -1) P2 ' (lb) __ q_A (M) dp(N) -h(M) p i ' (lc) q ---A (M) dp(N + 1) -h(M) ' '2 ' (1d) (3.37) The partial derivative of the current density with respect to the potential, 'I' is more complicated and requires use of the product rule and chain rule for derivatives since the lambda terms are functions of the potential themselves. The results are tabulated below. The lambda terms calculated as a part of the current density equations are repeated here for reference in equation (3.38) . The partial derivatives for (2a-b) are tabulated below. A (M)= (M)'If(N)-'!f(N+l) p i J1p 1 P < M ) ' -e A (M)= (M)'If(N)-'!f(N+l) p2 J1p 1 P ( M ) ' -e A (M -1) = J1 (M -1) 'lf(N -l)'lf(N) (3.38) p i p l-eP < M l ' A (M -1) = J1 (M -1) 'lf(N -l)'lf(N) where p2 p 1 _ eP< M ) ' {J(M) = 8 [ 'lf(N)-'lf(N + 1)] {J(M -1) = O['lf(N -1)'lf(N)] 44 PAGE 54 q { aA-P 1(M-1) aA-P 2(M-1)} af/I(N -1) h(M -1) p(N -1 ) af/I(N -1) + p(N) af/I(N -1) (2a) where, aA.P 1 (M -1) _ f.lp ( M -1) [ 1 /](M -1)eP < M I > ] alf/(N-l) (1-e P < M 1 ) ) (1-eP < M -1)) aA.P2 (M -1) = f.lp (M -1) [I+ /](M -1)eP < M t > ] a lf/(N -1) ( 1eP< M 1 ) ) ( 1eP< M 1 ) ) (M -1) (M -1) =-'----af/I(N) alf(N -1) (2b) where, aA-p1 (M -1) aA-p, (M -1) = a lf/(N) a lf/(N -1) aA-p2 (M -1) = _ aA.p2 (M -1) a lf/(N) a lf/(N1) =-q-{ (N) +p(N+I) aA-p2(M)} (2c) a If( N) h( M) P a If( N) a If( N) where, aA-P 1 (M) _ f.1/M) [1 -fJ(M)e P < M l ] alf/(N) (t-eP < M > ) (t-eP < M l ) p 2 -p 1 + ..:....,-.:-___:__----,-aA. (M) f.1 (M) [ fJ(M)e P < M > ] alf/(N) (1-eP < M > ) (t-eP< M l ) 45 PAGE 55 d/o(M) = p (2d) alfi(N + 1) alfi(N) where, aA"1(M) = aA"1 (M) al!f(N + 1) al!f(N) aA " 2 (M) = al!f(N + 1) al!f(N) The procedure for evaluating the electron related partials derivatives in 3a-d and 4a-d are exactly the same. The partial derivatives for these terms are tabulated below. -1) = q A (M -1) an(N -1) h(M -1) "1 ' (3a) -1)-q A (M -1) an(N) -h(M -1) " 2 ' (3b) __ q_A M an(N) -h(M) "1 ( ), (3c) _ _ q_A (M) an(N +1) -h(M) "2 ' (3d) at,: (M -1) q { (N 1 ) (M -1) + (N) a-1,12 (M -1)} ( 4a) alfi(N-1) h(M-1) n alfi(N-1) n alfi(N-1) where, aA111 (M -1) f.111 (M -1) aA"2 (M -1) ----"..:......_ __ = ----"---____;_ __ _ a 1/f(N -1) f.1" (M -1) a 1/f(N -1) aA"2 (M -1) f.1" (M -1) aAP, (M -1) ----"-=---= . ____:. __ _ al/f(N-1) f.l"(M-1) al/f(N-1) -1) = q [n(N _ 1 ) aA111 (M -1) + n(N) aA"2 (M -1)] ( 4 b) al!f(N) h{M-1) al!f(N) al!f(N) 46 PAGE 56 where, aA-,1 (M -1) aA-,1 (M -1) = JL, (M -1) = aA.p2 (M -1) a lf/(N) alf/(N -1) JLP(M -1) alf/(N -1) aA-,2 (M -1) = aA-,2 (M -1) = JL, (M -1) aA."1 (M -I) alf/(N) alf/(N -1) JLP (M -1) alf/(N -1) ai,:(M) =-q-{ n(N) +n(N+1) (4c) a If!( N) h( M) a If!( N) a If!( N) where , aA.,1 (M) = JL, (M) aA.p2 (M) alf/(N) JLp (M) a lf/(N) aA.,2 (M) = JL, (M) aA."1 (M) alf/(N) JL"(M) alj/(N) =-q-[n(N) aA.,1(M) +n(N +1 ) aA.,2(M)] (4d) alf/(N+l) h(M) alf/(N+l) alf/(N+1) where , aA-,1 (M) = aA.,1 (M) = JL, (M) aA.p2 (M) alf/(N + 1) a lf/(N) JLp(M) alf/(N) aA.,2 (M) = aA.,2 (M) = JL, (M) aA.p1 (M) alf/(N + 1) a lf/(N) JL"(M) a lf/(N) The amount of partial derivatives for the current density equations may seem excessive, but as will be explained in Chapter 5 , only one function for the lambda coefficients and partial derivative lambda expressions are required to calculate all of the above which allow s the calculation for each of the matrix coefficients to be more efficient. The partial derivatives for the recombination and generation terms are also given below . Their derivations are similar and briefly di cussed in Appendix B as well . 47 PAGE 57 Recombination terms: dU(N) n(N)-r,p(N) dp(N) r, ( p(N) + n ; ) + r, (n(N) + n ; ) (Sa) dU(N) p(N)+r,U(N) dn(N) r, (p(N) + n;) + r, ( n(N) + n;) (Sb) Generation terms: dG(N) = + rna A (M -1) a (N) (6a) dp(N -1) -h(M -1) "1 " dG(N)=+[ rna A (M)]a (N) (6b) dp(N) -h(M -1) "2 h(M) pi " dG(N) (M)a N dp(N + 1) -h(M) P 2 " ( ) (6c) where, h(M) rn = ------:-----'---" 2h .(N)' h(M -1) rn = ----'-.,------'b 2h . (N) a,(N)=a0exp[ IE,, 1 ] and E(N)=rn"E(M-l)+rnbE(M) " E(N) The choice of the sign for the partial derivative terms depends upon the sign of the current density at the given mesh point, that is J , (N) > 0 or JP (N) < 0 determines +or-respectively when, 48 PAGE 58 (JG(N) =+ rna A (M -1) a (N) (7a) dn(N -1) -h(M -1) "1 " dG(N) [ rna rnb ] -dn_(_N_) = h(M -1) An2(M -1)+-h(_M_) Ani(M) a"(N) (7b) (JG(N) (M)a (N) (7c) dn(N + 1) -h(M) "2 " h(M) rn =-,--a 2h .(N)' h(M -1) rn =---b 2h .(N) The choice of the sign for the partial derivative terms again depends upon the sign of the current density at the given mesh point, that is i" (N) > 0 or i " (N) < 0 determines +or-respectively when , The choice of the sign or + for the remaining partial derivatives below depends upon whether i" (N), i P (N), E(N) > 0 or i" (N),i,, (N), E(N) < 0 respectively. For example, G"'11 is positive is E(N) > 0 but negative if E(N) < 0, however G"'3 1 is negative is E(N) > 0 but positive if E(N) < 0. The term that governs the sign follows each Gvrij below. 49 PAGE 59 aG(N) a Vf(N -1) = G(/111 + G(/112 + G(/11 3 + G(/114 (8a) aG(N) avr(N) = G(/121 + G(/122 + G(/123 + G(/124 (8b) aG(N) avr(N + 1) = G(/131 + G(/132 + G(/133 + G(/134 (8c) where, G(/111=+ rna a(N) EpO .11p(N)I h(M -1) P E(N)2 q , E(N) G =+ m a ap(N) aJP(M -1) (/112 -q . avr(N -1)' JP(N) G(/113 = rna a (N) Eno 1111 (N)I h(M -1) " E(N)2 . q , E(N) G = + m"a" (N) a1, (M -1) (1114 . J (N q a vrCN -t) , , ) G -m b E II (N)I (1131 =+--a (N) P 0 â€¢ P h(M) P E(N)2 -q------', E(N) G _ -+ m b ap(N) aJ (M) (1132 -p J q . avr(N)' P(N) G(/133=+____!!!!!__a,(N) E,o .II,(N)I h(M) E(N)2 q , E(N) G =+mb a , (N) aJ, (M) (1134 q . a Vf(N) ' 1 , (N) 50 PAGE 60 Jn(N) With aJI of the coefficients for the matrix-vector form of the continuity equations and Poisson's equation tabulated and the partial derivatives for those coefficients calculated, we can u se them to develop an algorithm which calculates the thermal equilibrium or applied DC bias potential and carrier densities throughout the device. 3.6. Numerical Solution This section will discuss the numerical solution process used to solve the matrix vector equation (3.35), including development of the trial values, or initial guess, and finally Newton's iteration principle to update the trial values and check for convergence of the s olution . A solution process is also discussed which was used to develop the program given in Appendix C. This s olution algorithm involves solving the two continuity equations and Poisson's equation as previou s ly discu ssed. The proces s of transforming the se equations into finite difference equations allows the use of the implicit method as a solution procedure which is si mply a multidimensional Newton iteration problem. Equations (3.32-3 . 34) have been developed s uch that full coupling exists between the three fundamental variables, n , p , and 'I' or that is each fundamental variable is allowed to change during each iteration. In thi s manner the equations can be solved simultaneously in a straight forward fashion. Other methods exist to solve the numerical equation developed including decoupled method s (Gummel method [3]) but are not di sc ussed here. They have their various differences and benefits for s pecialized situations but for a generic solution algorithm for most po ss ible device conditions (low injection vs. high injection for example) the couple Newton iteration 51 PAGE 61 procedure is preferred. It will increase the overall number of computations per iteration but decrease the overall computation time for convergence. The solution to equation (3.35), repeated below for reference and show in matrix form, can involve a rather significant amount of matrix elements. An efficient solution procedure is required or the computation time, even for simple devices, can become excessive and limits the use of such simulator tools to any practical sense. Use of the coupled Newton iteration principle with the recursive solution procedure will suit the needs of this program adequately. A(N)o y( N -1) + B(N)o y(N) + C(N)o y(N + 1) = F(N), for 2 N L -1 where, [p(N)] [ Jp(N)] y(N) = n(N) , 8 y(N) = on(N) 1/f(N) 01/f(N) B(2) C(2) 0 0 o y(2) F(2) A(3) B(3) C(3) 0 o y(3) F(3) 0 A(4) B(4) C(4) oy(4) F(4) = 0 0 0 A(L-2) B(L-2) C(L-2) oy(L-2) F(L-2) 0 0 A(L-1) B(L-1) o y( L -1) F(L -1) Newton's iteration principle requires that we have an initial guess or starting value to calculate the error s tate for the first iteration . The calculated error is then added to the initial value and the process repeated until convergence is reached. All devices, whether simple 1 D P-N junctions or elaborate 2-D bipolar devices , should have the thermal equilibrium condition solved for first, or that i s no external applied bias . Not only does it allow for faster computation speeds on the applied bias conditions but allows the user to obtain a good feel for whether the device's equilibrium conditions are being met as intended. 52 PAGE 62 For the thermal equilibrium condition, trial potentials are simply provided based on the charge neutral assumption that the boundary conditions were established under. 2 n . Carrier densities are given by PT = -r, = __ , ' where r is the effective doping r concentration throughout the device. The zero-bias trial values for the potential are given by equation (3.39) below . '1fT = 1 [ n. ] BIn , p-region 1 [r(x)] . -In --, n-regron B n ; (3.39) Now that the trial potentials are defined we can conceptualize a solution procedure, which is given below . A flow chart representation is provided in [4] as well as the author's version in Chapter 5 , which represents this program's solution. (1) Define physical and numerical parameter (2) Define iteration limit and set current iteration number to zero (3) Use the matrix-vector coefficients, A, B, C, and F to calculate the error for each iteration u ing the recursive algorithm discus ed in Section 3.4. This is equivalent to solving equation (33) using direct Gaussian elimination for example. (4) Improve the initial guess or trial potentials by adding the calculated error. (5) Check for convergence against defined convergence tolerance (6) End program if solution has converged otherwise update current iteration limit by 1 and return to step (3) repeating the process until solution has converged or iteration limit is reached . Assuming that a converged solution for the equilibrium condition has been reached, extending the solution to an external biased condition only requires updating the equilibrium condition potential with a simple proportional relationship and then repeating the process defined above . Technically updating the potential is not necessary but will lower the total iteration number for the DC condition. This relationship is given in equation (3.40) and illustrated in Figure 3.10 below. 53 PAGE 63 .2 G) 0 Q._ x=O 'I', ( I ) )( )(:Q ( =L) (3.40) Figure 3.10-UPDATED POTENTIAL BASED ON PREVIOUS POTENTIAL [ 4] A tabulated list of phy ical constants and user defined constants as well as a tabulated list of numerical constants required to develop the program in Appendix C is given below. Table 3.2-Numerical constants Nomenclature Program Variable Default Value Units Device Length length 10 urn Convergence Tolerance etol I x 10-6 -# Grid Points L 1001 Iteration Number it 0 -Iteration Limit itl 15 -54 PAGE 64 Table 3 . 3 Phy sical constants Nomenclature Program Variable Default Value Units Dielectric Permittivity eo 8.854 X 10-14 F/cm Constant Electronic Charge q 1.602x 10-19 C (A-sec) Boltzmann's Constant k 8 . 62 X 10-5 eV/K Thermal Voltage Vt 0.02586 v Boltzmann Factor INt 38.6 IN (Thermal Voltage) Temperature T 300 K Donor Concentration Nd I x 1015 cm-3 Acceptor Concentration Na I x 1015 cm-3 Left Contact VI 0 v Right Contact Vr 0 v 55 PAGE 65 4. SDS Results Investigated In this chapter, we will review the output of the SDS program and compare it to analytical results a s well as to other simulators, notable SimWindows . [7] Sim Windows is also a 1-D semiconductor device simulator capable of simulating the same devices as the SDS program. The first section will compare the SDS results with respect to standard analytical models, for example the built-in potential with respect to using the depletion approximation as opposed to solving the numerical solution which makes no such assumptions. The purpose of this is to verify the working operation of the program versus general semiconductor operation that is learned at the classroom level. 4.1. SDS Output Compared to Analytical Assumptions Most basic PN Junction analysis i s performed using the analytical device expressions as defined in Chapter 2. To make these differential equations mathematically tractable , assumptions must be made in order to arrive at a solution. One of the main assumptions u se d in this approximation is the Depletion Approximation. [1] The Depletion Approximation is a way to simplify Pois s on's equation such that the potential throughout the device can be calculated without the need of the carrier concentrations , n and p. The approximation is performed by assuming that the transition region between the n and p bulk regions, or that is depletion region , is totally depleted of mobile carriers and that the mobile majority carrier densities abruptly become equal to corresponding dopant concentrations at the edges of the depletion region. The charge density is therefore zero everywhere except in the depletion region , where it equals the ionized dopant concentration. Additionally , we generally assume that the dopant atoms are complete (100%) ionized and free to contribute to current conduction . This reduces Poisson's equation to what is shown in equation (4.1) below . 56 PAGE 66 (4.1) If we further only evaluate Poisson's equation with respect to then-region of the PN junction, equation (4.1) reduces to equation (4.2) below . This reflects no acceptor atoms in then-region and can easily be integrated from an arbitrary point in thenregion depletion region to the edge of the region resulting in equation (4.3). Similar results are obtained in the p-region in equation (4.4). a2 a qNd --=--=---(4 . 2) (4 . 3) (4.4) If we further integrate equation (4.3) and (4.4) one more time, we arrive at expressions for the potential in each respective nand p-region shown in equations (4.5) and (4.6) . qNd 2 t/J(x)= " --2-(xn -x) O PAGE 67 The total potential change from one region to the other is referred to as the built-in potential and is the potential whereby at equilibrium it balances the drift and diffusion components across the junction. The expression for the built-in potential is shown below in equation (4.8). (4.8) With these expressions developed from an analytical approach using the Depletion Approximation we can evaluate the SDS program solution with respect to the analytical solution. The equations derived above generally result in the graphs shown in Figure 4.1. These graphs apply to an abrupt step junction which is the default device loaded into SDS. Plot (a) shows the step junction profile, with one dopant type per region. This creates a step junction trial potential for use in the numerical solution . Plot (b) is simply showing the depletion approximation, or that is , the depletion region carrier density is negligible. Plot (c) is the charge density associated with the dopant atoms. Plot (d) is the electric field in the device, while plot (d) is the potential in the device. 58 PAGE 68 n-p n-type X X p-type (a) (b) p E X X X n (c) (d) qln X (e) Figure 4.1 -ABRUPT STE P JUNCTION PROPERTIES From plot (d) in Figure 4.1 the potential is only zero at t h e transition region w h en the doping profiles on each side are equal. When they are not , the depletion region will push further into the region that has less doping. 5 9 PAGE 69 Referencing plot (d) and equations (7) and (8), we can calculate the built in potential and potential in then and p-regions only to see if they compare to our numerical solution. Let us assume that the dopant profile in each region is equal which gives us , = kT In N d = 1.38x10-23 . 300 . In lx1015 = 0 . 2878 V and q n; 1.602xl0 1 9 (1.45x1010 t = _ kT In N a = _ 1.38xl0 -23 300 . In lxl 0 1 5 = -Q_2878 V P q n ; 1.602xl0 1 9 (1.45x1010 t f/J; = , _ = kT In N d N a = 1.38xl o -23 300 In lxl01 5 lxl 0 1 5 = 0 _ 5756 V P q n;2 1.602x10 1 9 (1.45xl Olo ) 2 Figure 4.2 shows the output for the potential of the SDS program when run with an abrupt step junction of equal doping in both regions. The cursor in the graph is displaying the pertinent information. potential-02 0 1 .01 -02 I 1-----------------/ (j)n I( 7.S1e-008 v 02818 ..03 IC.1â€¢3e-OO& y .. ,.,. (j)j= PAGE 70 We can also observe in Figure 4.2 that the potential is zero at the junction transition as we expect since the doping densities are equal in the bulk n and p-regions. Figure 4.3 shows the equilibrium electric field throughout the device. It is also symmetrical about the junction as we expect. Energy-euilibrium llllllo!---..,Jo ,,.----="o2,...--0,L-3----,oLâ€¢ -----,-o'o-s ----='o s,...--""'o 7------,o"'"s --------oo'=-9 ---', ONce length (1t) 1r1 urn x 1o-5 Figure 4.3 -EQUILIBRIUM ELECTRIC FIELD Figures 4.4 and 4 . 5 show the potential for a PN junction where the acceptor dopant density is greater than the donor density. This shows the depletion region extending further into then-region as it should, as well as a nonzero value at the transition region . Figure 4.5 is taken from Sim Windows as verification the correct solution was obtained. The only difference in SimWindows is that all plots are scaled to zero, meaning there are no negative potential values . SDS has an option to do this as well and the result is shown in Figure 4.6. The electric field is always max at the junction transition ; however when the dopant densities are unequal across the junction, the field will be push further into the region with less doping atoms just as the depletion region does. 61 PAGE 71 potential equilibrium 02 0 I .{)2 .()3 .{)4 Oevu:e Length (x) 1n m.crons I( 10 .. Figure 4.4 DEVICE POTENTIAL WITH UNEQUAL DOPING DENSITIES PotM 1.00 0 . 90 0 . 80 0 . 70 0 . 60 0 . 50 0 . 4 0 0 . 30 0 . 20 0 .10 2 . 0e-16 0 . 00 1.00 2 . 00 3 . 00 4 . 00 5.00 6 . 00 7 . 00 8 . 00 G r id Position (micron s ) Figure 4.5 SIMWINDOWS COMPARISON TO FIGURE 4.4 62 PAGE 72 09 08 07 06 i OS !i oâ€¢ 0 3 0 2 0 I 00 0 I potential-equilibrium 0 2 03 O A 0 5 0 6 0 7 08 Figure 4.6 SDS OUTPUT SHIFfED TO ZERO oâ€¢ , â€¢ 10-d Figure 4.7 below i s the SimWindow s equilibrium electric field with unequal doping den s ities. This compare s to the SDS output in Figure 4 .8. Field (V/cm) 0 .00--,--------------------., -5.0e+OJ 2 .0e+04 2 . 5e+04 3 .0e+04 3 .5e+04 4 .0e+04 0 .00 1 .00 2 .00 3 .00 4.00 5 .00 6 .00 7 .00 8 .00 9 .00 10.00 Grid Position ( m icrons ) Figure 4.7SIMWINDOWS ELECTRIC FIELD 63 PAGE 73 Energy euilibrium -1 3 08\'lce length (11) ttl m 1cron s â€¢ 104 Figure 4.8 SDS ELECTRIC FIELD Finally the abrupt step junction is evaluated with an external applied bias of 0.5 volts at the left contact. This forward biases the PN junction as we would expect. A 0.5 volt signal is also applied at the right contact to reverse bias the junction. Alternatively, a negative voltage could be applied to the left contact to achieve the same affect. The doping densities are kept unequal as in the previous examples. Figure4 . 9 shows the device potential when the device is forwarded biased. From classic analytical theory we expect the built-in potential to be lowered as is shown in the graph. Figure 4.10 is a Sim Windows comparison of the same output. While the device is forwarded biased we can also calculate the current density through the device and compare it to Sim Windows output as well. One final experiment not shown here is to reverse bias the device and slowly increase the bias 64 PAGE 74 until a breakdown condition is achieved. This better shows the effect of the generation dependent models as a high bias tends to achieve impact ionization more readily then lower level electric fields. potential applied a9 a8 a7 as 05 t c aâ€¢ a3 a2 a I oa 01 02 a3 aâ€¢ as ao a7 08 09 length (x) tn mtcrons X 10-6 Figure 4.9-SDS DEVICE POTENTIAL WITH EXTERNAL APPLIED BIAS 65 PAGE 75 PotM 1.00 0 .90 0 . 80 0.70 0.60 0.50 0 . 40 0 . 30 0 . 20 0.10 0 . 00 0 . 00 1.00 2 .00 3.00 4 . 00 5 . 00 6 .00 7 . 00 8 . 00 9.00 10 . 00 Grid Position [microns) Figure 4 .10-SIMWINDOWS POTENTIAL WITH EXTERNAL APPLIED BIAS 66 PAGE 76 5. Semiconductor Device Simulator In this chapter, a brief overview of the SDS program is given with explanations for key aspects of the program that are not explicitly covered in any other section. We will begin by looking at the "main" program and then investigate the functions that are implemented in support of the main program. As we have seen throughout the project, the main purpose is to write a program that solves the matrix-vector equation which was established from the classic differential equations in Chapter 2 . In most aspects the program can be developed however the developer likes, but in certain instances poor coding techniques will lead to longer computation times and ultimately slower convergence speeds for a given device setup. The main important goal is to manipulate the matrix coefficients efficiently and to solve the matrix-vector equation with as little iteration as possible. On key element to this is establishing the proper boundary conditions and trial values. The fust few blocks simply establish any constants, or numerical parameters required to solve the device. The second block establishes any user defined inputs which also include device information for the device to be simulated. A note should be made that the code provided in Appendix III , and referenced in sections here, is heavily commented and the reader should reference the comments for any confusion that may arise. clear all; information close all; clc; TSTART = tic; format long eng; %Clears all previously stored variable %Closes all previously open plots %Clears the command window to blank %Starts a timer to record solution time %Formats the data to engineering mode This block is a simple setup that is used before any new simulation is run. All stored variable information in the Matlab workspace is cleared , all open plots are closed, and the command window is cleared. This section also starts a timer, which is used to investigate the speed of the computation algorithm for various device setups. 67 PAGE 77 Every calculated value that i s s tored in a variable i s acce ss ible after the program i s run in the main Matlab workspace . For example , the last iterated value for the error values is stored in an error vector. The program can be ea s ily modified to run through one iteration at a time allowing the u s er to inve stig ate the change in error over each iteration. Every opportunity for the u s er to acce ss the calculated data ha s been provided so that any insight into tho s e variables can be gained if nece s sary. One of the main rea s on s the author chose to u s e Matlab wa s due to the environments natural capability to work with matrice s . Everything in Matlab , even a one dim e n s ional variable or con s tant i s treated a s a matrix . One of the key feature s in SDS is the use of Cell arrays in Matlab. This allows the user to store various types of data in each cell. In fact , the frrst cell could have a t a ndard 3x 3 m a trix , while the s econd cell ha s a s tring text in it for example . The cell array s are treated like matrice s a s well and indexed s imilarly. All s tandard matrix operation s work , provided that matrix dimen s ion s a g ree , etc . Thi s work s elegantly with the tri-diagonal block matrices. E a ch matri x coefficient i s a 3 x 3 matrix but the calculated v alues are s tored in a c ell array where each element i s a 3x3 matrix. The cell array i s u s ed in conjunction with the re c ur s ive algorithm for the s olution proce ss . The delta _ y array , yo, andy array s are cell array s as well and the in s tantiation proce s i s s hown below . delta y = cell(1,L); delta-y(:) = {[0; 0 ;0]}; matrix yo = cell ( 1, L ) ; yo ( : ) = { [ 0 ; 0 ; 0] } ; matrix y = cell ( 1 , L ) ; y(: ) = { [ 0 ; 0 ;0]}; matrix Ap r c ell(l,L); Bp r cell(1,L); C p r cell( 1,L); Fpr cell (1,L); %Defines an empty cell from 1 to L %Intializes each cell to a 3x3 empty %Defines an empty cell from 1 to L %Initializes each cell to a 1x3 empty %Defines an empty cell from 1 to L %Initializes each cell to a lx3 empty delta_y{N } = inverse _func(Bpr{N})*Fpr{N} ... inverse_func(Bpr{N})*Cpr{N}* d elta_y{N+1}; 68 PAGE 78 This calculates the error during each iteration. Following the matrix dimensions gives delta_y{N} = 3x3 X 3xl-3x3 X 3x3 X 3xl = 3xl, which agrees. This method of using cell arrays is efficient and compact and is highly recommended. Even the update process is compact as shown below. In one line, the original trial potential, yo { N} is updated with the calculated error in delta_y { N} at each division point N. Remember the boundary conditions have no error associated with them as they are fixed. y{N} = yo{N} + delta_y{N}; After each iteration, the solution needs to be checked against the error tolerance. There is effectively two main ways to accomplish this. The first method requires that the absolute value of each calculated error value for the three fundamental variables, n, p, and sai, are checked against the convergence tolerance at each division point. If every single point has converged then the solution has been reached. This is somewhat inefficient, however is accurate for the whole device. A second method checks the normalized value of the entire error vector (normalized value of the error at each division point) against the convergence tolerance. If the normalized value is small enough the solution has converged. The process for one error vector is shown below. templ = abs(norm(p_error)/norm(po)); if (templ > etol) flagl 1; else flagl 0; end Once all three flags have been cleared to zero (they are initially set to one), then the device has converged. The solution algorithm loops until this occurs or the iteration limit set by the user has been reached. The default is 15 iterations. The final part of the "main" code or function simply plots the results, whether the device solution converged or not. The electric field, carrier densities, and device potential are all plotted against the spatial variable, x. 69 PAGE 79 The rest of the SDS program is a series of functions that are called by the "main" function or sub-functions. This is where the matrix coefficients are calculated and the partial derivatives as well. The code for one element of the B_func is shown below . It references several functions itself. b(1,1)= ( 1/hpr(N) )*((mu_func(M,h,po,no,Ntot,saio,1)/h(M)) ... *lambda_func(M,saio,1)-(mu_func(M-1,h,po,no,Ntot,saio,1)/h(M-1)) . . . *lambda_func(M-1,saio,-1))+ diffsrh_func(N,po,no,Ntot,1); Referencing back to Chapter 3 Section 3.5, you can see that there are several functions to support the calculation of the matrix coefficient. Mobility is calculated through a function so that models of temperature and effective doping can be used , and the lambda terms and recombination terms are supported by function calls as well. This allows for efficient storage of the calculated data. We don't need the lambda value stored across the whole device because unless you are troubleshooting , there is no interest in these values. The lambda function code is shown below. Recalling the various lambda parameters from Chapter 3, the beta term is also referenced as a function. This way beta(M) and beta(M-1) can be calculated in one function by passing in the current division point. The other important note for the lambda function is the use of the variable, I. This variable is passed in so that the proper sign is used for the exponential terms. In this manner, one function can support all of the lambda terms from Chapter 3 in a compact fashion. Vt = 0.02586; theta = 1/Vt; temp= beta_func(M,saio); if(I == 1) if(abs(temp) < 1E-5) lambda out (1/theta)*(1/(1-(1/2)*temp+(1/6)*tempA2)); else lambda out (1/theta)*(temp I (1-exp(-temp))); end elseif(I == -1) if(abs(temp) < 1E-5) lambda out= -(1/theta)*(1/(1+(1/2)*temp+(1/6)*tempA2)); else 70 PAGE 80 lambda out (1/theta}*(temp I (l-exp(temp))); end end end One final note on this function, which also applies to the difflmabda function as well, is that for special cases, the actual value of lambda needs to be approximated or some of the resulting matrix coefficient calculations will become undefined due to a division by zero. For the default device in this program, a step junction profile is used. This results in equal values of the device potential in each bulk region which will cause the output of the lambda function to go to zero or undefined. Therefore, depending upon the value of beta (which is the difference of the device potential between two division points) , instead of using the direction equation and approximated equation is used. In e sence a series expansion of the lambda terms are performed so that the exponentials are approximated. The technique used here was known as Bernoulli's function. It is also applied to the difflambda function which is used to calculate the partial derivative terms of the lambda coefficient. Several other functions exist but won ' t be discussed here. As mentioned , each block of code is commented in Appendix C. The remaining functions are similar in nature in that they support the overall solution process, many being to support the calculation of the matrix-vector coefficients , like generation and recombination . The current function will calculate the current density equations and so on. A basic flowchart of the program is provided in Figure 5 . 1 . This should be used in conjunction with Appendix C to better understand the flow of the program and where various calculation are made in reference to the discussions throughout the report. 71 PAGE 81 Set Constants / Numerical Parameters (T,VT , eo , q ,it,itl, L , etol) Set Device Parameters (Na,Nd , er , es , length, Vl , Vr) Define memory space (for faster computation) Set boundary conditions and trial potentials (n =nO, p=pO,sai = saio) Solve matrixvector equation u s ing the recursive algorithm Update trial potentials with calculated errors NO 0 It = itl + I Solution Converged ? YES Plot results (E-field, n, p , potential , etc.) YES Figure 5.1 -SDS BASIC FLOWCHART 72 PAGE 82 6. Summary Throughout this report an attempt has been made to thoroughly document the procedure used to develop a numerical solution to the simple PN junction diode. The basic device equations were transformed in to numerical equation that could then be solved in any desired manner. This program implemented the implicit Newton method as the chosen solution procedure to the system of numerical equations . A recursive algorithm was used to solve the system of numerical equations which had been placed in matrix form resulting in a tri-diagonal block matrix . The program is intended to allow for different PN junction diode characteristics to be studied by giving the user the capability to alter the default device already built into the program. Parameters like doping profile, device length , and device materials can all be changed to suit the simulation desired. Numerical parameters like convergence tolerance, mesh points, and mesh spacing can also be modified to suit the users need or even to help solidify the understanding of the numerical parameter's role in defining the obtained solution. The solution to the default PN junction was used to compare with the theoretical solutions obtained from making certain assumptions, like the depletion approximation, when solving the standard differential equations. The work done to develop this program lays the groundwork necessary to develop a more in-depth simulator capable of solving multiple semiconductor devices over a wide range of applications. The base program can easily be modified and updated to include further enhancements and capabilities without drastically changing the nature of the program. The program was written modularly and only portions need be modified at a time . Applications like zener breakdown, reverser recovery time, and more can be studied with the use of this simulator, with some suitable setup to the program. 73 PAGE 83 APPENDIX A Integral Form of the Current Density Equation As was seen in Chapter 3, Section 3.5, the current density equations were not tran formed directly from the difference method as was the continuity and Poisson equations . The main reason behind this is to generate stability within the tri-diagonal block matrix to ensure proper convergence when using the implicit method. It is possible however to employ the standard difference technique on the current density equations and substitute them directly into the continuity equations, but this can lead to numerical instability if the potential changes drastically between mesh points . To examine why this occurs let us re-define the current density equation for holes, repeated from Chapters 2 and 3 in equation (A. I) and then transform the equation into standard difference form as shown in equation (A.2a and A.2b). [4] dp df/1 J =-qD --qp p" "dx " dx (A.1) J(M)= qf.L"(M)[p(N+l)(N)] " Bh(M) p ( qf.L" (M) . [lf/(N + 1) -lj/(N)]J. ( p(N) + p(N + 1)) h(M) 2 (A.2a) qf.L"(M -1) [ ] J (M-1)=p(N)-p(N-1) " Bh(M -1) -(qf.L" (M -1). [ (N) _ (N -l)]J. ( p(N -1) + p(N)) h(M -1) If/ If/ 2 (A . 2b) The continuity equation for holes is again given in equation (A.3) in difference form. Equations (A.2a and A.2b) are then substituted into equation (A.3) and rearranged with like terms and simplified to obtain equation (A.4) . 74 PAGE 84 1 J (M)-J (M -1) P . '' -G(N)+U(N)=O q h (N) (A.3) Expanding J P (M) and J P (M -1) J (M)=qf.lp(M) (N+1)+ qf.l p(M) (N)qf.lp(M) (N+1) (N) P Bh(M) p Bh(M) p 2h(M) f// p qf.lp(M) (N+l) (N+l)+ qJlp(M) (N) (N) 2h(M) f// p 2h(M) f// p + qf.lp (M) (N) (N + 1) 2h(M) f// p and .. . J (M -1) =qf.l/M -1 ) (N) + qf.lp (M 1 ) (N -1)qj.lP (M 1 ) (N) (N -1) p Bh(M -1) p Bh(M -I) p 2h(M -1) f// p _ qf.l p (M 1 ) (N) (N) + qf.lp (M 1 ) (N -1) (N -1) 2h(M -1) f// p 2h(M -1) f// p + qf.lp (M 1 ) (N -1) (N) 2h(M -1) f// p Subtracting both current densities . . . _!_. J (M)-_!_ 1 (M-1)= f.lp(M) (N+l)+f.lp(M) (N) q p q p Bh(M)P Bh(M)P f.l (M) f.l (M) f.l (M) P f//(N + 1) p(N)P f//(N + 1) p(N + 1) + P f//(N) p(N) ... 2h(M) 2h(M) 2h(M) f.l " (M) (N) (N 1) f.l" (M 1 ) (N)f.l " (M -1) (N -1) + 2h(M) f// P + + Bh(M -1) P Bh(M 1) P ... f.l p (M 1 ) (N) (N -1) f.l " (M -1 ) (N) (N) + 2h(M -1) f// p + 2h(M -l) f// p ... f.lp (M -1) f.l p (M -1) f//(N-1)p(N-1)-f//(N-1)p(N) 2h(M -1) 2h(M -1) 75 PAGE 85 Combining like terms ... [ f.lp (M -1) f.lp (M -1) f.Lp (M -1) l . Bh(M -1) -2h(M -1) lf/(N -1) + 2h(M -1) lf/(N) . p(N -1 ) f.lp (M -1) f.l p (M -1) f.Lp (M -1) f.lp (M) Bh(M -1) + 2h(M -1) lf/(N)2h(M -1) lf/(N -1 ) + Bh(M) + p(N) f.lp (M) (N 1) f.lp (M) (N) -2h(M) If/ + + 2h(M) If/ +[-f.Lp(M)_f.LP(M) (N+1)+f.L"(M) (N)] (N+l) Bh(M) 2h(M) If/ 2h(M) If/ p Finally , simplifying ... f.Lp( M -1) [2 ] ___:__ __ -+ lf/(N -1) -ljf(N) p(N -I)+ 2h(M -1) e f.L"(M -1) [2 ] ___:__ __ -ljf(N -1) + lf/(N) ... 2h(M -I) e . p(N) f.1 (M) [2 ] + P â€¢ -+lj/(N)-Ij/(N +1) 2h(M) e (A.4) f.Lp(M) [ . 2] I + lf/(N) -ljf(N + 1)-p(N + 1) 2h(M) e =h. (N)[ G(N)-U(N)] If we assume that the mobility and mesh spacing is con tant and that the potential variation between the three points , N-1, N , N+l is linear , that is lj/(N -I)+ lj/(N) = lj/(N) -ljf(N + 1) = V0 , then we obtain equation (A.5) 76 PAGE 86 [(! + V0 } p(N p(N)+( V0 !} p(N +I)] (A.5) = h ' (N)[ G(N)U(N)] Examining equation (A.5) shows that unless the generation or recombination terms grow to proportions that deform the matrix coefficients on the left, the diagonal dominance of the equation is only lost if IV0 I which would allow the suband B super-diagonal elements to dominate the equation . The middle element p(N) is not affected by the change in potential from mesh points and so the only way to combat this issue when using the standard difference form is to decrease the mesh spacing to sufficient size to lower the potential change between mesh points. However, this can lead to an extensive amount of device grid points leading to excessive computation times. [1] This issue can be re olved by tran forming the current density equation into integral form and then applying the difference technique. This allows for deviation from diagonal dominance to be negligible to the overall solution algorithm , according to Sharfetter and Gummel. [8] To do this , we assume that in the range of x , between two division points Nand N+ 1 , the electric field, mobility , and current density are constant. Then equation (A.l) can be rearranged to integrate the hole density with respect to x to obtain equation (A.6). fap =BE r pax-+Jp [ax , where q , , kT D P = -J.Lp, q E=_alf/ B=.!L ax' kT qJ.LP(M)E(M) eBE< M l h < M > qJ.LP (M)E(M) lp(M)=-1 8E( M ) h ( M J p(N)+ 1 8E( M )h(M) p(N+1) (A.6) -e e 77 PAGE 87 To arrive at a more simplified solution, as shown in Chapter 3, we must do some more algebraic manipulation of equation (A.6). First , multiply the first term in equation (A.6) by -e8f: ( M ) h ( M ) --8f: ( M ) h ( M ) . . -e , which gives, qj1 (M)E(M). efE( M ) h ( M l . e 8f: < M l h < M > qfl (M)E(M) J p(M) = P (1 8f: ( M ) h ( M ) ) ( 8f: ( M ) h ( M ) ) p(N) + 1 P 8f:(M) h ( M ) p(N + 1) -e -e -e Combing the exponential terms gives, J (M)= qfl,,(M)E(M) p(N)+ qf1r(M)E(M)p(N+1) then, p ( 1 _ e -8( M ) h ( M ) ) }e8(M)h( M ) J p(M) = _q_[ Ar1 (M)p(N) + A r2(M)p(N + l)] h(M) 1 (M)= (M)If/(N)-1f(N+1) /Lpl f.lp J -/}( M ) -e A (M)= (M)If/(N)-!f(N+l) p2 Jip 1 /} ( M ) -e {J(M) = B[ lf/(N) -!f(N + 1)] where, (A.7) Equation (A.7) is exactly as shown in Chapter 3 and the same process can be used to arrive at a stable equation for the electron current density. This process is omitted here but is also the same and for both J P (M -1) and 111 (M -1) as well. In Chapter 5 , we will see that only one lambda equation function must be coded, but depending upon the variables passed into the function we can calculate the desired lambda terms . The only difference between Ar1 (M) and A r2 (M) is the sign of the beta term. A111 (M) and A112 (M) beta terms are simply swapped. So with one equation we can evaluate any of the four lambda terms. Al s o when calculations for M-1 division point are required, we simply calculate the function at the current division point less 1. 78 PAGE 88 Jp(M) =-q-[ A p1(M)p(N)+Ap2(M)p(N +1)] and h(M) J P (M -1) = q [ A P1 (M -1) p(N -1) +)., 2 (M -1) p(N)J h(M -1) p where, )., (M) = (M) 1/f(N) -llf(N + 1) p i f.l p 1 -/J(M) ' -e )., (M) = (M) 1/f(N) -llf(N + 1 ) and p 2 f.l p 1 /J( M ) ' -e )., (M -1) = II (M -1) 1/f(N -1 ) -1/f(N) p i r p 1 _ e -fJ( M I l ' A (M -1) = 11 (M -1) 1/f(N -1 ) -llf(N) and p2 r p 1 _ efJ< M 1 ) ' /J(M) = B[ 1/f(N) -llf(N + 1)) {J(M -1) = B[llf(N -1) -1/f(N)) J"(M)= and J" (M -1) = q [ A111 (M -l)n(N -1) + A112 (M -l)n(N)) h(M -1) where, )., (M)= (M)I/f(N)-11f(N+1) 111 f.l" 1 /J(M ) ' -e (M) = (M) 1/f(N) -llf(N + 1 ) and /L112 f.l" 1 -/J(M ) ' -e )., (M -1) = II (M -1) 1/f(N -1) -llf(N) nl t-"11 1-efJ( M 1 ) ' A (M -1) = 11 (M -1) 1/f(N -1 ) -llf(N) and 112 t-"11 1 /J ( M 1 ) ' -e /J(M) = B[llf(N) -llf(N + 1)) /J(M -1) = B[ 1/f(N -1) -llf(N)) Figure A.l-INTEGRAL FORM OF CURRENT DENSITY EQUATIONS 79 PAGE 89 Finally, if we take the new integral forms of the current densities in equation (A.7) just derived and plug them back into the continuity equation (A.3) and again assume that the mobility and mesh spacing are constant and that the potential variation between the three points, N-1, N, N+ 1 is linear we obtain equation (A.8). 1 -1--e_-P . p(N -1) + ( 1 +-. _1_) p(N) 1-eP 1-eP 1 + 1-eP . p(N + 1) = h[G(N)-U(N)] With this form, even if IV0 I 0 ..!_ 0 V0 , deviation from the diagonal dominant e (A.8) condition is nonessential. Using the integral form for the current densities greatly improves the matrix property, in that the nonsingularity condition is at least approximately satisfied even for high potential difference between two neighboring division points. 80 PAGE 90 APPENDIXB Derivation Techniques for the Partial Derivative Coefficients The following appendix will discuss the technique s used to derive the partial derivative s of the matrix-vector coefficients, A, B , and C. The derivation requires use of the product and chain rule from ba sic calculus . One example derivation is given for the hole current density and one example for the recombination terms. The rest of the derivations are left to the reader , however the partial derivatives for the coefficients A, B , and C are tabulated below so that they can be compared to the relationships given in Chapter 3 Section 3.5. Doing so will help develop the understanding required to program one function that supports the various differential lambda term that greatly reduces the number of computations. Repeating the current density equation from Appendix A , J p(M) = _q_[A. I (M)p(N) + A.p2(M)p(N + 1)] h(M) '' ') (M)= (M) 1/f(N)-'If(N+l) /l..p l Jlp l -{J(M) -e A. (M)=Jl (M)I/f(N)-1/f(N+l) p2 P 1 _ efl PAGE 91 To take the partial derivative of the lambda term with respect to the potential at mesh point N we can start by first applying the product rule since the lambda terms is not a straightforward function of the potential. f ( f//(N)) = f//(N)-f//(N + 1) g ( f//(N)) = ( 1-e P < M l f Using the product rule: where, However , before going further the partial derivative of the newly defined g term requires use of the chain rule, so that is evaluated fust. to take d [ g ( f//(N)) J, let g ( f//(N)) be a function of u then , a f//(N) g ( f//(N)) = u -1 and dg ( f//(N)) = -u-2 where, du u = 1-e P < M l and du = BeP PAGE 92 The same procedure is followed for the lambda term 2. aA-p2(M) a [ J a 'lf(N) = Jip (M) a 'lf(N) f ( 1/f(N)) g ( 1/f(N)) where, f ( 1/f(N)) = 'lf(N)-'lf(N + 1) g ( 1/f(N)) = ( 1e P C M > r Using the product rule: = .U, (M) { f ( \fi(N)) [ g ( \fi(N)) J + g ( \fi(N)) a:rN) [ f ( \fi(N)) J} However, before going further the partial derivative of the newly defined g term requires use of the chain rule, so that is evaluated first. to take a [ g ( 1/f(N)) J, let g ( 1/f(N)) be a function of u then, a 'lf(N) g ( 1/f(N)) = u -1 and ag ( 1/f(N)) = -u-2 where, au u = 1e P < M > and au = -eeP < M > a'lf(N) Using the chain rule: ag('lf(N)) = ag('lf(N)). au =-(l-ePCM ) ( (-BeP < M > ) a 'lf(N) au a 'lf(N) finall y , aA.p2 (M) = J1 (M){( 1/f(N)-'lf(N + 1) )( 1ePCM l )-2 ( ee P < M l ) + ( 1eP< M l ) -1} a 'lf(N) p 83 PAGE 93 Essentially the same procedure is used for the recombination and generation partial derivatives. An example of the recombination derivation is given below. U(N) = p(N)n(N)-ni2 r" ( p(N) + n i ) + r , ( n(N) + n i ) au (N) = a [( p(N)n(N)n ; ) ( r" ( p(N) + n i ) + r , ( n(N) + n i ) f J ap(N) ap(N) _aU---'-(N---'-) = a [( p(N)n(N)-n 2 ) ( r , ( p(N) + n i ) + r , ( n(N) + nJ f J an(N) an(N) ' ' Using the product and chain rule: a [(p(N)n(N)n ; ) ( r " (p(N) + ni) + r, J (n(N) + ni )f J ap(N) = (p(N)n(N)-n i2 ) ( -r" ( r " (p(N) + n i ) + r , (n(N) + n i )f) ... + n(N)( r " (p(N) + n i ) + r , (n(N) + ni) f so, aU(N) n(N)-r"U(N) ap(N) r"(p(N)+ni)+r,(n(N)+ni ) 84 PAGE 94 APPENDIXC Semiconductor Device Simulator Implemented %-------------------------------------------------------------------%Title: 1D PN Step Junction Numerical Device Simulator %Description: Using pre-defined material parameters for Silicon, a %device can be described numerically along a grid at which the %physical differential equations are solved for the three basic %device parameters, potential, and carrier concentrations, herein %referred to as sai, n, and p. %Written By: Shawn Pace %Date: 2011 UCD Thesis Project %Contributions: Many thanks to Dr. Hamid Fardi whose knowledge and %previous solution program fardev, without which, would have made %this nearly impossible. %-------------------------------------------------------------------clear all; information close all; clc; TSTART = tic; format long eng; %Input Parameters Nd = 1E15; Na = 1E15; T = 300; ni nieff func( T ) ; Vt 0 .02586; eo 8 .854E-14; er 11.7; es eo * er; length = 10E-6; q = 1.602E-19; theta = 1/Vt; %Clears all previously stored variable %Closes all previously open plots %Clears the command window to blank %Starts a timer to record solution time %Formats the data to engineering mode %Donor Concentration %Acceptor Concentration %Temperature in Kelvin %Intrinsic Carrier Concentration %Thermal Voltage %Relative Permittivity Constant %Relative Permittivity of Silicon %Permittivity of Silicon %Device C %Electron Charge %Boltzmann Factor 85 PAGE 95 alpha_pO = 2.25E7; %Generation defined constants; see Chapter E_pO = 3 .26E6; %Generation defined constants; see Chapter alpha_no = 3.80E6; %Generation defined constants; see Chapter E nO = 1.75E6; %Generation defined constants; see Chapter m=1; %Generation defined constants; see Chapter %Initial applied voltage: Vl 0; %Applied voltage to left contact Vr = 0; %Applied voltage to right contact %Numerical Parameters: it = 0; %Iteration number initialized to zero 5 5 5 5 5 itl = 15; %Max iteration limit to find solution within %a given tolerance L = 1001; etol = 1E-6; %Numerical Device Solution Points %Convergence Tolerance; calculated and %normalized error must be less than etol to %converge %Distance between each device point in em %10A2 is a conversion factor for um>cm delta (length/ (L-1)) *10A2; %Define uniform mesh points: x(1:L)=O; %Memory space defined for device points for(N=2:1:L) x(N) = x(N-1)+delta; end %Mesh points are defined from 1 to L in %increments of delta %Memory space defined for mesh spacing; h is uniform mesh spacing %and is equal to delta for a uniformly spaced device. hpr (h-prime) %is the auxiliary defined mesh spacing but is also equal to h for a %uniform device h(1:L-1)=0; hpr(1:L-1)=0; for(M=1:1:L-1) h(M) = x(M+1)-x(M); end for(N=2:1:L-1) M=N; hpr(N) = (1/2)*(h(M-1) + h(M)); end %Define Memory Space: x1(1:L)=O; gamma(1:L)=O; %Not currently used %Device doping profile 86 PAGE 96 Ntot(1:L)=O; saio(1:L)=0; sai(1:L)=0; po(1:L)=O; no(1:L)=O; E(1:L)=0; %Device doping total concentration (Nd + Na) %Device potential previous solution delta y = cell(1,L); delta=y(: ) = { [0;0;0]}; yo = cell ( 1, L); yo ( : ) = { [ 0 ; 0 ; 0 l } ; y = cell (1,L); y(:) = {[0;0;0]}; %Device potential final solution %Device hole concentration %Device electron concentration %Device energy field %Defines an empty cell from 1 to %Initializes each cell to a 3x3 %Defines an empty cell from 1 to %Initializes each cell to a 1x3 %Defines an empty cell from 1 to %Initializes each cell to a 1x3 L empty matrix L empty matrix L empty matrix %y and yo above are used to store the updated values after adding %the error for each iteration. They store sai, n, and p. %Boundary Conditions: %The solution of the device equations require a two-point boundary %condition and below are the boundary conditions for the device %length, doping profile, carrier concentrations, and built-in %potential x(1) = 0; x(L) = length*10A2; gamma(1)= -Na; gamma(L)= Nd; Ntot(1)= Na; Ntot(L)= Nd; po(1)=(-gamma(1)/2)*(1+sqrt(1+((2*ni)/gamma(1))A2)); no(1)=niA2/po(1); no(L)=(gamma(L)/2)*(1+sqrt(1+((2*ni)/gamma(L))A2)); po(L}=niA2/no(L); phi_p = -Vt*log(po(1)/(ni)); saio(L) = Vr + Vt*log(no(L)/(ni)) -phi_p; saio(1) = Vl; yo { 1} [po ( 1) ; no ( 1) ; saio ( 1) J ; yo{L} = [po(L) ;no(L) ;saio(L)]; %Numerical Solution %Trial Values: %Before the device equations can be calculated, trial values for the %carrier concentrations and device potential must be provided across %the whole device. for(N=2:1:L-1) %#ok<*N04LP> if(x(N) <= (length*10A2)/2) gamma(N) = -Na; Ntot(N)= Na; else gamma(N) = Nd; Ntot (N) = Nd; 87 PAGE 97 end if(sign(gamma(N)) == -1) po(N)=-gamma(N); no(N)=-(niA2/gamma(N)); else po(N)=niA2/gamma(N); no (N) =gamma (N) ; end saio(N) = Vt*log(abs(gamma(N))/ni)*sign(gamma(N)) -phi_p; if(saio(N) < 1E-9) saio(N) = 0; end yo{N} = [po(N) ;no(N) ;saio(N)]; end %Figure 1 plots the trial potential calculated above. This should be %a step profile for a step junction. figure(1) plot(x.*10A-2,saio) xlim ( [0 length]) ylim([O. O 1.0]) title('saio-initial trial potential', 'FontSize' ,16) %This pre-defined memory space is to store the error calculated %after each iteration for device potential and carrier %concentrations. This was mainly used for troubleshooting purposes %during development of the program, however these values can be %plotted to give an understanding of how the error changes during %each iteration. p_error(1:L)=0; n_error(1:L) =O; v _error(1:L) =O; %This pre-defined memory space is defining an empty cell from 1 to %L; it is not necessary to intialize each cell to 3x3 empty matrices %as above, because each cell can hold whatever the programmer %chooses. These cells hold the 3x3 matrices used to calculate the %error for each device point in the recursive equation that is %employed. Apr cell(1,L); Bpr cell(1,L) ; Cpr cell(1,L) ; Fpr cell(1,L) ; %Since the boundary conditions are pre-defined, there is no error %associated with grid points 1 and L, so only 2 through L-1 must be %solved; please refer to the appendices in the Kurata book N=2; Bpr{2} Cpr{2} b_func(N,h,hpr,po,no,Ntot,saio); c func(N,h,hpr,po,no,Ntot,saio); 88 PAGE 98 Fpr{2} = f_func(N,h,hpr,po,no,Ntot,gamma,saio); for(N=3:1:L-1) end Apr{N} = a_func(N,h,hpr,po,no,Ntot,saio); Bpr{N} = b_func(N,h,hpr,po,no,Ntot,saio) Apr{N}* ... inverse_func(Bpr{N-1})*Cpr{N-1}; Cpr{N} = c_func(N,h,hpr,po,no,Ntot,saio); Fpr{N} = f_func(N,h,hpr,po,no,Ntot,gamma,saio) Apr{N}* ... inverse_func(Bpr{N-1})*Fpr{N-1}; %Once Apr (A prime), Bpr, Cpr, and Fpr are calculated across the %whole device, the recursive equation can be used to calculate the %error at each device point. This error is then added to the trial %potential and if the convergence criteria are satisfied, the %solution has been found. for(N=L-1:-1:2) delta_y{N} = inverse_func(Bpr{N})*Fpr{N} ... inverse_func(Bpr{N})*Cpr{N}*delta_y{N+1}; end %Error for each variable is extracted for review if necessary for(N=1:1:L) end p_error(N) n_error(N) v_error(N) delta_y{N} (1); delta_y{N} (2); delta_y{N} (3); %This section adds the calculated error to the trial values; recall %that y is a cell matrix with each cell containing a 1x3 matrix and %since delta_y is of the same type they can easily be added %together. Then each updated value is extracted into its own vector %for convergence determination and also replaced in the trial value %cell matrix yo. for(N=1:1:L) y{N} = yo{N} + delta_y{N}; po(N) = y{N} (1); no(N) = y{N}(2); saio(N) = y{N}(3); if (saio (N) < 0) saio (N) = 0; end yo{N} = [po(N) ;no(N) ;saio(N)]; end 89 PAGE 99 %Variable temp1, 2, and 3 contain the normalized calculation for the %updated (trial + error) value which is then compared to the %convergence tolerance, etol and a flag is set to signify %convergence or not. temp1 = abs(norm(p_error)/norm(po)); if (temp1 > etol) flag1 1; else flag1 0; end temp2 = abs(norm(n_error)/norm(no)); if (temp2 > etol) flag2 1; else flag2 0; end temp3 = abs(norm(v_error)/norm(saio)); if (temp3 > etol) flag3 1; else flag3 0; end %Whether the solution was found, converged, or not, the iteration %number is still increased by 1 to signify one more iteration was %completed. it = it + 1; %This while loop repeats the solution algorithm above until the %iteration limit is reached or the solution found. while(((flag1 II flag2 II flag3) == 1) && (it< itl)) N=2; Bpr{2} Cpr{2} Fpr{2} b_func(N,h,hpr,po,no,Ntot,saio); c_func(N,h,hpr,po,no,Ntot,saio); f func(N,h,hpr,po,no,Ntot,gamma,saio); for(N=3:1:L-1) Bpr{N} = b_func(N,h,hpr,po,no,Ntot,saio) .. . -a func(N,h,hpr,po,no,Ntot,saio)* .. . inverse func(Bpr{N-1}l*Cpr{N-1}; Cpr{N} = c_func(N,h,hpr,po,no,Ntot,saio); 90 PAGE 100 Fpr{N} = f_func(N,h,hpr,po,no,Ntot,gamma,saio) ... -a func(N,h,hpr,po,no,Ntot,saio)* ... end for(N=L-1:-1:2) delta_y{N} = inverse_func(Bpr{N})*Fpr{N} ... -inverse_func(Bpr{N})*Cpr{N}*delta_y{N+1}; end for(N=1:1:L) p_error(N) n_error(N) v_error(N) delta_y{N} (1); delta_y{N} (2); delta_y{N} (3) ; end for(N=1:1:L) y{N} = yo{N} + delta_y{N}; po(N) = y{N}(1); no(N) = y{N} (2); saio(N) = y{N}(3); if(saio(N) < 0) saio (N) = 0; end yo{N} = [po(N) ;no(N) ;saio(N)]; end temp1 = abs(norm(p_error)/norm(po)); if (temp1 > etol) flag1 1; else flag1 0; end temp2 = abs(norm(n_error)/norm(no)); if (temp2 > etol) flag2 1; else flag2 0; end temp3 = abs(norm(v_error)/norm(saio)); if (temp3 > etol) flag3 1; 91 PAGE 101 else flag3 0; end it it + 1; end %This loop calculates the energy field at each device point once the %final device potential is found. for(M=1:1:L-1) E(M) = energy_func(M,h,saio); end %Figures 2-4 plot the zero-applied bias device potential and carrier %concentrations for a step junction PN device. figure(2) plot(x.*10A-2,po,x.*10A-2,no) xlim ( [ 0 length] ) % ylim([O 2E17]) title ('carrier densisties equilibrium', 'FontSize', 16) figure(3) plot(x.*10A-2,saio) xlim( [0 length]) title('potential-equilibrium', 'FontSize' ,16) ylim( [0 1. 0]) figure(4) plot(x.*10A-2,E) xl im ( [ 0 length] ) title('Energy-euilibrium', 'FontSize' ,16) time = toc(TSTART); time TSTART = tic; time %Stops timer for the equilibrium solution %Starts a timer for applied bias solution %Applied external Vl 0 .5; voltage (forward bias) Vr 0 .0; V1 0.0; %Applies a new external bias of 0.5 volts %Right contact is still grounded %Previous applied external voltage 92 PAGE 102 %Re-define Trial Values for applied external voltage; this equation %is a simple proportional relation to move the new "trial" values %with respect to the old applied voltage. The equilibrium solution %is used as the trial values for the applied voltage solution. for(N=2:l:L) %#ok<*N04LP> saio(N) = (1-((Vl-Vl)/(saio(L)-saio(l))))*saio(N) ... +((Vl-Vl)/(saio(L)-saio(l)))*saio(L); yo{N} = [po(N) ;no(N) ;saio(N)]; end % Re-define Boundary Conditions for applied external voltage; The %doping profile hasn't changed, and so therefore only the device %potential boundary conditions need updating. saio(l) = Vl; saio(L) = Vr + Vt*log(no(L)/(ni)) -phi_p; yo{l} [po(l) ;no(l) ;saio(l)]; yo{L} = [po (L) ;no (L); saio (L)]; g. 0 % % % % % % % This was used file for(N=l:1:L) out (N, 1) out (N, 2) out(N,3) out(N,4) end during development to extract the values to a text x(N); po(N); no(N); saio (N) ; %dlmwrite('output.txt' ,out) ; %This plots the updated device trial potential; it was mainly used %during development to ensure the starting trial potential was %correct. figure(S) plot(x.*l0A-2,saio) xlim ( [0 length]) title('saio-applied voltage trial 'Font8ize',l6) % ylim( [-1.2 1.2]) %Reset iteration number it = 0; itl = 15; %Iteration number re-initialized to zero %Max iteration limit to find solution within given %tolerance; this can be changed if solution is not %converging in the allotted iterations 93 PAGE 103 %The remaining portion of this code, down through the while loop is %exactly as above, and is performing the recursive algorithm on the %tri-diagonal matrix to calculate the error at each grid point and %continuing until the solution has converged or limit is reached. N=2; Bpr{2} Cpr{2} Fpr{2} b_func(N,h,hpr,po,no,Ntot,saio); c_func(N,h,hpr,po,no,Ntot,saio); f_func(N,h,hpr,po,no,Ntot,gamma,saio); for(N=3:1:L-1) end Apr{N} = a_func(N,h,hpr,po,no,Ntot,saio); Bpr{N} = b func(N,h,hpr,po,no,Ntot,saio) Apr{N}* ... Cpr{N} = c func(N,h,hpr,po,no,Ntot,saio); Fpr{N} = f=func(N,h,hpr,po,no,Ntot,gamma,saio) Apr{N}* ... inverse_func(Bpr{N-1})*Fpr{N-1}; for(N=L-1:-1:2) delta y{N} =inverse func(Bpr{N})*Fpr{N} ... --inverse_func(Bpr{N})*Cpr{N}*delta_y{N+1}; end for(N=1:1:L) p_error(N) n_error(N) v_error(N) end for(N=1:1:L) delta_y{N} (1); delta_y{N} (2); delta_y{N} ( 3) ; y{N} = yo{N} + 1 .0*delta_y{N}; po(N) = y{N} (1); no(N) = y{N} (2); saio(N) = y{N} (3); yo{N} = [po(N) ;no(N) ;saio(N)]; end temp1 = abs(norm(p_error)/norm(po)); if (temp1 > etol) flag1 1; else flag1 0; end temp2 abs(norm(n_error)/norm(no)); 94 PAGE 104 if (temp2 > etol) flag2 1; else flag2 0; end temp3 = abs(norm(v_error)/norm(saio)); if (temp3 > etol) flag3 1; else flag3 0; end it it + 1; while(((flag1 I I flag2 II flag3) == 1) && (it< itl)) N=2; Bpr{2} Cpr{2} Fpr{2} b func(N,h,hpr,po,no,Ntot,saio); c_func(N,h,hpr,po,no,Ntot,saio); f_func(N,h,hpr,po,no,Ntot,gamma,saio); for(N=3:1:L-1) end Bpr{N} = b_func(N,h,hpr,po,no,Ntot,saio) . . . -a func(N,h,hpr,po,no,Ntot,saio)* .. . Cpr{N} = c_func(N,h,hpr,po,no,Ntot,saio); Fpr{N} = f_func(N,h,hpr,po,no,Ntot,gamma,saio) ... -a func(N,h,hpr,po,no,Ntot,saio)* ... for(N=L-1:-1:2) delta_y{N} = inverse_func(Bpr{N})*Fpr{N} ... -inverse_func(Bpr{N})*Cpr{N}*delta_y{N+1}; end for(N=1:1:L) p_error(N) n_error(N) v_error(N) end for(N=1:1:L) delta y{N} ( 1) ; delta=y{N} (2); delta_y{N} (3); y{N} = yo{N} + 1.0*delta_y{N}; 95 PAGE 105 end end po(N) = y{N} (1); no(N) = y{N} (2); saio(N) = y{N} (3); yo{N} = [po(N) ;no(N) ;saio(N)]; temp1 = abs(norm(p_error)/norm(po)); if (temp1 > etol) flag1 1; else flag1 0; end temp2 = abs(norm(n_error)/norm(no)); if (temp2 > etol) flag2 1; else flag2 0; end temp3 = abs(norm(v_error)/norm(saio)); if (temp3 > etol) flag3 1; else flag3 0; end it it + 1; %Defined memory space for the SRH terms used to calculate current U(1:L)=0; %Energy field is re-calculated based on applied voltage for(M=1:1:L-1) E(M) = energy_func(M,h,saio); end %SRH Recombination is calculated across the device in order to %calculate the device current for the give applied voltage. for(N=2:1:L-1) U(N) = srh_func(N,po,no,Ntot); end 96 PAGE 106 %Based on the device equations, the current can be calculated by %numerically integrating the SRH Recombination values at each grid %point. A simple trapezoidal integration technique was used. This is %only useful for forward bias situations. Generation terms need to %be included for reverse bias current draw calculations. Jtot = q*((delta/2)*(U(1)+U(L)) + delta*sum(U(2:L-1))); %Finally figures 6 8 plots the device potential and carrier %concentrations for the given applied bias. figure(6) plot(x.*10A-2,saio) xlim ( [0 length]) ylim ( [0 1. 0]) title('potential-applied voltage','FontSize',16) figure(?) plot(x.*10A-2,E) xlim ( [0 length]) title('Energy applied voltage', 'FontSize' ,16) figure(B) plot(x.*10A-2,po,x.*10A-2,no) xlim ( [0 length]) % ylim( [0 2E17]) title('carrier dens1ties-applied voltage','FontSize',16) time1 = toc(TSTART); %Stops the clock for the applied bias %solution time 97 PAGE 107 function mu_out = mu_func(M,h,po,no,Ntot,saio,I) %This function returns the carrier concentration mobility based on %the index, I and grid point M, passed into the function. %Currently this function only returns the low E -field mobility as a %constant and for the hole mobility an I of 1 is used and for %electron mobility an I of -1 is used. The rest of the equations and %constants below are currently not used, but could be to develop the %mobility based on total concentration and E-field. % M -Device grid point of interest % h -Grid mesh spacing % po -Hole concentration % no -Electron concentration % Ntot -Total concentration % saio -Device potential % I -Index used to define which mobility, hole or electron, to % output % N=M; %Mobility Parameters Nrefp = 6 .3E16; % cmA-3 alphap = 0.76; %fitting paramter betap = 1; %fitting parameter mu_maxp = 495; % cmA2/V-sec mu_minp = 47.7; % cmA2/V-sec Ecp = 1.95E4; % V/cm Nrefn = 8.5E16; alphan = 0.72; betan = 2; mu_maxn = 1330; mu_minn = 65; Ecn = 8E3 ; % crn"'-3 %fitting paramter %fitting parameter % cmA2/V-sec % cmA2/V-sec % V/cm % fp % fn ((po(N)*no(N))A(1/2)/(2.04*Nrefp))Aalphap; ((po(N)*no(N) )A(1/2)/(2.04*Nrefn))Aalphan; fp fn 0 i 0 i % if(I == 1) % mu out =(((mu maxp-u minp)/((1+(Ntot(M)/Nrefp)Aalphap+fp))) ... % % elseif(I == -1} % mu out =(((mu maxn-u minn)/((1+(Ntot(M)/Nrefn)Aalphan+fn))) ... % % end 98 PAGE 108 if (I == 1) mu_out = 495; elseif(I == -1) mu out = 1330; end function nieff out = nieff func(T) %Returns the effective intrinsic carrier concentration based on the %provided device temperature. k = 8.62E-5; EG300 = 1.08; EG_alpha = 4.73E-4; EG_beta = 636; NC300 2.8E19; NV300 = 1.04E19; EG EG300 + EG alpha*((300A2/(300+EG beta)) -(TA2/(T+EG_beta))) i NC (T/300)A(372)*NC300; -NV (T/300)A(3/2)*NV300; assignin ( 1 base 1 1 1 EG' I EG) ; assignin ( 1 base 1 I 1 NC 1 INC) ; assignin ( 1 base 1 I 1 NV 1 1 NV); nieff out sqrt(NC*NV)*exp(-EG/(2*k*T)); 99 PAGE 109 function a = a_func(N,h,hpr,po,no,Ntot,saio) %This function calculates a 3x3 matrix based on the grid point %index, M, passed into the function. %This matrix, along with the b,c, and f_func makes up the %Tri-diagonal matrix for each device point a when combined can be %solved for the error using the recursive tri-diagonal solution %algorithm. % N -Device grid point of interest % h -Grid mesh spacing % hpr -Auxillary grid spacing; same as h for uniform device % po -Hole concentration % no -Electron concentration % Ntot -Total concentration % saio -Device potential %The main 3x3 matrices (a,b,c,f) are calculated upon the main grid points N, but the sub-equations are calculated on the auxiliary grid points. This is why M is set equal toN. M=N; a(1,1)= -(mu_func(M,h,po,no,Ntot,saio,1)/hpr(N)) ... *(lambda_func(M-1,saio,1)/h(M-1)); a(1,2)= 0; a(l,3)= -(1/(hpr(N)*h(M-1)))*(mu_func(M-1,h,po,no,Ntot,saio,l)* ... po(N-1)*difflambda func(M-1,saio,1)+ mu func(M--1,h,po,no,Ntot,saio,1)* ... po(N)*difflambda_func(M-l,saio,-1)); a(2,1)= 0; a (2, 2) = -(mu_func (M, h, po, no,Ntot, saio, -1) /hpr (N)) ... *(lambda_func(M-1,saio,-1)/h(M-1)); a(2,3)= -(1/(hpr(N)*h(M-1)))*(mu_func(M-1,h,po,no,Ntot,saio,-1) *no (N-1) ... *difflambda_func(M-1,saio,-1)+ mu_func(M-1,h,po,no,Ntot,saio,-1) ... *no(N)*difflambda_func(M-1,saio,1)); a(3,1)= 0 ; a(3,2)= 0 ; a(3,3)= 1/(h(M-1)*hpr(N)); end 100 PAGE 110 function b = b_func(N,h,hpr,po,no,Ntot,saio) %This function calculates a 3x3 matrix based on the grid point index, M, passed into the function. %This matrix, along with the a,c, and f func makes up the tri%diagonal matrix for each device point a when combined can be solved %for the error using the recursive tri-diagonal solution algorithm. % N Device grid point of interest % h Grid mesh spacing % hpr Auxiliary grid spacing; same as h for uniform device % po Hole concentration % no Electron concentration % Ntot Total concentration % saio -Device potential %The main 3x3 matrices (a,b,c,f) are calculated upon the main grid points N, but the sub-equations are calculated on the auxiliary grid points. This is why M is set equal to N. M=N; eo 8.854E-14; er 11.7; es eo * er; q = 1.602E-19; %Permittivity Constant %Relative Permittivity of Silicon %Permittivity of Silicon %Electron Charge b(1,1)= ( 1/hpr(N) )*((mu_func(M,h,po,no,Ntot,saio,1)/h(M)) ... *lambda_func(M,saio,1)-(mu_func(M-1,h,po,no,Ntot,saio,1)/h(M-1)) ... *lambda_func(M-1,saio,-1))+ diffsrh_func(N,po,no,Ntot,1); b(1,2)= diffsrh_func(N,po,no,Ntot,-1); b(1,3)= ( 1/hpr(N) )*((mu_func(M,h,po,no,Ntot,saio,1)/h(M)) ... *(po(N)*difflambda_func(M,saio,1)+po(N+1) ... *difflambda_func(M,saio,-1))+ ... (mu_func(M-1,h,po,no,Ntot,saio,1)/h(M-1))* (po(N-1) ... *difflambda_func(M-1,saio,1)+po(N)*difflambda_func(M-1,saio,-1))) i b(2,1)= -diffsrh_func(N,po,no,Ntot,1); b(2,2)= ( 1/hpr(N) )*((mu_func(M,h,po,no,Ntot,saio,-1)/h(M)) ... *lambda_func(M,saio,-1)-(mu_func(M-1,h,po,no,Ntot,saio,-1)/h(M-1)) ... *lambda_func(M-1,saio,1))-diffsrh_func(N,po,no,Ntot,-1); b(2,3)= ( 1/hpr(N) ) *((mu_func(M,h,po,no,Ntot,saio,-1)/h(M)) ... *(no(N)*difflambda_func(M,saio,-1)+no(N+1) ... *difflambda_func(M,saio,1))+ ... (mu_func(M-1,h,po,no,Ntot,saio,-1)/h(M-1))*(no(N-1) ... 101 PAGE 111 *difflambda_func(M-1,saio,-1)+no(N)*difflambda_func(M-1, saio, 1) ) ) ; b(3,1)= q/es; b(3,2)= -q/es; b(3,3)= -(1/hpr(N))*(1/h(M-1) + 1/h(M)); end 102 PAGE 112 function c = c_func(N,h,hpr,po,no,Ntot,saio) %This function calculates a 3x3 matrix based on the grid point %index, M, passed into the function. %This matrix, along with the a,b, and f func makes up the tri%diagonal matrix for each device point a when combined can be solved %for the error using the recursive tri-diagonal solution algorithm. % N -Device grid point of interest % h -Grid mesh spacing % hpr -Auxiliary grid spacing; same as h for uniform device % po Hole concentration % no -Electron concentration % Ntot -Total concentration % saio Device potential %The main 3x3 matrices (a,b,c,f) are calculated upon the main grid %points N, but the sub-equations are calculated on the auxiliary %grid points. This is why M is set equal toN. M=N; c(1,1)= (mu_func(M,h,po,no,Ntot,saio,1)/hpr(N)) ... *(lambda_func(M,saio,-1)/h(M)); c(1,2)= 0; c(1,3)= -(1/(hpr(N)*h(M)))*(mu_func(M,h,po,no,Ntot,saio,1) ... *po(N)*difflambda_func(M,saio,1)+mu_func(M,h,po,no,Ntot,saio,1) ... *po(N+1)*difflambda_func(M,saio,-1)); c(2,1)= 0; c(2,2)= (mu_func(M,h,po,no,Ntot,saio,-1i/hpr(N)) ... *(lambda_func(M,saio,1)/h(M)}; c(2,3)= -(1/(hpr(N)*h(M)))*(mu_func(M,h,po,no,Ntot,saio,-1) ... *no(N)*difflambda_func(M,saio,-1)+mu_func(M,h,po,no,Ntot,saio,-1) ... *no(N+1)*difflambda_func(M,saio,1)) c(3,1)= 0; c(3,2)= 0; c(3,3)= 1/(h(M)*hpr(N)); end 103 PAGE 113 function f = f_func(N,h,hpr,po,no,Ntot,gamma,saio) %This function calculates a 1x3 matrix based on the grid point %index, M, passed into the function. %This matrix, along with the a,b, and c func makes up the tri%diagonal matrix for each device point a when combined can be solved %for the error using the recursive tri-diagonal solution algorithm. % N -Device grid point of interest % h -Grid mesh spacing % hpr -Auxiliary grid spacing; same as h for uniform device % po -Hole concentration % no -Electron concentration % Ntot -Total concentration % gamma -Doping concentration % saio -Device potential %The main 3x3 matrices (a,b,c,f) are calculated upon the main grid %points N, but the sub-equations are calculated on the auxiliary %grid points. This is why M is set equal to N. M=N; eo er es q 8.8548-14; 11.7; %Permittivity Constant eo * er; 1.6028-19; %Relative Permittivity of Silicon %Permittivity of Silicon %Electron Charge 1/(h(M-1)*hpr(N)); gamma1 gamma2 gamma3 -(1/hpr(N)) * (1/h(M-1) + 1/h(M)); 1/(h(M)*hpr(N)); f(1,1)= -(1/(q*hpr(N)))*(current_func(M,h,po,no,Ntot,saio,1) -current func(M-1,h,po,no,Ntot,saio,1)) -srh_func(N,po,no,Ntot); f(2,1)= -(1/(q*hpr(N)))*(current_func(M,h,po,no,Ntot,saio,-1) ... -current func(M-1,h,po,no,Ntot,saio,-1)) + srh_func(N,po,no,Ntot); f(3,1)= (-q/es)*(gamma(N)+po(N)-no(N)) -gamma1*saio(N-1) ... -gamma2*saio(N) -gamma3*saio(N+1); assignin('base', 'gamma1' ,gamma1); assignin('base', 'gamma2' ,gamma2); assignin('base', 'gamma3' ,gamma3); end 104 PAGE 114 function current_out = current func(M,h,po,no,Ntot,saio,I) %This function calculates and returns the hole or electron current %density based on which index, I is passed into the function. %The current and hole density equations come from the numerical %equations established in order to solve the physical differential %equations. It uses the lambda and mobility functions as well. % M -Device grid point of interest % h Grid mesh spacing % po Hole concentration % no Electron concentration % Ntot Total concentration % saio Device potential % I -Index used to define which current hole or electron, % to output q 1.602E-19; %Electron Charge %The main 3x3 matrices (a,b,c,f) are calculated upon the main grid %points N, but the sub-equations are calculated on the auxiliary %grid points. This is why N is set equal to M. N=M; if(I == 1) current_out = (q/h(M))*(mu_func(M,h,po,no,Ntot,saio,1) ... *lambda_func(M,saio,1)*po(N)+ mu_func(M,h,po,no,Ntot,saio,1) ... *lambda func(M,saio,-1)*po(N+1)); elseif(I == -1) current_out = (q/h(M))*(mu_func(M,h,po,no,Ntot,saio,-1) ... *lambda_func(M,saio,-1)*no(N)+ mu_func(M,h,po,no,Ntot,saio,-1 ) ... *lambda_func(M,saio,1)*no(N+1)); end function beta_out = beta_func(M, saio) %This function is only a called function, usually by another %function, in order to calculate the beta parameter which is found %in some of the numerical equations, for example, a,b,c,d_func, %current func, etc. Physically beta has no meaning and is only a %variable. Vt = 0.02586; theta = 1/Vt; beta out= theta*(saio(M)-saio(M+1)); end 105 PAGE 115 function lambda_out = lambda_func(M,saio,I) %This function is only a called function, usually by another %function, in order to calculate the lambda parameter which is found %in some of the numerical equations, for example, a,b,c,d_func, %current_func, etc. Physically lambda has no meaning and is only a %variable. For calculation of the lambda parameter the Bernoulli %function was used. For very small numbers, the physical equation %breaks down and therefore must be expanded so that larger terms can %be neglected. This is basically one use of a series expansion. Vt = 0.02586; theta = 11Vt; temp= beta_func(M,saio); if(I == 1) if(abs(temp) < 1E-5) lambda out (1ltheta)*(11(1-(112)*temp+(116)*tempA2)); else lambda out (1ltheta)*(temp I (1-exp(-temp))); end elseif(I == -1} if(abs(temp) < 1E-5) lambda out -{1ltheta)*(11(1+(112)*temp+(116)*tempA2)); else lambda out (1ltheta)*(temp I (1-exp(temp))); end end function difflambda_out = difflambda_func(M,saio,I) %This function is only a called function, usually by another %function, in order to calculate the derivative of the lambda %parameter which is found in some of the numerical equations, for %example, a,b,c,d_func. temp= sign(I)*beta_func(M,saio); if(abs(temp) < 1E-5) difflambda out 0.5 + sign(I)*(templ6); else difflambda out (11(1-exp(-temp)}}*(1-templ(exp(temp)-1)); end 1 06 PAGE 116 function E_out = energy_func(M,h,saio) %This function calculates the energy field at a given device point %based on the given grid point, mesh spacing, and device potential. %The E-field is simply the difference between two grid point %potentials divided by the mesh spacing. N=M; E out ( saio (N+l)-saio (N) ) /h (M) ; end function invBpr_out = inverse func(Bpr) %This function is used to calculate the inverse of a 3x3 matrix; the %built in Matlab function det is used to calculate the determinate %of the matrix invBpr_out(l,l) invBpr_out(l,2) invBpr_out(l,3) invBpr_out(2,1) invBpr_out(2,2) invBpr_out(2,3) invBpr_out(3,1) invBpr_out(3,2) invBpr out(3,3) end (Bpr(2,2)*Bpr(3,3) Bpr(3,2)*Bpr(2,3))/det(Bpr); (Bpr(3,2)*Bpr(l,3) Bpr(l,2)*Bpr(3,3))/det(Bpr); (Bpr(l,2)*Bpr(2,3) Bpr(l,3)*Bpr(2,2))/det(Bpr); (Bpr(3,l)*Bpr(2,3)-Bpr(2,l)*Bpr(3,3))/det(Bpr); (Bpr(l,l)*Bpr(3,3) Bpr(3,l)*Bpr(l,3))/det(Bpr); (Bpr(2,l)*Bpr(l,3) Bpr(l,l)*Bpr(2,3))/det(Bpr); (Bpr(2,l)*Bpr(3,2) Bpr(3,l)*Bpr(2,2))/det(Bpr); (Bpr(3,l)*Bpr(l,2) Bpr(l,l)*Bpr(3,2))/det(Bpr); (Bpr(l,l)*Bpr(2,2) Bpr(2,l)*Bpr(l,2))/det(Bpr); 107 PAGE 117 function gen_out gen_func(N,h,hpr,po,no,Ntot,saio) M = N; q = 1 .602E-19; %Electron Charge alpha_pO = 2.25E7; E_pO = 3.26E6; alpha_no = 3 .80E6; E nO = 1.75E6; m=1; rna h(M)/(2*hpr(N)); mb h(M-1)/(2*hpr(N)); E = ma*energy_func(M-1,h,saio) + mb*energy_func(M,h,saio); alpha_p alpha_pO*exp(-(E_pO/abs(E))Am); alpha_n = alpha_nO*exp(-(E_nO/abs(E))Am); J_p = ma*current_func(M-1,h,po,no,Ntot,saio,1) + mb*current_func(M,h,po,no,Ntot,saio,1); J_n = ma*current_func(M-1,h,po,no,Ntot,saio,-1) + mb*current func(M,h,po,no,Ntot,saio,-1); gen_out = (1/q)*(alpha_p*abs(J_p) + alpha_n*abs(J_n)); end 1 08 PAGE 118 function gen_out = diffgen_func(N,h,hpr,po,no,Ntot,saio,I) q = 1.602E-19; %Electron Charge alpha_pO = 2 .25E7; E_pO = 3.26E6; alpha_no = 3.80E6; E no = 1 .75E6; m=1; M=N; rna h(M)/(2*hpr(N)); mb h(M-1)/(2*hpr(N)); E = ma*energy func(M-1,h,saio) + mb*energy func(M,h,saio); alpha_p alpha_pO*exp(-(E_pO/abs(E))Am);-alpha_n = alpha_nO*exp(-(E_nO/abs(E))Am); J_p = ma*current_func(M-1,h,po,no,Ntot,saio,1) + mb*current_func(M,h,po,no,Ntot,saio,1); J_n = ma*current_func(M-1,h,po,no,Ntot,saio,-1) + mb*current_func(M,h,po,no,Ntot,saio,-1); G11 = (ma/h(M-1))*alpha_p*(E_p0/EA2)*(abs(J_p)/q)*sign(E); G12 = ((ma*alpha_p)/h(M-1))*(mu_func(M-1,h,po,no,Ntot,saio,1)*po(N-1)*difflambda_func(M-1,saio,1) ... + mu_func(M-1,h,po,no,Ntot,saio,1)*po(N)*difflambda_func(M-1,saio,2))*sign(J_p); G13 = (ma/h(M-1))*alpha_n*(E_n0/EA2)*(abs(J_n)/q)*sign(E); G14 = ((ma*alpha_n)/h(M-1))*(mu_func(M-1,h,po,no,Ntot,saio,-1)*no(N-1)*difflambda_func(M-1,saio,3) ... + mu_func(M-1,h,po,no,Ntot,saio,-1)*no(N)*difflambda_func(M-1,saio,4))*sign(J_n); G31 = -(mb/h(M))*alpha_p*(E_p0/EA2)*(abs(J_p)/q)*sign(E); G32 = ((mb*alpha_p)/h(M))*(mu_func(M,h,po,no,Ntot,saio,1)*po(N)*difflambda func(M,saio,1) ... + mu_func(M,h,po,no,Ntot,saio,1)*po(N+1)*difflambda_func(M,saio,2))*si gn(J_p); G33 = -(mb/h(M))*alpha_n*(E_n0/EA2)*(abs(J_n)/q)*sign(E); G34 = ((mb*alpha_n)/h(M))*(mu_func(M,h,po,no,Ntot,saio,1)*no(N)*difflambda func(M,saio,3) ... 109 PAGE 119 + rnu_func(M,h,po,no,Ntot,saio,l)*no(N+l)*difflambda_func(M,saio,4))*si gn (J_n); G21 = -Gll G31; G22 = alpha_p*((mb/h(M))*(rnu_func(M,h,po,no,Ntot,saio,l)*po(N)*difflambda_ func(M,saio,l) . . . + rnu_func(M,h,po,no,Ntot,saio,l)*po(N+l)*difflambda_func(M,saio,2)) ... -(rna/h(M-l))*(rnu_func(M-l,h,po,no,Ntot,saio,l)*po(Nl)*difflambda_func(M-l,saio,l) . . . + rnu_func(M-l,h,po,no,Ntot,saio,l)*po(N)*difflambda_func(Ml,saio,2)))*sign(J_p); G23 = -Gl3 G33 ; G24 = alpha_n*((mb/h(M))*(rnu_func(M,h,po,no,Ntot,saio,l)*no(N)*difflambda_func(M,saio,3) ... + rnu_func(M, h ,po,no,Ntot,saio, l)*no(N+l)*difflambda_func(M,saio,4)) ... -(rna/h(M-l))*(rnu_func(M-l,h,po,no,Ntot,saio,-l)*no(Nl)*difflambda_func(M-l,saio,3) . . . + rnu_func(M-l,h,po,no,Ntot,saio,-l)*no(N)*difflambda_func(Ml,saio,4)))*sign(J_n); assignin ('base', 'rna', rna); assignin ('base', 'mb', mb) ; assignin( 'base', 'E' ,E); assignin('base', 'alpha_p',alpha_p); assignin('base', 'alpha_n',alpha_n); assignin ('base', â€¢ J_p â€¢ , J_p); assignin('base', 'J_n',J_n); assignin('base', 'Gll',Gll); assignin ('base' , 'Gl2' , Gl2) ; assignin('base', 'Gl3' ,Gl3); assignin ( 'base' , 'Gl4 ' , Gl4) ; assignin('base', 'G21' ,G21); assignin ( 'base', 'G22' , G22) ; assignin ('base', 'G23', G23) ; assignin ('base' , 'G24' , G24) ; ass ignin ( 'base ' , 'G31 ' , G31) ; assignin( 'base', 'G32' ,G32); assignin ('base', 'G33', G33); assignin ('base' , 'G34' , G34) ; if (I == 1) gen_out = (rna/h(M-l))*rnu_func(Ml,h,po,no,Ntot,saio,l)*lambda_func(M-l,saio,l)*alpha_p*sign(J_p); elseif(I == 2) gen_out = ((rna/h(M-l))*rnu_func(Ml,h,po,no,Ntot,saio,l)*lambda_func(M-l,saio,-1) ... 110 PAGE 120 + (mb/h(M))*mu_func(M,h,po,no,Ntot,saio,l)*lambda_func(M,saio,l))*alph a_p*sign (J_p) ; elseif(I == 3) gen_out = (mb/h(M))*mu_func(M,h,po,no,Ntot,saio,l)*lambda_func(M,saio,l)*alpha_p*sign(J_p); elseif(I == 4) gen_out = (ma/h(M-l))*mu_func(M-l,h,po,no,Ntot,saio,l)*lambda_func(M-l,saio,-l)*alpha_n*sign(J_n); elseif (I == 5) gen_out = ((ma/h(M-l))*mu_func(M-l,h,po,no,Ntot,saio,l)*lambda_func(M-l,saio,l) ... + (mb/h(M))*mu_func(M,h,po,no,Ntot,saio, l)*lambda_func(M,saio,-l))*alpha_n*sign(J_n); elseif(I == 6) gen_out = (mb/h(M))*mu_func(M,h,po,no,Ntot,saio, l)*lambda_func(M,saio,l)*alpha_n*sign(J_n); elseif(I == 7) gen_ out = Gll + G12 + G13 + G14; elseif(I --8) gen_ out = G21 + G22 + G23 + G24; elseif(I --9) gen_ out = G31 + G32 + G33 + G34; end end 111 PAGE 121 function srh_out = srh_func(N,po,no,Ntot) %This function calculates the Shockley Read-Hall Recombination term %at a given device grid point. Only the current carrier %concentrations are required to calculate this parameter. %The SRH term is used in order to calculate the current in a forward %biased device. Once a solution is found, the SRH terms can be re %calculated across the device and a simple numerical integration %method can be used to calculate the current. Employing these terms %in the overall tri-diagonal matrix gives a more accurate solution. ni = 1.45E10; %Intrinsic Carrier Concentration tau_pO tau no Nsrh_p Nsrh n tau_p tau n 0 .2E-7; 0.2E-7; 5E15; 5E15; tau_pO/(l+(Ntot(N)/Nsrh_p)); tau_nO/(l+(Ntot(N)/Nsrh_n)); srh_out = (po(N)*no(N) niA2)/(tau_n*(po(N)+ni) + tau_p*(no(N)+ni)); end 112 PAGE 122 function diffsrh_out = diffsrh_func(N,po,no,Ntot,I) %This function calculates the derivative of the SRH terms, as when %the numerical equations are expanded, it requires derivative of the %SRH terms as well. If I = 1, then the partial derivative of the %Recombination term is calculated with respect to the hole %concentration and if I = -1, then it is calculated with respect to %the electron concentration. ni = 1.45E10; tau_pO tau nO Nsrh_p Nsrh n 0 .2E-7; 0 .2E-7; 5E15; 5E15; %Intrinsic Carrier Concentration %Hole lifetime constant %Electron lifetime constant %Hole SRH Carrier concentration constant %Electron SRH Carrier concentration constant tau_p tau n tau_p0/(1+(Ntot(N)/Nsrh_p)); %Calculated hole lifetime tau_n0/(1+(Ntot(N)/Nsrh_n)); %Calculated electron lifetime if (I == 1) diffsrh_out = (no(N) -tau_n*srh_func(N,po,no,Ntot))/ ... (tau_n*(po(N)+ni)+ tau_p*(no(N)+ni)); elseif(I == -1) end diffsrh_out = (po(N) -tau_p*srh_func(N,po,no,Ntot))/ ... (tau_n*(po(N)+ni)+ tau_p*(no(N)+ni)); 113 PAGE 123 REFERENCES [1] R. Muller et al., Device El ec tronics for Int egrated Circuits, 3rd ed. New York , Wiley , 2003 , ch. 1-5. [2] G. Neudeck and R Pierret, The PN Jun ction Diode , 2nd ed., Reading , Addison Wesley , 1989 , ch. 1-3. [3] M. Pinto et al., " PISCES-II ," Stanford Electronics Lab., Dept. of Elec. Engr. , Stanford Univ. , California , Sept. 1984. [4] M. Kurata , Numerical Anal ysis for Semiconductor De v i ces, Lexington , D .C. Heath and Company , 1982 , ch . 1-4. [5] J . Stewart, Calculus Earl y Transcendentals, 41h ed. Pacific Grove, Brooks/Cole, 1999 , ch 11 & 14. [6] A. Kharab and R. Guenther , An Introdu ction to Numerical Methods , A MATLAB Approach, 2nd ed. Boca Raton , Chapman & Hall, 2006, ch. 9-12. [7] D. Winston , " www-ocs . colorado.edu/SimWindows/simwin.html" , 1999. [8] D.L. Scharfetter and H.K. Gummel, "Large-Signal Analysis of a Silicon Read Diode Oscillator ," IEEE Transaction s on Electron Devic es, no. ED-16 , January 1969 , pp , 64-77. 114 |