DriftDiffusion Modeling Prepared by Dragica Vasileska Professor Arizona

Drift-Diffusion Modeling Prepared by Dragica Vasileska Professor Arizona State University

Outline of the Lecture l l l l Classification of PDEs Why Numerical Analysis? Numerical Solution Sequence Flow-Chart of Equilibrium Poisson Equation Solver Discretization of the Continuity Equation Numerical Solution Techniques for Sparse Matrices Flow-Chart of 1 D Drift-Diffusion Simulator

Classification of PDEs

Classification of PDEs Different mathematical and physical behaviors: Ø Elliptic Type Ø Parabolic Type Ø Hyperbolic Type System of coupled equations for several variables: Ø Time : first-derivative (second-derivative for wave equation) Ø Space: first- and second-derivatives

Classification of PDEs (cont. ) General form of second-order PDEs ( 2 variables)

PDE Model Problems Ø Hyperbolic (Propagation) • Advection equation (First-order linear) • Wave equation (Second-order linear )

PDE Model Problems (cont. ) Ø Parabolic (Time- or space-marching) • Burger’s equation (Second-order nonlinear) (Diffusion / dispersion) • Fourier equation (Second-order linear )

PDE Model Problems (cont. ) Ø Elliptic (Diffusion, equilibrium problems) • Laplace/Poisson (second-order linear) • Helmholtz equation

Well-Posed Problem Numerically well-posed Discretization equations Ø Auxiliary conditions (discretized approximated) Ø the computational solution exists (existence) • the computational solution is unique (uniqueness) • the computational solution depends continuously on the approximate auxiliary data • the algorithm should be well-posed (stable) also •

Boundary and Initial Conditions ØInitial conditions: starting point for propagation problems ØBoundary conditions: specified on domain boundaries to provide the interior solution in computational domain R R s n

Numerical Methods Ø Complex geometry Ø Complex equations (nonlinear, coupled) Ø Complex initial / boundary conditions Ø No analytic solutions Ø Numerical methods needed !!

Why Numerical Analysis?

Coupling of Transport Equations to Poisson and Band-Structure Solvers D. Vasileska and S. M. Goodnick, Computational Electronics, published by Morgan & Claypool , 2006.

Drift-Diffusion Approach Constitutive Equations l Poisson l Continuity l Current Equations Density Equations S. Selberherr: "Analysis and Simulation of Semiconductor Devices“, Springer, 1984.

Poisson/Laplace Equation Solution Poisson/Laplace Equation No knowledge of solving of PDEs With knowledge for solving of PDEs Method of images Theoretical Approaches Poisson Green’s function method Laplace Numerical Methods: finite difference finite elements Method of separation of variables (Fourier analysis)

Numerical Solution Sequence

Numerical Solution Details Governing Equations ICS/BCS Continuous Solutions Discretization Finite-Difference Finite-Volume System of Algebraic Equations Equation (Matrix) Solver Discrete Nodal Values Tridiagonal φi (x, y, z, t) SOR p (x, y, z, t) Finite-Element Gauss-Seidel Spectral Krylov Boundary Element Multigrid Approximate Solution n (x, y, z, t) Hybrid D. Vasileska, EEE 533 Semiconductor Device and Process Simulation Lecture Notes, Arizona State University, Tempe, AZ.

What is next? l MESH l Finite Difference Discretization l Boundary Conditions

MESH TYPE The course of action taken in three steps is dictated by the nature of the problem being solved, the solution region, and the boundary conditions. The most commonly used grid patterns for two-dimensional problems are Common grid patterns: (a) rectangular grid, (b) skew grid, (c) triangular grid, (d) circular grid.

Mesh Size

Example for Meshing

Finite Difference Schemes Before finding the finite difference solutions to specific PDEs, we will look at how one constructs finite difference approximations from a given differential equation. This essentially involves estimating derivatives numerically. Let’s assume f(x) shown below: Estimates for the derivative of f (x) at P using forward, backward, and central differences.

Finite Difference Schemes We can approximate derivative of f(x), slope or the tangent at P by the slope of the arc PB, giving the forward-difference formula, or the slope of the arc AP, yielding the backward-difference formula, or the slope of the arc AB, resulting in the central-difference formula,

Finite Difference Schemes We can also estimate the second derivative of f (x) at P as or Any approximation of a derivative in terms of values at a discrete set of points is called finite difference approximation.

Finite Difference Schemes The approach used above in obtaining finite difference approximations is rather intuitive. A more general approach is using Taylor’s series. According to the wellknown expansion, and Upon adding these expansions, where O( x)4 is the error introduced by truncating the series. We say that this error is of the order ( x)4 or simply O( x)4. Therefore, O( x)4 represents terms that are not greater than ( x)4. Assuming that these terms are negligible,

Finite Difference Schemes Subtracting from We obtain and neglecting terms of the order ( x)3 yields This shows that the leading errors of the order ( x)2. Similarly, the forward and backward difference formula have truncation errors of O( x).

Poisson Equation

Poisson Equation Linearization l The 1 D Poisson equation is of the form:

Renormalized Form LD=sqrt(qni/εVT)

Finite Difference Representation Equilibrium: Non-Equilibrium: n calculated using PM coupling and p still calculated as in equilibrium case (quasi-equilibrium approximation)

Criterion for Convergence There are several criteria for the convergence of the iterative procedure when solving the Poisson equation, but the simplest one is that nowhere on the mesh the absolute value of the potential update is larger than 1 E-5 V. This criterion has shown to be sufficient for all device simulations that have been performed within the Computational Electronics community.

Boundary Conditions

1 D Discretization

