ECE539 - Advanced Theory of Semiconductors and Semiconductor Devices Numerical Methods and Simulation / Umberto Ravaioli Review of Conventional Semiconductor Device Models Based on Partial Differential Equations 1 Boltzmann equation The standard semi–classical transport theory is based on the Boltzmann equation X ∂f eE + v · ∇r f + · ∇k f = ∂t h̄ k0 [S(k0 , k)f (r, k0 , t)[1 − f (r, k, t)] −[S(k, k0 )f (r, k, t)[1 − f (r, k0 , t)] (1) where r is the position, k is the momentum, f (k, t) is the distribution function (for instance, a Fermi–Dirac in equilibrium), v is the group velocity, E is the electric field, S(k, k0 ) is the transition probability between the momentum states k and k0 , and [1 − f (k0 , t)] is the probability of non– occupation for a momentum state k0 . The summation on the right hand side is the collison term, which accounts for all the scattering events. The terms on the left hand side indicate, respectively, the dependence of the distribution function on time, space (explicitly related to velocity), and momentum (explicitly related to electric field). The Boltzmann equation is valid under assumptions of semi–classical transport: effective mass approximation (which incorporates the quantum effects due to periodicity of the crystal); Born approximation for the collisions, in the limit of small perturbation for the electron–phonon interaction and instantaneous collisions; no memory effects, i.e. no dependence on initial condition terms. The phonons are usually treated as in equilibrium, although the condition of non–equilibrium phonons may be included through an additional equation. Analytical solutions of the Boltzmann equation are possible only under very restrictive assumptions. Direct numerical methods for device simulation have been limited by the complexity of the equation, which in the complete 3–D time–dependent form requires seven independent variables for time, space and momentum. In recent times, more powerful computational platforms have spurred a renewed interest in numerical solutions based on the spheroidal harmonics expansion of the distribution function. To–date, most semiconductor applications have been based on stochastic solution methods (Monte Carlo), which involve the simulation of particle trajectories rather than the direct solution of partial differential equations. The vast majority of device simulations are normally based on the numerical solution of approximate models which are related to the Boltzmann equation, coupled to Poisson’s equation for self–consistency. In the simplest approach, the collision term on the right hand side of (1) is substituted with a phenomenological term feq − f (r, k, t) (2) τ where feq indicates the (local) equilibrium distribution function, and τ is a microscopic relaxation time. It is very useful to express the distribution function in terms of velocity, rather than momentum, since it will be easier to calculate electrical currents. In equilibrium we may use the Maxwell–Boltzmann distribution function 1 feq (r, v) = n(r) 2πkB To m∗ −3/2 m∗ |v|2 exp − 2kB To ! (3) where n(r) is the carrier density, To is the lattice temperature, kB is the Boltzmann constant, and m∗ is the effective mass. The use of (3) for semiconductors is justified in equilibrium as long as degeneracy is not present. The carrier density n(r) is directly related to the distribution function according to Z n(r) = dvf (r, v) (4) which is of general applicability. The significance of the momentum relaxation time can be understood if the electric field is switched off instantaneously and a space–independent distribution is considered. The resulting Boltzmann equation is then ∂f feq − f = (5) ∂t τ which shows that the relaxation time is a characteristic decay constant for the return to the equilibrium state. 2 Drift–diffusion model The popular drift-diffusion current equations can be easily derived directly from the Boltzmann equation. Let’s consider a steady state situation and for simplicity a 1–D geometry. With the use of a relaxation time approximation as in (2) the Boltzmann equation becomes eE ∂f ∂f feq − f (v, x) +v = (6) ∗ m ∂v ∂x τ Here, we have used the relation m∗ v = h̄k, which is valid for a parabolic energy band. Note that the charge e has to be taken with the proper sign of the particle (positive for holes and negative for electrons). A general definition of current density is given by Z J(x) = e vf (v, x)dv (7) where the integral on the right hand side represents the first moment of the distribution function. This definition of current can be related to (6) after multiplying both sides by v and integrating over v. From the left hand side we get 1 J(x) vfeq dv − vf (v, x)dv = − (8) τ eτ The equilibrium distribution function is symmetric in v, and its integral is zero. Therefore, we have from (6) Z eτ J(x) = −e ∗ E m Integrating by parts we have Z Z Z ∂f d v dv − eτ ∂v dx ∂f v dv = [vf (v, x)]∞ −∞ − ∂v | {z } 0 2 Z Z v 2 f (v, x)dv (9) f (v, x)dv = −n(x) (10) and we can write Z D v 2 f (v, x)dv = n(x) v 2 E (11) where v 2 is the average of the square of the velocity defined as D v 2 E 1 = n Z v 2 f (v, x)dv (12) For a purely 1–D treatment, the − 23 exponent in (3) may be replaced with − 12 , while the appropriate thermal kinetic energy becomes kB2T instead of 3k2B T . 2 eτ The drift–diffusion equations are derived introducing the mobility µ = m ∗ and replacing v BT with its average equilibrium value km ∗ , therefore neglecting thermal effects. The diffusion coefficient D = µkBe To (Einstein’s relation) is also introduced, and the resulting drift–diffusion current is dn dx dp Jp = qp(x)µp E(x) − qDp dx Jn = qn(x)µn E(x) + qDn (13) (14) where q is used to indicate the absolute value of the electronic charge. Although no direct assumptions on the non–equilibrium distribution function f (v, x) was made in the derivation of (13) and (14), in effect, the choice of equilibrium (thermal) velocity means that the drift–diffusion equations are only valid for very small perturbations of the equilibrium state (low fields). The validity of the drift–diffusion equations is empirically extended by introducing field–dependent mobility µ(E) and diffusion coefficient D(E), obtained from empirical models or detailed calculations. 3 Physical Limitations on Numerical Drift–Diffusion Schemes The complete drift–diffusion model is based on the following set of equations 1. Current equations Jn = qnµn E + qDn ∇n (15) Jp = qpµp E − qDp ∇p (16) 2. Continuity equations ∂n 1 = ∇·Jn + Un ∂t q ∂p 1 = − ∇·Jp + Up ∂t q 3. Poisson’s equation ∇2 V = + q(n − p + NA− − ND ) (17) (18) (19) where Un and Up are the net generation–recombination rates. The continuity equations (17) and (18) are the conservation laws for the carriers. A numerical scheme which solves the continuity equations should 3 1. Conserve the total number of particles inside the device simulated. 2. Respect local positivity of carrier density. Negative density is unphysical. 3. Respect monotonicity of the solution (i.e. it should not introduce spurious space oscillations). Conservative schemes are usually achieved by subdivision of the computational domain into patches (boxes) surrounding the mesh points. The currents are then defined on the boundaries of these elements, thus enforcing conservation (the current exiting one element side is exactly equal to the current entering the neighboring element through the side in common). For example, on a uniform 2–D grid with mesh size ∆, the electron continuity equation could be discretized in an explicit form as follows n(i, j; k + 1) − n(i, j; k) ∆t = J x (i + 21 , j; k) − J x (i − 12 , j; k)) q∆ 1 y J (i, j + 2 ; k) − J y (i, j − 21 ; k) + + U (i, j; k) q∆ (20) This simple approach has certainly practical limitations, but is sufficient to exemplify the idea behind the conservative scheme. With the present convention for positive and negative components, it is easy to see that in the absence of generation-recombination terms, the only contributions to the overall device current arise from the contacts. Remember that, since electrons have negative charge, the particle flux is opposite to the current flux. The actual determination of the current densities appearing in (20) will be discussed later. When the equations are discretized, using finite differences for instance, there are limitations on the choice of mesh size and time step: 1. The mesh size ∆x is limited by the Debye length. 2. The time step is limited by the dielectric relaxation time. A mesh size must be smaller than the Debye length where one has to resolve charge variations in space. A simple example is the carrier redistribution at an interface between two regions with different doping levels. Carriers diffuse into the lower doped region creating excess carrier distribution which at equilibrium decays in space down to the bulk concentration with approximately exponential behavior. The space decay constant is the Debye length s LD = kB T q2N (21) where N is the doping. In GaAs and Si, at room temperature the Debye length is approximately 400 Å when N ≈ 1016 cm−3 and decreases to about only 50 Å when N ≈ 1018 cm−3 . The dielectric relaxation time is the characteristic time for charge fluctuations to decay under the influence of the field that they produce. The dielectric relaxation time may be estimated using tdr = qN µ (22) To see the importance of respecting the limitations related to the dielectric relaxation time, imagine to have a space fluctuation of carrier concentration which is going to relax to equilibrium with a law 4 ∂∆n ∆n(t = 0) =− ∂t tdr (23) A finite difference discretization of this equation gives at the first time step ∆n(∆t) = ∆n(0) − ∆t∆n(0) tdr (24) Clearly, if ∆t > tdr , the value of ∆n(∆t) is negative, which means that the actual concentration is evaluated to be less than the equilibrium value, and it is easy to see that the solution at higher time steps will decay oscillating between positive and negative values of ∆n. An excessively large ∆t may lead, therefore, to nonphysical results. In the case of high mobility the dielectric relaxation 2 time can be very small. For instance, a sample of GaAs with a mobility of 6, 000 Vcm−s and doping 1018 cm−3 has approximately tdr ≈ 10−15 s, and in a simulation the time step ∆t should be chosen to be considerably smaller. 4 Steady State Solution of Bipolar Semiconductor Equations The general semiconductor equations may be rewritten as ∇·(∇V ) = q(n − p + NB ) ∂n ∇·Jn = qU (n, p) + q ∂t ∂p ∇·Jp = −qU (n, p) + q ∂t kB T Jn = qµn −n∇V + ∇n q kB T Jp = qµp −p∇V − ∇p q (25) (26) (27) (28) (29) with NB = NA − ND . Note that the above equations are valid in the limit of small deviations from equilibrium, since the Einstein relations have been used for the diffusion coefficient, normally valid for low fields or large devices. The generation–recombination term U will be in general a a function of the local electron and hole concentrations, according to possible different physical mechanisms, to be examined later in more detail. We will consider from now on steady state, with the time dependent derivatives set to zero. The semiconductor equations constitute a coupled nonlinear set. It is not possible, in general, to obtain a solution directly in one step, but a nonlinear iteration method is required. The two more popular methods for solving the discretized equations are the Gummel’s iteration method and the Newton’s method. It is very difficult to determine an optimum strategy for the solution, since this will depend on a number of details related to the particular device under study. There are in general three possible choices of variables. 1. Natural variable formulation (V, n, p) 2. Quasi–Fermi level formulation (V, φn , φp ), where the quasi–Fermi levels derive from the following definition of carrier concentration out of equilibrium (non–degenerate case) 5 q(V − φn ) kB T q(φp − V ) p = ni exp kB T n = ni exp (30) (31) 3. Slotboom formulation (V, Φn , Φp ) where the Slotboom variables are defined as −qφn = ni exp k T B qφp = ni exp kB T Φn Φp (32) (33) The Slotboom variables are therefore related to the carrier definitions (30) and (31), and the extension to degenerate conditions is cumbersome. Normally, there is a preference for the quasi–Fermi level formulation in steady state simulation, and for the natural variables n and p in transient simulation. 4.1 Normalization and Scaling For the sake of clarity, all formulae have been presented without the use of simplifications or normalization. It is however common practice to perform the actual calculation using normalized units to make the algorithms more efficient, and in cases to avoid numerical overflow and underflow. It is advisable to input the data in M.K.S. or practical units (the use of centimeters is for instance very common in semiconductor practice, instead of meters) and then provide a conversion block before and after the computation blocks to normalize and denormalize the variables. It is advisable to use consistently scaling, rather than set certain constants to arbitrary values. In computational physics, it is common practice to use h̄ = 1, for instance. This simplifies a lot the derivation of formulae, but if this is not done carefully, when denormalization is applied to recover quantitative results in M.K.S. or other appropriate units, it may be very difficult to verify which terms are multiplied by h̄, h̄2 , h̄3 , etc. Problems may arise when formulae are taken from the published literature, since c.g.s. and M.K.S. units are practically interchangeable, unless the dielectric constant is involved. In c.g.s. units, the vacuum permittivity is o = 1 or 1/4π but it is about 8.85 × 10−12 [F/m] in M.K.S. units. A simple change of centimeters into meters and grams into kilograms, will leave many orders of magnitude unaccounted for when the dielectric permittivity is involved. According to our experience, the majority of mistakes in computational applications prepared by beginners are caused by improper scaling or choice of units. Great care and systematic dimensional analysis of the formulae is suggested, at least until a good feeling is developed for the orders of magnitudes of the variables involved in a computation. The most common scaling factors for normalization of semiconductor equations are listed in Table I. 4.2 Gummel’s Iteration Method Gummel’s method solves the equations with a decoupled procedure. If we choose the quasi–Fermi level formulation, we solve first a nonlinear Poisson’s equation. The potential obtained is substituted into the continuity equations, which are now linear, and are solved directly to conclude the 6 TABLE I - Scaling Factors Space Potential Carrier concentration Diffusion coefficient Intrinsic Debye length (N = ni ) Extrinsic Debye length (N = Nmax ) Thermal voltage Intrinsic concentration Maximum doping concentration Practical unit Maximum diffusion coeff. L= q kB T q2 N V ∗ = kBq T N = ni N = Nmax 2 D = 1 cm s D = Dmax M = VD∗ R = DN L2 L2 T = D Mobility Gen–Recomb Time iteration. The result in terms of quasi–Fermi levels is then substituted back into Poisson’s equation until convergence is reached. In order to check for convergence, one can calculate the residuals obtained by positioning all the terms to the left hand side of the equations and substituting the variables with the iteration values. For exact solution values, the residuals should be zero. Convergence is assumed when the residuals are smaller than a set tolerance. The rate of convergence of the Gummel method is faster when there is little coupling between the different equations. The computational cost of one Gummel iteration is one matrix solution for each carrier type plus one iterative solution for the linearization of Poisson’s equation. Note that in conditions of equilibrium (zero bias) only the solution of Poisson’s equation is necessary, since the equilibrium Fermi level is constant and coincides with both quasi–Fermi levels. We give some examples of the quasi–linearization of Poisson equation, as necessary when Gummel’s method is implemented. Let’s consider the 1–D case in equilibrium first. As said earlier, one has to solve only Poisson’s equation, since exact expressions for the carrier concentrations are known. In the non–degenerate case, (30) and (31) are substituted into Poisson’s equation to give d2 V q qV qV = ni exp(−qφn ) exp − ni exp(qφp ) exp − + NA − ND (34) 2 dx kB T kB T In equilibrium it is φn = φp = 0 (Fermi level taken as reference for the potential energy). Furthermore, the equation may be scaled by using the (minimum) extrinsic Debye length for the space coordinate x, and the thermal voltage kB T /q for the potential V . Indicating with V̄ and x̄ for the normalized potential and space coordinates d2 V̄ ni NA − ND = exp(V̄ ) − exp(−V̄ ) + (35) dx̄2 N ni The equilibrium Poisson’s equation (35) can be solved with the following quasi–linearization procedure 1. Set initial guess for the potential V̄ . 2. Write the potential at the next iteration as V̄new = V̄ + δV , and write Poisson’s equation for V̄new with the above substitution d2 V̄ d2 δV ni NA − ND + = exp(V̄ ) exp(δV ) − exp(−V̄ ) exp(−δV ) + 2 2 dx̄ dx̄ N ni 7 (36) 3. Use the linearization exp(±δV ) ≈ 1 ± δV and discretize (36). This equation has a tridiagonal matrix and is readily solved for δV (i). ni 2 δV (i − 1) − 2 + ∆ x[exp(V̄ (i)) − exp(−V̄ (i))] δV (i) + δV (i + 1) = N (37) n1 2 NA − ND ∆ x exp(V̄ (i)) − exp(−V̄ (i)) + N ni −V̄ (i − 1) + 2V̄ (i) − V̄ (i + 1) + 4. Check for convergence. The residual of (35) is calculated and convergence is achieved if the norm of the residual is smaller than a preset tolerance. If convergence is not achieved, return to step 2. In practice one might simply check the norm of the error ||δV ||2 ≤ T ol or ||δV ||∞ ≤ T ol (38) Note that for the solution of Poisson’s equation the boundary conditions are referenced to the equilibrium Fermi level. One may use the separation between the Fermi level and the intrinsic Fermi level at the contacts for the boundary conditions. After the solution in equilibrium is obtained, the applied voltage is increased in steps ∆V ≤ kB T /q. Now the scaled nonlinear Poisson equation becomes d2 V ni NA − ND = exp(−φn ) exp(V ) − exp(φp ) exp(−V ) + 2 dx N ni (39) where the quasi–Fermi levels are also normalized. The continuity equations, as long as the Einstein’s relations are valid, may be written as ∂V kB T ∂ q(V − φn ) −qµn n + qµn ni exp ∂x q ∂x kB T ∂V kB T q ∂V ∂φn −qµn n + qµn n − ∂x q kB T ∂x ∂x ∂φn −qµn n ∂x q(V − φn ) ∂φn −qµn ni exp kB T ∂x qV −kB T ∂ −qφn −qµn ni exp exp kB T q ∂x kB T Jn = = = = = (40) which can be compacted, including quasi–Fermi level normalization as Jn = an (x) ∂ exp(−φn ) ∂x (41) ∂ exp(φp ) ∂x (42) A similar formula is obtained for the holes Jp = ap (x) and the continuity equations are 8 ∂ ∂ an (x) exp(−φn ) = qU (x) ∂x ∂x ∂ ∂ ap (x) exp(φp ) = qU (x) ∂x ∂x (43) (44) The continuity equations may be discretized with a straightforward finite difference approach (here for simplicity with uniform mesh) an (i+ 12 )[Φn (i+1)−Φn (i)] ∆x an (i− 1 )[Φn (i)−Φn (i−1)] 2 − ∆x =U (45) ∆x where the Slotboom variables have been used for simplicity of notation. Note that the inner derivative has been discretized with centered differences around the points (i ± 12 ) of the interleaved mesh. Variables on the interleaved mesh must be determined very carefully, using consistent interpolation schemes for potential and carrier density, as discussed later. The continuity equations give the tridiagonal system 1 1 1 1 an (i − )Φn (i − 1) − [an (i + ) + an (i − )]Φn (i) + an (i + )Φn (i + 1) = 2 2 2 2 ∆2 xU (i) 1 1 1 1 ap (i − )Φp (i − 1) − [ap (i + ) + ap (i − )]Φp (i) + ap (i + )Φp (i + 1) = 2 2 2 2 −∆2 xU (i) (46) (47) The decoupled iteration solves now Poisson’s equation (39), initially with a guess for the quasi– Fermi levels. The voltage distribution obtained for the previous voltage considered is normally a good initial guess for the potential. Since the quasi–Fermi levels are inputs for Poisson’s equation, the quasi–linearization procedure for equilibrium can be used again. The potential is then used to update the an (i) and ap (i), and (46) and (47) are solved to provide new quasi–Fermi level values for Poisson’s equation, and the process is repeated until convergence is reached. The generation– recombination term depends on the electron and hole concentrations, therefore it has to be updated at each iteration. It is possible to update the generation–recombination term also intermediately, using the result of (46) for the electron concentration. The examples given here to illustrate the Gummel’s approach are limited to the nondegenerate case. If field dependent mobility and diffusion coefficients are introduced, minimal changes should be necessary, as long as it is still justified to use of Einstein’s relations. Extension to nonuniform mesh is left as an exercise for the reader. In the 2–D case, the quasi–linearized Poisson’s equation becomes ni − 4+h [Φn (i, j) exp(V (i, j)) + Φp (i, j) exp(−V (i, j))] δV (i, j) N +[δV (i − 1, j) + δV (i + 1, j) + δV (i, j − 1) + δV (i, j + 1)] = 2 4V (i, j) − V (i − 1, j) − V (i + 1, j) − V (i, j − 1) − V (i, j + 1) + NA + NB 2 ni h Φn (i, j) exp(V (i, j)) + Φp (i, j) exp(−V (i, j)) + N ni 9 (48) The normalized mesh size is h = ∆x = ∆y. As before, the thermal voltage kB T /q has been used to normalize the potential V and the quasi–Fermi levels φn and φp included in the Slotboom variables Φn,p = exp(±φn,p ). The continuity equations with the form ∇·(a(x, y)∇Φ) = ±U (x, y) are discretized as 1 1 1 1 1 − a(i + , j) + a(i − , j) + a(i, j + ) + a(i, j − ) Φ(i, j) + a(i + , j)Φ(i + 1, j) 2 2 2 2 2 1 1 1 +a(i − , j)Φ(i − 1, j) + a(i, j + )Φ(i, j + 1) + a(i, j − )Φ(i, j − 1) = ±h2 U (i, j) 2 2 2 4.3 (49) Newton’s Method Newton’s method is a coupled procedure which solves the equations simultaneously, through a generalization of the Newton–Raphson method for determining the roots of an equation. We rewrite (25)–(27) in the residual form WV (V, n, p) = 0 Wn (V, n, p) = 0 WV (V, n, p) = 0 (50) Starting from an initial guess Vo , no , and po , the corrections V , ∆n, and ∆p are calculated from the Jacobian system ∂WV ∂V ∂Wn ∂V ∂Wp ∂WV ∂n ∂Wn ∂n ∂Wp ∂n ∂V ∂WV ∂p ∂Wn ∂n ∂Wp ∂p W ∆V V ∆n = − Wn Wp ∆p (51) which is obtained by Taylor expansion. The solutions are then updated according to the scheme V (k + 1) = V (k) + ∆V (k) n(k + 1) = n(k) + ∆n(k) (52) p(k + 1) = p(k) + ∆p(k) where k indicates the iteration number. In practice, a relaxation approach is also applied to avoid excessive variations of the solutions at each iteration step. The system (51) has 3 equations for each mesh point on the grid. This indicates the main disadvantage of a full Newton iteration, related to the computational cost of matrix inversion (one may estimate that a 3N × 3N matrix takes typically 20 times longer to invert than an analogous N × N matrix!). On the other hand convergence is usually fast for the Newton method, provided that the initial condition is reasonably close to the solution, and is in the neighborhood where the solution is unique. There are several viable approaches to alleviate the computational requirements of the Newton’s method. In the Newton–Richardson approach, the Jacobian matrix in (51) is updated only when the norm of the error does not decrease according to a preset criterion. In general, the Jacobian matrix is not symmetric positive definite, and fairly expensive solvers are necessary. Iterative schemes have been proposed to solve each step of Newton’s method by reformulating (51) as ∂W V ∂V ∂Wn ∂V ∂Wp ∂V 0 ∂Wn ∂n ∂Wp ∂n 0 0 ∂Wp ∂p | W V = − Wn Wp {z } ∆V ∆n ∆p k+1 10 0 − 0 0 ∂WV ∂n 0 0 ∂WV ∂p ∂Wn ∂p 0 | {z } ∆V ∆n ∆p k (53) Since the matrix on the left hand side is lower triangular, one may solve decoupling into three systems of equations solved in sequence. First, one solves the block of equations (again, one for each grid point) ∂WV ∂WV ∂WV (∆V )k+1 = −WV − (∆n)k − (∆p)k ∂V ∂n ∂p and the result is used in the next block of equations (54) ∂Wn ∂Wn ∂Wn (∆n)k+1 = −Wn − (∆V )k+1 − (∆p)k ∂n ∂V ∂p Similarly, for the third block (55) ∂Wp ∂Wp ∂Wp (∆p)k+1 = −Wp − (∆V )k+1 − (∆n)k+1 (56) ∂p ∂V ∂n The procedure achieves a decoupling of the equations as in a block Gauss-Seidel iteration, and can be intended as a generalization of the Gummel method. A block–SOR method is obtained if the left hand sides are premultiplied by a relaxation parameter. This iteration procedure has better performance if the actual variables are (V, φn , φp ). In general, Gummel’s method is preferred at low bias because of its faster convergence and low cost per iteration. At medium and high bias the Newton’s method becomes more convenient, since the convergence rate of Gummel’s method becomes worse as the coupling between equations becomes stronger at hogher bias. But since Gummel’s method has a fast initial error reduction, it is often convenient to couple the two procedures, using Newton’s method after several Gummel’s iterations. Remember that it is very important for the Newton’s iteration to start as close as possible to the true solution. Close to convergence, the residual in Newton’s iteration should decrease quadratically from one iteration to the other. 5 Generation and Recombination The Shockley–Reed–Hall model is very often used for the generation–recombination term due to trap levels USRH = np − n2i h τp n + ni exp q(Et −Ei ) kB T i h + τn p + ni exp q(Ei −Et ) kB T i (57) where Et is the trap energy level involved and τn and τp are the electron and hole lifetimes. Surface 1 rates may be included with a similar formula, in which the lifetimes are substituted by sn,p where sn,p is the surface recombination velocity. The Auger recombination may be accounted for by using the formula UAug = Cn [pn2 − nn2i ] + Cp [np2 − pn2i ] (58) where Cn and Cp are appropriate constants. The Auger effect is for instance very relevant in the modeling of highly doped emitter regions in bipolar transistors. The generation process due to impact ionization can be included using the field–dependent rate a∞ n exp crit βn − EnE |Jn | UI = + q 11 a∞ p exp E crit − pE βp |Jp | (59) 6 Time-dependent simulation The time–dependent form of the drift–diffusion equations can be used both for steady–state and transient calculations. Steady–state analysis is accomplished by starting from an initial guess, and letting the numerical system evolve until a stationary solution is reached, within set tolerance limits. This approach is seldom used in practice, since now robust steady–state simulators are widely available. It is nonetheless an appealing technique for beginners since a relatively small effort is necessary for simple applications and elementary discretization approaches. If an explicit scheme is selected, no matrix solutions are necessary, but it is normally the case that stability is possible only for extremely small timesteps. The simulation of transients requires the knowledge of a physically meaningful initial condition, which can be obtained from a steady-state calculation. The same time–dependent numerical approaches used for steady-state simulation are suitable, but there must be more care for the boundary conditions, because of the presence of displacement current during transients. In a transient simulation to determine the steady–state, the displacement current can be neglected because it goes to zero when a stationary condition is reached. Therefore, it is sufficient to impose on the contacts the appropriate potential values provided by the bias network. In a true transient regime, however, the presence of displacement currents manifests itself as a potential variation at the contacts, superimposed to the bias, which depends on the external circuit in communication with the contacts. Neglect of the displacement current in a transient is equivalent to the application of bias voltages using ideal voltage generators, with zero internal impedance. In such a situation, the potential variations due to displacement current drop across a short circuit, and are therefore cancelled. In this arrangement, one will observe the shortest possible switching time attainable with the structure considered, but in practice an external load and parasitics will be present, and the switching times will be normally longer. A simulation neglecting displacement current effects may be useful to assess the ultimate speed limits of a device structure. When a realistic situation is considered, it is necessary to include a displacement term in the current equations. It is particularly simple to deal with a 1–D situation. Consider a 1–D device with length W and a cross–sectional area A. The total current flowing in the device is ∂E(x, t) (60) ∂t The displacement term makes the total current constant at each possition x. This property can be exploited to perform an integration along the device ID (t) = In (x, t) + A 1 ID (t) = W Z W In (x, t)dx + 0 A ∂V ∗ W ∂t (61) where V ∗ (t) is the total voltage drop across the structure, with the ground reference voltage applied at x = W . The term A W is called cold capacitance. The 1–D device, therefore, can be studied as the parallel of a current generator and of the cold capacitance which is in parallel with the (linear) load circuit. At every time step, V ∗ has to be updated, since it depends on the charge stored by the capacitors. To illustrate the procedure, consider a simple Gunn diode in parallel with an RLC resonant load containing the bias source. Calling Co the parallel of cold and load capacitance, it is I(t) = Co dV ∗ (t) + Io (t) dt 12 (62) where Io (t) is the particle current given by the first term on the right hand side of (61), calculated at the given time step with drift–diffusion (or any other suitable scheme). It is also V ∗ (t) V ∗ (t) − Vb − dt R L Upon time differencing this last equation, with the use of finite differences we obtain I(t) = − Z ∆t Co V ∗ (t + ∆t) − V ∗ (t) ∆t I(t + ∆t) = I(t) − − [V ∗ (t) − Vb ] R L V ∗ (t + ∆t) = V ∗ (t) + [I(t) − Io (t)] (63) (64) (65) This set of difference equations allows one to update the boundary conditions for Poisson’s equation at every time step to fully include displacement current. A robust approach for transient simulation should be based on the same numerical apparatus established for purely steady-state models. It is usually preferred to use fully implicit schemes, which require a matrix solution at each iteration, because the choice of the timestep is more likely to be limited by the physical time constants of the problem rather than by stability of the numerical scheme. In order to estimate the timestep limits, let’s assume a typical electron velocity v = 107 cm/s and a spatial mesh ∆x = 0.01µm. The C.F.L. condition necessary to resolve correctly a purely drift process on this mesh requires ∆t ≤ ∆x/v = 10−15 s. As calculated earlier, this value is not too far from typical values of the dielectric relaxation time in practical semiconductor structures. When dealing with unipolar devices, as often used in many microwave applications, it is possible to formulate very simple time–dependent drift–diffusion models, which can be solved with straightforward finite difference techniques and are suitable for small student projects. If we can neglect the generation–recombination effects, the 1–D unipolar drift–diffusion model is reduced to the following system of equations ∂n ∂t 2 d V dx2 d d d [nvd (E)] + D(E) n dx dx dx q(n − ND ) = − (66) = (67) where vd (E) = −µn (E)E is the drift velocity. There are two physical processes involved: drift (advection) expressed by the first term on the right hand side of (66), and diffusion described by the second term. The continuity equation (66) is an admixture of competing hyperbolic and parabolic behavior whose relative importance depends on the local electric field strength. The system (66) and (67) can be used for both transient or steady state conditions if the simulation is run ∂n/∂t = 0. A basic simple algorithm could consists of the following steps 1. Guess the carrier distribution n(x). 2. Solve Poisson’s equation to obtain the field distribution. 3. Compute one iteration of the discretized continuity equation with time step ∆t. v(E) and D(E) are updated according to the local field value. 4. Check for convergence. If convergence is obtained, stop. Otherwise, go back to step (2) updating the charge distribution. 13 This is an uncoupled procedure, since (66) and (67) are not solved simultaneously. Usually, explicit methods are used for computational speed. The time step must respect the limitations due to the C.F.L. condition (related to the advective component) and to the dielectric relaxation time. A simple discretization scheme could employ an explicit finite difference approach n(i; k + 1) = n(i; k) + + ∆t {[vd (i − 1; k)n(i − 1; k) − vd (i; k)n(i; k)] ∆x 1 D(i; k)[n(i − 1; k) − 2n(i; k) + n(i + 1; k)]}; ∆x n(i; k + 1) = n(i; k) + vd < 0 (68) ∆t {[vd (i; k)n(i; k) − vd (i + 1; k)n(i + 1; k)] ∆x 1 D(i; k)[n(i − 1; k) − 2n(i; k) + n(i + 1; k)]}; vd > 0 (69) ∆x where we have introduced upwinding for the drift term and we have assumed that the diffusion coefficient is slowly varying in space. There are of course many other possible explicit and implicit discretizations. Such simple finite difference approaches are in general a compromise which cannot provide at one time an optimal treatment of both advective and diffusive components. Because of spatially varying drift velocity, spurious diffusion and dispersion are present. This could be mitigated by using a nonuniform grid discretization, where the mesh size is locally adapted to achieve vd = ∆x/∆t everywhere, which would involve interpolation to the new grid-points. The discretization for a diffusive process is better behaved with a fully implicit scheme (if the CrankNicholson approach is used, one needs to make sure that spurious oscillations in the solution do not develop). On the other hand, the fully implicit algorithm for advection is not conservative. From these conflicting requirements, it emerges that it would be beneficial to split the drift and diffusion processes, and apply an optimal solution procedure to each. There are 1–D situations where this is known to be nearly exact. In well known experiments, a small concentration of excess carriers is generated in a semiconductor sample with a uniform electric field, and the motion of the centroid of the carrier envelope can be studied independently of the diffusive spread of the spatial distribution around the centroid itself. For an initial Gaussian distribution in space, a simple analytical solution shows that drift and diffusion can be treated as a sequential process, each using the total duration of the observation as simulation time. In analogy with this, the 1–D continuity equation can be solved in two steps, for instance + ∆t ; vd < 0 (70) ∆x n(j, i + 1) = n∗ (j, i + 1) + D(j)[n(j − 1, i + 1) − 2n(j, i − 1) ∆t +n(j + 1, i + 1)] 2 (71) ∆ x where again a simple explicit upwinding scheme is used for the drift, while a fully implicit scheme is used for the diffusion. n∗ (j, i + 1) = n(j, i) + vd (j)[n(j − 1, i) − n(j, i)) 7 Scharfetter-Gummel approximation The discretization of the continuity equations in conservation form, requires the determination of the currents on the mid–points of mesh lines connecting neighboring grid nodes. Since the solutions are accessible only on the grid nodes, interpolation schemes are needed to determine the currents. 14 For consistency with Poisson’s equation, it is common to assume that the potential varies linearly between two neighboring nodes. This is equivalent to assume a constant field along the mesh lines, and the field at the mid–point is obtained by centered finite differences of the potential values. In order to evaluate the current, it is also necessary to estimate the carrier density at the mid–points. The simplest approximation which comes to mind is to also assume a linear variation of the carrier density, by taking the arithmetic average between to neighboring nodes. This simple approach is only acceptable for very small potential variation between the nodes, and indeed is exact only if the field between two nodes is zero, which implies the same exact carrier density on the two points. In order to illustrate this, let’s consider a 1–D mesh where we want to discretize the electron current dψ dn + qDn (72) dx dx Here, the field is explicitly expressed by the derivative of the potential. The discretization on the mid–point of the mesh line between nodes xi and xi+1 is given by Jn = qµn n − ψi+1 − ψi ni+1 − ni + qDn 2 2 ∆x ∆x In the simple approach indicated above, the carrier density ni+ 1 is expressed as Ji+ 1 = − qµn ni+ 1 (73) 2 ni+1 + ni (74) 2 In (73), the assumed linearity of the potential between meshes, is implied by the use of the centered finite differences to express the field on the mid–point. We can now rewrite (73) including the approximation in (74) as ni+ 1 ≈ 2 Ji+ 1 = 2 ni+1 −q µn ψi+1 − ψi Dn + q 2 ∆x ∆x − µn ψi+1 − ψi Dn ni + q q 2 ∆x} | {z∆x } | {z a (75) b If we assume a condition where Jn = 0 (equilibrium) and a >> b (negligible diffusion) it is easy to see that positivity of the carrier density is not guaranteed, since the solution oscillates as ni+1 ≈ −ni . Also, it can be shown that for stability we need to have ψi+1 − ψi > 2kB T /q, which requires very small mesh spacing to be verified. The approach by Scharfetter and Gummel (1969) has provided an optimal solution to this problem, although the mathematical properties of the proposed scheme have been fully recognized much later. We consider again a linear potential variation between neighboring mesh points, which is consistent with the use of finite differences to express the field. We express the current in the interval [xi ; xi+1 ] as a truncated expansion about the value at the mid–point xi+ 1 2 ∂ Jn (x) (76) 2 2 ∂x From (76) we obtain a first order differential equation for Jn which can be solved to provide n(x) in the mesh interval, using as boundary conditions the values of carrier density ni and ni+1 . We obtain Jn (x) = Jn (xi+ 1 ) + (x − xi+ 1 ) 15 n(x) = [1 − g(x, ψ)] ni + g(x, ψ)ni+1 ; x ∈ [xi ; xi+1 ] (77) where g(x, ψ) is the growth function ψi+1 − ψi x − xi 1 − exp kB T /q ∆x g(x, ψ) = ψi+1 − ψi / 1 − exp kB T /q (78) The result in (77) can be used to evaluate n(xi+ 1 ) for the discretization of the current in (73). It 2 is easy to see that only when ψ(i + i) − ψ(i) = 0 we have 1 1 ni + ni+1 n1+ 1 = (1 − )ni + ni+1 = (79) 2 2 2 2 The continuity equation can be easily discretized on rectangular uniform and nonuniform meshes using the above results for the currents, because the mesh lines are aligned exactly. For a more in depth reading on advanced aspects of drift–diffusion simulation, we recommend to consult the book by S. Selberherr, cited below. Some classic references: H.K. Gummel, “A self–consistent iterative scheme for one–dimensional steady state transistor calculation,” IEEE Transactions on Electron Devices, vol. ED-11, pp.455–465, 1964. A. DeMari, “An accurate numerical steady state one–dimensional solution of the p–n junction,” Solid–state Electronics, vol. 11, pp. 33–59, 1968. A. DeMari, “An accurate numerical one–dimensional solution of the p–n junction under arbitrary transient conditions,” Solid–state Electronics, vol. 11, pp. 1021–2053, 1968. S. Selberherr, Analysis and Simulation of Semiconductor Devices, Springer–Verlag, 1984. D. L. Scharfetter and D. L. Gummel, “Large signal analysis of a Silicon Read diode oscillator,” IEEE Transaction on Electron Devices, vol. ED-16, pp.64–77, 1969. J. W. Slotboom, “Iterative scheme for 1 and 2–dimensional d.c. transistor simulation,” Electronics Letters, vol. 5, pp. 677–678, 1969. J.W. Slotboom, “Computer–aided two–dimensional analysis of bipolar transistors,” IEEE Transactions on Electron Devices, vol. ED-20, pp. 669–679, 1973. 16