A Modified Characteristics-mixed Finite Element for Semiconductor Device of Heat Conduction and Numerical Analysis

Numerical simulation of a three-dimensional semiconductor device of heat conduction is a fundamental problem in modern information science. The mathematical model is formulated by a nonlinear system of initial-boundary problem, which is interpreted by four partial differential equations: an elliptic equation for electrostatic potential, two convection-diffusion equations for electron concentration and hole concentration, a heat conduction equation for temperature. The electrostatic potential appears within the latter three equations, and the electric field strength controls the concentrations and the temperature. The electric field potential is solved by a mixed finite element method, and the electric field strength is obtained simultaneously. The first order of the accuracy is improved for the latter. The concentrations and temperature are computed by the characteristics-finite element method, where the characteristic approximation is adopted for the hyperbolic term and finite element method is use to treat the diffusion. The composite computational scheme can solve the convection-dominated diffusion equations well because it can cancel numerical dispersion and nonphysical oscillation. The temperature is computed by finite element method, and an interesting simulation tool is proposed for solving semiconductor device problem numerically. By using the technique of a priori estimates of differential equations, an optimal order error estimates is obtained. A theoretical work is shown for numerical simulation of information science, and the actual problem is solved well.


Introduction
In this paper we discuss numerical simulation of a three-dimensional semiconductor device of heat conduction, a fundamental problem in information science.Its mathematical model is formulated by four nonlinear partial differential equations: 1) an elliptic equation for electric potential, 2) a convection-diffusion equation for electron concentration, 3) a convection-diffusion equation for hole concentration, 4) a heat conduction equation for temperature.The electric potential appears within the latter three equations, and the electric field strength controls the concentrations and the temperature.The nonlinear partial differential system with initial-boundary conditions on a three-dimensional domain Ω is defined as follows (Bank, Coughran, Fichtner, Grosse, Rose & Smith, 1985;Jerome, 1994;Lou, 1995;Yuan, 1996), −∆ψ = α(p − e + N(X)), X = (x, y, z) T ∈ Ω, t ∈ J = (0, T ], (1) ρ ∂T ∂t − ∆T = { (D p (X)∇p + µ p (X)p∇ψ) − (D e (X)∇e − µ e (X)e∇ψ) } • ∇ψ, (X, t) ∈ Ω × J. (4) The electric potential, electron concentration, hole concentration and temperature are the objective functions, denoted by ψ, e, p and T , respectively.All the coefficients of ( 1)-( 4) are bounded.α = q/ε, where q and ε are positive constants denoting the electronic load and the permittivity, respectively.U T is the thermal voltage.The diffusion D s (X) depends on the mobility µ s (X), i.e., D s (X) = U T µ s (X), (s = e for the electron and s = p for the hole).N D (X) and N A (X) are the donor impurity concentration and acceptor impurity concentration, respectively.N(X), defined by N(X) = N D (X) − N A (X), changes rapidly as X approaches nearby the P-N junction.R 1 (e, p, T ) and R 2 (e, p, T ) are the recombination rates of the electron, hole and temperature.ρ(X) is the heat transfer coefficient.A nonuniform partition is adopted usually in numerical simulation (He, 1989;Shi, 2002;Yuan, 2009Yuan, , 2013)).
In this paper we consider the second type boundary condition (Neumann boundary condition) mainly: where ∂Ω is the boundary of Ω, and γ is the unit outer normal vector of ∂Ω.
A compatibility condition is added ∫ and the following condition is given to avoid the ambiguous solution Numerical simulation of a semiconductor device is important and valuable in manufacturing modern semiconductor (He, 1989;Shi, 2002;Yuan, 2009Yuan, , 2013)).Gummel proposes the sequence iteration to compute the semiconductor problem in 1964 and states a new problem of numerical simulation in semiconductor device (Gummel, 1964).Douglas and Yuan put forward a simple but useful finite difference method and discuss the application and numerical analysis first for the onedimensional and two-dimensional preliminary problems (constant coefficients and without temperature effect) (Douglas & Yuan, 1987;Yuan, Ding & Yang, 1982), and the research becomes basic theoretical work in numerical simulation of semiconductor device problem.Yuan discusses the characteristic finite element method for variable coefficient problem (Yuan, 1993).Since the diffusion only includes the electric-field strength −∇ψ, Yuan presents the characteristics-mixed finite element where the concentration is obtained by characteristic finite element and the potential is solved by mixed finite element, and derives optimal-order error estimates in H 1 -norm and L 2 -norm for a semidiscrete scheme and a fullydiscrete method (Yuan, 1991, 19912).There the authors only consider a two-dimensional problem without heat factor, and the index k of mixed finite element space and l of finite element space are restricted by k ≥ 1 and l ≥ 1.These features should be improved for actual applications (Jerome, 1994;Yuan, 2009Yuan, , 2013)).Then a characteristic finite element and a characteristic finite difference method are proposed on a uniform partition for a three-dimensional semiconductor device problem of heat conductor (Yuan, 19912, 2000).In this paper the authors put forward a mixed finite elementcharacteristic mixed finite element for solving a three-dimensional semiconductor device problem, where the potential, concentrations and temperature are computed by a mixed finite element, characteristics-finite element and finite element approximation, respectively.Suppose that k ≥ 0 and l ≥ 1.By applying a priori estimates theory and special techniques of differential equations, we obtain optimal-order error estimates in L 2 norm.This composite numerical method shows important suggestions in solving semiconductor problem such as numerical method, software design, actual applications and theoretical and physical study (He, 1989;Jerome, 1994;Shi, 2002;Yuan, 2009Yuan, , 2013)).
Equations ( 10) and ( 11) are solved by using MMOC.For convenience to obtain the approximations at the boundary, we suppose that Ω is a cube and the problem of ( 9)-( 14) is Ω-periodic, i.e. all the functions are Ω-periodic.This assumption is reasonable in physical science because mirror reflection is used according to non-permeation condition.Furthermore, the boundary condition (13) can be omitted because it has a small effect on the interior flow of oil reservoirs in common numerical simulations (Ewing, 1984;He, 1989;Shi, 2002;Yuan, 2009Yuan, , 2013)).
Introduce Sobolev space and norms on Ω, and The time-dependent space are given.Let [a, b] ⊂ J and let X denote a space above.For a function f T ] and X = Ω, we omit J and Ω, and replace L ∞ (0, T ; Suppose that the problem of ( 9)-( 14) is regular, (R) Here l ≥ 1 and k ≥ 0 are integers, and they denote the degrees of polynomials approximating e, p, T and ψ.

The Procedures of MFEM-MMOC
Numerical scheme of ( 9)-( 14) is constructed in this section.A mixed finite element method (MFEM), a modified method of characteristics (MMOC) and a finite element method (FEM) are used to solve the potential, the electron and hole concentrations, and the temperature, respectively.

The MFEM for Potential
Eq. ( 1) is reformulated by ( 9a) and (9b).Test both sides of (9b) by v ∈ H(div), obtain an inner product equation on Ω and apply the divergence theorem for the term dependent on ∇ψ.Eq. ( 9a) is treated similarly.Then, we have a saddle-point problem of (9 Introduce the time-dependent partition and finite element space.Let Ω = ∪ Ω e be a quasiuniform partition with diameter

The MMOC for Concentrations
Generally, the strength changes slower than the concentrations and the computational cost in space and time is more expensive when the MFEM (18) used at each time level.Thus, a larger time step is adopted for the potential equation.
For convenience, the partition for concentrations is obtained by refining the meshes of the potential, 0 During the computation of MMOC, we should define an extrapolation of u h (X, t p m−2 ) and u h (X, t p m−1 ), denoted by Eu h (X, t c n ), to approximate the value at The MMOC is defined now for ( 10) and ( 11).Let denote a unit vector τ e of (−µ e Eu, 1) in the composite space Ω × [0, T ], and let ϕ e (X) It is approximated by using the characteristic finite difference, where Xn−1 The hole concentration equation ( 11) is discussed similarly.τ p is a unit vector of (µ p Eu, 1).Let ϕ p (X) = ( µ p Eu(X) 2 + 1 ) 1/2 , then we have The characteristic derivative is treated by where Xn−1 During the computations, the space step h p for the potential is larger than that for the concentrations and temperature, i.e., h p > h c .S l c ⊂ W 1 ∞ (Ω) denotes a finite element space for the concentrations and temperature consisting of all the piecewise polynomial functions of degree at most l, l ≥ 1. e h (X, 0), p h (X, 0) and T h (X, 0) are approximations of e 0 (X), p 0 (X) and T 0 (X), where Ritz projection and interpolations are used generally.Then MFEM-MMOC for solving the problem of ( 9)-( 14) follows.
(1) Initial approximations are given by e h (X, 0), p h (X, 0), T h (X, 0) X ∈ Ω. (26) (2) Given e h (X, t p m−1 ), p h (X, t p m−1 ) and T h (X, t p m−1 ), numerical solutions at where êh (X, where ph (X, (3) When e h (X, t p m ) and p h (X, t p m ) are computed, then the potential and strength at t p m , u h (X, MFEM-MMOC runs as follows.First ( 26) is determined by using the Ritz projection (see the following section) or interpolations.Then by using ( 27), e h (X, is used to obtain u h (X, t p m ), ψ h (X, t p m ) at t p m .The procedures of ( 27) and ( 28) are used repeatedly and all the numerical solutions are obtained.According to (C), numerical solutions exist and are unique.

Preliminary Estimates
Suppose that the following approximation property and inverse property hold where A 1 and K 1 are positive constants independent of h c .
Similarly, {V k , S k ψ } has the following properties where A 2 and K 2 are positive constants independent of h p .
The Ritz projection of e(X, t), Πe(X, t) ∈ S l c , t ∈ (0, T ], is defined by According to the discussions (Ewing, 1984;Russell, 1980;Wheeler, 1973), we have where A 1 is independent of e and h c .
The Ritz projection of p(X, t), Πp(X, t) ∈ S l c , t ∈ (0, T ], is defined by and it holds where A 1 is independent of p and h c .
The Ritz projection of T (X, t), ΠT (X, t) ∈ S l c , t ∈ (0, T ], is defined by Similarly, where A 1 is independent of T and h c .
The projection of (u, Then it follows from the discussion of Raviart-Thomas (Brezzi, 1974;Douglas, 1983) where A 2 is independent of h p , u, ψ and e, p.
Define an extrapolation of u, where A 3 is a positive constant.

Convergence Analysis
For simplicity, ∆t c and ∆t p are supposed to be independent of the time.Then, MFEM-MMOC is reformulated as follows to compute numerical solutions {e h , p h , where ên−1 Given e 0 h , p 0 h , T 0 h t, by using ( 45), u h,0 and ψ h,0 are obtained.Then from (44) we get {e } by (45).All the numerical solutions are obtained after repeated computations.
In the following discussions, the symbols K and ε denote a generic positive constant and a generic small positive number, respectively.They can take different values at different places.
The optimal error estimates for k ≥ 0 and l ≥ 1 follow.
Theorem 1 Suppose that ψ, u, e, p, T are exact solutions of ( 1)-( 8) and satisfy ).Let ψ h , u h , e h , p h , T h be numerical solutions of ( 43)-( 45).Suppose that k ≥ 0, l ≥ 1, and the following partition constrict holds Then, where K * is a positive constant independent of h c , h p , ∆t 1 p and ∆t p .Proof.The electric potential equation is considered first.Subtracting (39) at t = t m from (28), we get From the Brezzi's discussion of saddle-point problem (Brezzi, 1974;Douglas, 1983) and eq.( 40), we have then we get From ( 40) and ( 49), it is necessary to show optimal error estimates of ∑ s=e,p ∥s h − s∥ L ∞ (L 2 ) to complete the proof of (47).
Let 44a), we obtain error equation of electron concentration Take Z h = ξ n e to get an L 2 -norm result.The first term on the left-hand side is estimated by Then, we have The right-hand terms of (53) are considered.Applying ( 22) to estimate T 1 2 dτ e dX. Thus, Estimate T 2 by T 3 is bounded by T 4 -T 7 are estimated by T 8 , T 9 and T 10 are considered in a similar manner.Let f be defined on Ω. f denotes one of the three functions, e, ξ e or η e .Z denotes the unit vector of Eu n h − Eu n .Then, where Z ∈ [0, 1] and X − X = µ e E(u n h − u n ).Define then, we derive three results from ( 55) From ( 40) and ( 49), we have Since g e (X) is the mean value of the partial derivatives of e n−1 , it can be estimated by e n−1 W 1

∞
. By (55a), we obtain the estimate of T 8 To estimate g η and g ξ , we need to introduce the following induction hypothesis, Noting that and defining the transformation from ( 61), we have where J is a partition element of the potential equation.Using (60), we get It follows from (63), Therefore, To complete the proof, we need show that ξ e,m−i = O(h k+1 p + h l+1 c + ∆t c + (∆t 1 p ) 3/2 + (∆t p ) 2 ).From the discussions (Ewing, 1984;Russell, 1980) and (50), it holds obviously that E(u n − u . Then (66b) is obtained.