2 D Discretization

2 D Discretization (cont’d)

Flow-Chart of Equilibrium Poisson Equation Solver

Discretization of the Continuity Equation

Sharfetter-Gummel Discretization Scheme

(a) Linearized Scheme • Within the linearized scheme, one has that • This scheme can lead to substantial errors in regions of high electric fields and highly doped devices.

(b) Sharfetter-Gummel Scheme

Bernouli Function Implementation Others as Well …

Numerical Solution Techniques for Sparse Matrices

Numerical Methods Objective: Speed, Accuracy at minimum cost Ø Numerical Accuracy (error analysis) Ø Numerical Stability (stability analysis) Ø Numerical Efficiency (minimize cost) Ø Validation (model/prototype data, field data, analytic solution, theory, asymptotic solution) Ø Reliability and Flexibility (reduce preparation and debugging time) Ø Flow Visualization (graphics and animations)

Solution Methods

LU Decomposition (cont’d)

Iterative Methods l Iterative (or relaxation) methods start with a first approximation which is successively improved by the repeated application (i. e. the “iteration”) of the same algorithm, until a sufficient accuracy is obtained. l In this way, the original approximation is “relaxed” toward the exact solution which is numerically more stable. l Iterative methods are used most often for large sparse system of equations, and always when a good approximation of the solution is known. l Error analysis and convergence rate are two crucial aspects of theory of iterative methods.

Error Equation

Jacobi, Gauss-Seidel Methods

SOR Method

Convergence

Convergence (cont’d)

Multi-Grid Method

Coarsening Techniques

Coarsening Techniques (cont’d)

Restriction

Prolongation

The Coarsest Grid Solver l As shown by the algorithm description, the final coarsest grid has just a few grid-points. A typical grid has 3 up to 5 points per axis. l On this grid, usually called W 0, an exact solution of the basic equation Ae = r is required. l The number of grid points is so small on W 0 that any solver can be used without changing the convergence rate in a noticeable way. l Typical choices are a direct solver (LU), a SOR, or even a few iterations of the error smoothing algorithm.

Relaxation Scheme l The relaxation scheme forms the kernel of the multigrid method. l Its task is to reduce the short wavelength Fourier components of the error on a given grid. l The efficiency of the relaxation scheme depends sensitively on details such as the grid topology and boundary conditions. Therefore, there is no single standard relaxation scheme that can be applied. l Two Gauss-Seidel schemes, namely point-wise relaxation and line relaxation, can be considered. The correct application of one or more relaxation methods can dramatically imp-rove the convergence. The point numbering scheme plays also a crucial role.

Relaxation Scheme (cont’d)

Comparison of Relaxation Schemes

Complexity of linear solvers Time to solve model problem (Poisson’s equation) on regular mesh n 1/2 n 1/3 2 D 3 D Sparse Cholesky: O(n 1. 5 ) O(n 2 ) CG, exact arithmetic: O(n 2 ) CG, no precond: O(n 1. 5 ) O(n 1. 33 ) CG, modified IC: O(n 1. 25 ) O(n 1. 17 ) CG, support trees: Multigrid: O(n 1. 20 ) -> O(n 1+ ) O(n 1. 75 ) -> O(n 1. 31 ) O(n)

Validity of the Drift-Diffusion Model

Validity of the Drift-Diffusion Model (cont’d)

Time-Dependent Simulations l The time-dependent form of the drift-diffusion equations can be used both for ¡ steady-state, ¡ transient l Steady-state and calculations. 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.

Time-Dependent Simulations (cont’d) l 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 time-steps.

Time-Dependent Simulations (cont’d) The simulation of transients requires the knowledge of a physically meaningful initial condition. l 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. l

Displacement Current l In a true transient regime, 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.

Displacement Current (cont’d) ¡ 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 this arrangement, one observes the shortest possible switching time attainable with the structure considered. ¡ Hence, a simulation neglecting displacement current effects may be useful to assess the ultimate speed limits of a device structure.

Displacement Current (cont’d) l 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 crosssectional area A. The total current flowing in the device is (1) l The displacement term makes the total current constant at each position x. This property can be exploited to perform an integration along the device (2) 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. l At every time step, Eq. (1) has to be updated, since it depends on the charge stored by the capacitors. l

Example: Gunn Diode

Some General Comments l. A robust approach for transient simulation should be based on the same numerical apparatus established for purely steady-state models. l It is usually preferred to use fully implicit schemes, which require a matrix solution at each iteration, because the choice of the time-step is more likely to be limited by the physical time constants of the problem rather than by stability of the numerical scheme. l Typical calculated values for the time step are not too far from typical values of the dielectric relaxation time in practical semiconductor structures.

Solution of the Coupled DD Equations There are two schemes that are used in solving the coupled set of equations which comprises the Drift-Diffusion model: ¡ Gummel’s method ¡ Newton’s method

Gummel’s Method Gummel’s relaxation method, which solves the equations with the decoupled procedure, is used in the case of weak coupling: • Low current densities (leakage currents, subthreshold regime), where the concentration dependent diffusion term in the current continuity equation is dominant • The electric field strength is lower than the avalanche threshold, so that the generation term is independent of V • The mobility is nearly independent of E The computational cost of the Gummel’s iteration is one matrix solution for each carrier type plus one iterative solution for the linearization of the Poisson Equation

Gummel’s Method (cont’d)

Newton’s Method • The three equations that constitute the DD model, written in residual form are: • Starting from an initial guess, the corrections are calculated by solving:

Newton’s Method (cont’d)

Flow-Chart of 1 D Drift. Diffusion Simulator