Domain Decomposition Modified with Characteristic Finite Element Method for Numerical Simulation of Semiconductor Transient Problem of Heat Conduction

A characteristic finite element algorithm based on domain decomposition is structured in this paper to approximate numerically multi-dimensional semiconductor transient problems of heat conduction. Finite element approximation is presented for the electric field potential equation, and a domain decomposition discretization with characteristic finite element is put forward for the electron concentration equation, hole concentration equation and heat conductor equation. An optimal order error estimate in L2 norm is derived for the coupled system by using some techniques such as variation, domain decomposition, the method of characteristics, the principle of energy, negative norm estimates, induction hypothesis, prior estimates theory and other techniques of partial differential equations. Finally, experimental data consistent with theoretical convergence rate are shown. This type of numerical method is of high computational efficiency and can successfully solve this international problem.


Introduction
Because of the great development of semiconductor device industry described by an initial-boundary value system of diffusion type of nonlinear partial differential equations, numerical simulation must be of high order accuracy and quite efficient, and traditional numerical methods are not considered in actual computation.Then some new modern discretization techniques are introduced to solve multi-dimensional problems with complicated and irregular geometric regions (He & Wei,1989;Shi, 2002;Yuan, 2009Yuan, , 2013)).
Numerical simulation of semiconductor device begins at the use of the sequence iteration presented Gummel successfully in 1964 (Gummel, 1964).Considering actual computations, Douglas and Yuan put forward a useful finite difference method for simplified one-dimensional and two-dimensional models (constant coefficient, without consideration of temperature) and give the precise theoretical results firstly (Douglas & Yuan, 1987;Yuan, Ding & Yang, 1982).Then Yuan applies the method of characteristics respectively with finite element and mixed finite element to solve two-dimensional problems with variable coefficients and gives optimal order error in H 1 norm and L 2 norm (Yuan, 1991(Yuan, , 1993(Yuan, , 1991)).In view of applications, Yuan discusses the effect of heat conduction to semiconductor device and gives the characteristic finite difference and upwind fractional steps finite difference of three dimensional problem on uniform partition and the convergence rate in discrete l 2 norm (Yuan, 1996(Yuan, , 2000(Yuan, , 2005(Yuan, , 2008)).In numerical simulation of modern semiconductor device, the computation scale is huge, and the discretization consists of tens of thousand or several millions nodes, so a powerful parallel tool is introduced (He & Wei,1989;Shi, 2002;Yuan, 2009Yuan, , 2013)).For the simplest parabolic equation, Dawson, Dupont and Du take the lead in discussing Galerkin domain decomposition method and theoretical analysis (Dawson & Dupont, 1992, 1994;Dawson, Du & Dupont, 1991;Dawson & Du, 1990).On the discussion how to simulate threedimensional semiconductor transient problem of heat conduction efficiently, the authors state a type of domain decomposition with modified characteristic finite element and analyze the convergence and optimal order estimates in l 2 norm by using the techniques such as variation, domain decomposition, characteristics, energy norm, negative norm estimates, mathematical induction hypothesis and prior estimate theory and technique of partial differential equation.Numerical data of experimental example are shown consistent with theoretical convergence rate and by which it illustrates that this procedure is more efficient and applied successfully to solve international famous problems (He & Wei,1989;Shi, 2002;Yuan, 2009Yuan, , 2013;;Bank, et al., 1985;Jerome, 1994;Seidmann, 1986;Lou, 1995).It is most valuable in model analysis, numerical method, principle research, theory and application of semiconductor device simulation.This paper is organized as follows.In the first section, mathematical model, physical interpretation and some related research are introduced.Some notations and primary work are given in the second section.In the third section, the authors present the procedures of domain decomposition and modified characteristic finite element.In the fourth section, the authors attempt to give convergence analysis.We examine numerically the accuracy and parallelism of the scheme by an example.Some conclusions and discussions are stated in the last section.In this paper M and ε express general positive constant and general positive small constant, respectively, and they have different meanings at different places.

Mathematical Model and Physical Background
The mathematical model of semiconductor transient problem of heat conduction is described by an initial boundary value system of four partial differential equations (Yuan, 1996(Yuan, , 2000(Yuan, , 2005(Yuan, , 2008;;Bank, et al., 1985;Jerome, 1994;Seidmann, 1986;Lou, 1995), where the potential is defined by an elliptic equation, the concentrations of electron and hole are defined by convection-diffusion equations and the temperature is defined by a heat conduction equation.Electronic field potential is generated by the intensity in the electron equation, the hole concentration equation and the heat conduction equation, and becomes a system with initial and boundary conditions, which is defined on three variables spatial domain Ω as follows, (1) where the electronic field potential ψ, electron concentration e, hole concentration p, and temperature T are unknown functions.All the coefficients of (1)-( 4) are greater than a positive number and less than another positive number.The number α = q/ε denotes the quotient of two positive constants, electronic load q and dielectric coefficient ε.The relation of diffusion coefficient D s (X) and mobility ratio µ s (X) is formulated by , where U T means the quantity of heat (unit: volt).N(X), a given function, is the difference of the donor impurity concentration N D (X) and the acceptor impurity concentration N A (X).The values of N(X) vary quickly as X lies nearby semiconductor knot P-N.R i (e, p, T )(i = 1, 2) represents the recombination rate influenced by the electron concentration, the hole concentration and the temperature.ρ = ρ(X) is the coefficient of heat conduction.Let gradient operator and Laplace operator of a function be denoted by Initial-value conditions are stated as follows, e(X, 0) = e 0 (X), p(X, 0) = p 0 (X), T (X, 0) = T 0 (X), X ∈ Ω. (5) Boundary-value conditions are defined by where ∂Ω denotes the boundary of Ω, and γ is the outward unit normal vector of ∂Ω.
Neumann boundary value conditions determine a family of electric field potential functions differed by a constant.So an additional condition is supplied to get a unique formulation A compatibility condition is given, In general, the problem (1)-( 7) is supposed to be positive definite, where D * , D * , µ * , µ * , ρ * and ρ * are positive constants.R i (e, p, T )(i = 1, 2) are Lipschitz continuous in ε 0 -neighbor centered by (e, p, T ).

Notations and Preparations
For simplification, let the computational domain as in Fig. 1.We use approximation of derivatives of delta function at several points in this work, and for further reference we define two special functions Φ 2 and Φ 4 as follows (Dawson & Dupont, 1992, 1994;Dawson, Du & Dupont, 1991;Dawson & Du, 1990).
Note that if p(x 1 ) is a polynomial of degree at most one, then and if p(x 1 ) is a polynomial of degree at most three, then Definition 1: For any number H ∈ (0, 1 2 ), define a normalized function Let N h, j be a finite-dimensional finite element space H 1 (Ω j )( j = 1, 2), and let N h (Ω) be a finite-dimensional finite element subspace of l degree of L 2 (Ω).Moreover, if a function v comes from N h , then v| Ω j must belong to N h, j .Note that a given function v ∈ N h (Ω) has a well-defined jump [v] on Γ, Definition 2: Define a bilinear mapping Ds (u, v), where u, v ∈ H 1 (Ω j ), j = 1, 2, and D s (X)(s = e, p) are all positive definite functions, D T = 1, and λ s is a positive constant.
Definition 3: An integral operator approximating normal derivative at interior boundary is defined as follows where Φ formulated by ( 12) is used.
Consider a function ψ of H 1 (Ω 1 ) and H 1 (Ω 2 ), and give a special norm, The second term is formulated in another expression, where Continue to find a positive constant M 0 such that Similarly, we have the following estimates for 0 ≤ t ≤ T , where M 1 , M 2 , M 3 are positive constants, m = 2, 4, and ∂u ∂γ denotes the normal derivative of u across interior boundary Γ.

Domain Decomposition Scheme Modified with Characteristic Finite Element
The parallel procedures are illustrated in this section, which consists of a finite element scheme for the electric field potential and a domain decomposition scheme modified with characteristics for the electron concentration, hole concentration and temperature.Firstly, rewrite the potential equation ( 1), Secondly, ( 2)-( 4) are reformulated as follows where u = ∇ψ.
Let W h = W h ψ denote a finite element space of k degree, and let h ψ be the partition step, then we have Noting that the electric field potential changes very slow with respect to t, we can adopt a large step in its calculations.While in the computations of the electron concentration, hole concentration and temperature small steps are adopted.Some notations are given as follows, ∆t c , time step of the concentration equation, ∆t ψ , time step of the potential equation, j = ∆t ψ /∆t c , t n = n∆t c , t m = m∆t ψ , ψ n = ψ(t n ), and ψ m = ψ(t m ).For a potential function ψ(X, t), define a linear extrapolation Eψ n of the values at the closest two time levels to t n , where the subscripts correspond to potential time levels, superscripts to concentration time levels.The extrapolation can be applied for the vector function u = −∇ψ, that is to say Eu is defined similarly to Eψ and has some same properties.
Electric field potential equation ( 20) is approximated by the following finite element scheme Then the gradient of potential is computed by U h,m = −∇ψ h,m .
The electron concentration equation ( 2), hole concentration equation ( 3) and heat conduction equation ( 4) are considered.Multiply them by ω ∈ N h (Ω), apply integration by parts, then we have weak forms of domain decomposition as follows Noting that the flow moves essentially along the characteristics, we apply a modified procedure of characteristics for the first-order hyperbolic part of ( 2) and (3) to overcome numerical dispersion and oscillation.This kind of method has a high order accuracy and a fine stability in numerical simulation and a large time step can be adopted (Ewing 1983;Ewing, Russell & Wheeler, 1984;Douglas & Russell, 1982;Russell, 1985;Douglas & Yuan, 1986;Ewing, Yuan & Li, 1989).Let τ = τ(X, t) denote the unit vector in characteristic direction (−µ e u 1 , −µ e u 2 , −µ e u 3 , 1) and let τ p = τ p (X, t) be the unit vector in (µ p u 1 , µ p u 2 , µ p u 3 , 1)-direction.Let , s = e, p, and compute the characteristic directional derivative by Then (24a) and ( 24b) are changed into Approximate where Similarly, the heat conduction equation is discretized by the method of domain decomposition with finite element, The program works for a time step as follows.Firstly, given the initial approximation {e 0 h , p 0 h , T 0 h }, {ψ h,0 , U h,0 } is computed by the finite element scheme (23).Secondly,  (26).Finally, numerical solutions can be shown after a series of above computation and they exist and are sole by the positive definite condition (C).

Initial approximation is given by
) where Ω = Ω 1 ∪ Ω 2 , and h c denotes the spacial step of finite element space N h (Ω).
Theorem 1 Suppose that exact solutions of ( 1)-( 7) are suitably regular, Adopt the parallel procedure modified characteristic finite element ( 23) and ( 26) on Ω 1 and Ω 2 , and suppose that the relation of partition parameters is where M 1 is a positive constant and k ≥ 1 and l ≥ 1 are the indexes of finite element space.We have where ||g|| L∞ (J;X) = sup n∆t≤T ||g n || X and the constant M * is dependent on ψ, e, p, T and their derivatives.
Proof.The potential function ψ is considered firstly, and ψ h − ψ is estimated.By (20) (t = t m ), ( 27) (t = t m ) and ( 23), we have The electron concentration e is discussed secondly.Subtract (26a) from (24a) (t = t n+1 ), and use (29a The terms on the left-hand side are estimated as follows, ) The terms on the right hand side of (37) are estimated later.Using (19), we have for positive constants Take ∆t c sufficiently small, 2∆t An induction hypothesis is given, Other terms on the right hand side are estimated.Note that and Φ e is bounded because of (42), we have where It follows from negative norm estimate and induction hypothesis ( 42 Collecting ( 38)-( 44), we derive It continues to estimate the hole concentration.Subtract (26b) from (24b) (t = t n+1 ), and use (29b where pn = p n ( Xn p ) and Xn p = X + µ p EU n+1 h ∆t c .Test function ω h = ξ n+1 p is substituted in (46), and it follows from a similar discussion and analysis, Finally, error equation of temperature is derived from (4), Take ω h = π n+1 as test function in (48), Collecting ( 45), ( 47) and ( 49), multiplying both sides of the resulting equation by 2∆t c , summing on n (0 ≤ n ≤ L − 1), and noting that ξ 0 e = ξ 0 p = π 0 = 0, we have Using Gronwall Lemma, Therefore, it follows from ( 35) and ( 51), It remains to testify the induction hypothesis (42).It is right as n is equal to 0 because of ξ 0 e = ξ 0 p = π 0 = 0. Assume (42) holds for any positive integer n between 1 and a given positive integer L − 1.By (51) (52) and the restriction (40), it is easy to see that (42) holds for n = L. Based on the resulting estimates (51) and (52), and (31), Theorem 1 is proved.
The method discussed in this paper can be extended to three-dimensional case, such as Mespet model (see Fig. 2) (He & Wei,1989).The technique of domain decomposition is very important in numerical simulation of actual applications, as shown in Fig. 2, where computational region Ω is divided into five subdomains, Ω = 5 ∪ i=1 Ω i (He & Wei, 1989;Shi, 2002;Yuan, 2009Yuan, , 2013)).

Numerical Example
In this section one numerical example is presented to testify the efficiency of the parallel algorithm discussed above.Consider the following model problem,  The computational interval is decomposed into two subintervals, [0, 1] = [0, 0.5] ∪ [0.5, 1], and interior boundary is Γ = 0.5.Absolute errors of numerical data are shown in Table 1.   1, and error results of approximation to normal derivative at inner boundary, ∂c ∂x (0.5) = e T sin π = 0 are shown in Table 2. Table 2. Error results of approximation to norma derivative at interior boundary B h = 1/40 6.2728E − 16 h = 1/80 5.1736E − 15 h = 1/160 7.8249E − 14 Time costs (unit: second) in numerical computation of domain decomposition method (DDM) are illustrated in Table 3 comparing with another computation strategy, non-domain decomposition method (NDDM).From Table 3, it is easy to conclude that domain decomposition becomes more and more efficient as the partition becomes more refined and the system of algebraic equations becomes larger and larger.

Conclusion and Discussion
In this paper a domain decomposition of finite element is discussed to simulate semiconductor transient problems of heat conductor.Mathematical model, physical basic and international research are introduced in the first section.In the second section some primary elliptic projections, notations and basic properties are stated.Then the
are computed in parallel by the domain decomposition scheme with finite element (26).It continues to get {ψ h,1 , U h,1 } by the values {e j h , p j h , T j h } = {e h,1 , p h,1 , T h,1 } and (23), and as a result, the solutions {e j+1 h , p j+2 h , T j+1 h },• • • , {e h,2 , p h,2 , T h,2 } are obtained in parallel by