An approximation of semiconductor device by mixed finite element method and characteristics-mixed finite element method

Abstract

The mathematical model of a semiconductor device is described by a system of three quasi-linear partial differential equations for initial boundary value problem. The electric potential is governed by an elliptic equation. The electron and hole concentrations are governed by two equations of convection-dominated diffusion type. Standard mixed finite element is used for the electric potential equation. A characteristics-mixed finite element method is presented for the concentration equations. This scheme conserves mass locally. In order to derive the optimal $L^2$-norm error estimates, a post-processing step is included in the approximation to the scalar concentration. Numerical experiment is presented finally to validate the theoretical analysis.

Introduction

In this paper, we will consider the following mathematical model of two-dimensional semiconductor device, which is governed by a quasi-linear coupled system of three partial differential equations with initial and boundary values [1], [2]. One equation of the elliptic type is for the electric potential; two of convection-dominated diffusion type are for the conservation of electron and hole concentrations. The three equations, relevant initial condition and boundary condition make up a closed system.

$$\left\{\begin{array}{c}(\mathtt{a})-\mathrm{\Delta }\psi =\nabla ·u=\alpha (p-e+N(x))\text{,}(x\text{,}t)\in \mathrm{\Omega }\times [0\text{,}T]\text{,}\hfill \\ (\mathtt{b})\frac{\partial e}{\partial t}=\nabla ·[D_e(x)\nabla e+{\mu }_e(x)\mathit{eu}]-R(e\text{,}p)\text{,}(x\text{,}t)\in \mathrm{\Omega }\times J\text{,}\hfill \\ (\mathtt{c})\frac{\partial p}{\partial t}=\nabla ·[D_p(x)\nabla p-{\mu }_p(x)\mathit{pu}]-R(e\text{,}p)\text{,}(x\text{,}t)\in \mathrm{\Omega }\times J\text{,}\hfill \end{array}\right.$$

where $J=(0\text{,}T]$, and $\mathrm{\Omega }$ is a bounded domain in ${\mathcal{R}}^2$. The unknown functions are the electrostatic potential $\psi$, the electron concentration e and the hole concentration p. $u=-\nabla \psi$ is the electric field intensity. $\alpha =\frac{q}{∊}\text{,}q$ and $∊$ are constants(q is the electron charge, $∊$ is the dielectric permittivity). The diffusion coefficients $D_s(x)(s=e\text{,}p)$ are related to the mobilities ${\mu }_s(x)(s=e\text{,}p)$ by the relation $D_s(x)=U_T{\mu }_s(x)$, where $U_T$ is the thermal voltage. $N(x)=N_D(x)-N_A(x)$ is a given function, $N_D(x)$ and $N_A(x)$ being the donor and acceptor impurity concentrations, respectively. $R(e\text{,}p)$ is the recombination term.

We consider the following boundary condition

$$-{\left.\frac{\partial \psi }{\partial \nu }\right|}_{\partial \mathrm{\Omega }}=u·\nu =0\text{,}{\left.\frac{\partial e}{\partial \nu }\right|}_{\partial \mathrm{\Omega }}=0\text{,}{\left.\frac{\partial p}{\partial \nu }\right|}_{\partial \mathrm{\Omega }}=0\text{,}t\in J\text{,}$$

where $\nu$ is the unit outward normal vector on boundary $\partial \mathrm{\Omega }$.

The initial condition is

$$e(x\text{,}0)=e_0(x)\text{,}p(x\text{,}0)=p_0(x)\text{,}x\in \mathrm{\Omega }\text{,}$$

In order that a solution is possible, the compatibility condition [2]

$${\int }_{\mathrm{\Omega }}(p_0-e_0+N)\mathit{dx}=0$$

must be imposed on the data. Note that the electric potential is determined only up to an additive constant; thus, the conditions

$${\int }_{\mathrm{\Omega }}\psi \mathit{dx}=0\text{,}0⩽t⩽T$$

can be applied to suppress the ambiguity.

In reality $D_s(s=e\text{,}p)$ are quite small, so that (1.1b) and (1.1c) are strongly convection-dominated. It is well known that the standard Galerkin or difference scheme applied to the convection-dominated problems does not work well. A variety of numerical techniques have been introduced to obtain better approximations for (1.1), such as characteristic finite element method [3], characteristic finite difference method [4], upwind finite volume method [5], etc.

The modified method of characteristic finite element method (MMOC-Galerkin) [5], [6] has much smaller numerical diffusion, nonphysical oscillations and time-truncation than those of standard methods, and can be used with a large time step, with corresponding improvement in efficiency and without cost in accuracy. Unfortunately, it fails to preserve local mass balance. Since it uses a Galerkin spatial discretization, local constants are not in the space of test functions. Mixed finite element method has been proven to be an effective numerical method for solving fluid problems. It has an advantage of approximating the unknown variable and its diffusive flux simultaneously.

Arbogast and Wheeler [7] defined a characteristics-mixed method for approximating the solution to an advection-dominated transport problem. It was based on a space-time variational form of the advection-diffusion transport problem. It used a characteristic approximation that is similar to that of MMOC-Galerkin method to handle advection in time and a lowest-order mixed finite element spatial approximation for the diffusion term. Piecewise constants were in the space of test functions, so mass is conserved element by element. A post-processing step was included in the schemes to improve the rate of convergence of the method. It should be pointed out that the schemes of this characteristics-mixed method include many integrals of the test function’s mappings, which are difficult to be evaluated and used in practical computation.

Sun and Yirang Yuan [8] defined an approximation scheme for the incompressible miscible displacement in porous media. Standard mixed finite element is used for the Darcy velocity equation. A characteristics-mixed finite element method is used for the concentration equation. Characteristic approximation is applied to handle the convection part of the concentration equation, and a lowest-order mixed finite element spatial approximation is adopted to deal with the diffusion term. The scheme conserves mass globally and is much easier to compute.

The research of the transient behavior of semiconductor devices is an important field in modern computational mathematics, which is of great value in theory as well as in practice [9], [10]. In this paper, we will consider a combined numerical method for (1.1). Note that electrostatic potential $\psi$ does not appears explicitly in (1.1b) and (1.1c); however, the electric field intensity u does. Consequently, standard mixed finite element is used for the electrostatic potential equation to approximate the electric potential and the electric field intensity simultaneously. The direct evaluation of the electric field intensity can be expected to give improved accuracy for the same computational effort. A characteristics-mixed finite element method is presented for approximating the electron and hole concentrations. Characteristic approximation is applied to handle the convection terms, and a lowest-order mixed finite element spatial approximation is adopted to deal with the diffusion terms. The scheme of approximating the concentrations is locally conservative; in fact, on the discrete level, current is transported along the approximation characteristics. As in [7], [8], a post-processing step is included in the scheme in which the approximations to the scalar unknowns are improved by utilizing the approximate vector flux. This has the effect of improving the rate of convergence of the method.

In this paper, we use the following Sobolev spaces and the norms associated with these spaces:

$$\begin{array}{cc}\hfill & L^2(\mathrm{\Omega })=\left\{f:{\int }_{\mathrm{\Omega }}|f{|}^2\mathit{dx}<\infty \right\}\text{,}\Vert f\Vert ={\left[{\int }_{\mathrm{\Omega }}|f{|}^2\mathit{dx}\right]}^{1/2}\text{,}\hfill \\ \hfill & L^{\infty }(\mathrm{\Omega })=\{f:\mathtt{ess}{\displaystyle \underset{\mathrm{\Omega }}{\sup }}|f|<\infty \}\text{,}\Vert f{\Vert }_{L^{\infty }}=\mathtt{ess}{\displaystyle \underset{\mathrm{\Omega }}{\sup }}|f|\text{,}\hfill \\ \hfill & H^m(\mathrm{\Omega })=\left\{f:\frac{{\partial }^{|\alpha |}f}{\partial x^{\alpha }}\in L^2(\mathrm{\Omega })\text{,}|\alpha |⩽m\right\}\text{,}\Vert f{\Vert }_m={\left[{\displaystyle \sum _{|\alpha |⩽m}}{‖\frac{{\partial }^{|\alpha |}f}{\partial x^{\alpha }}‖}^2\right]}^{1/2}\text{,}m⩾0\text{,}\hfill \\ \hfill & W_{\infty }^m(\mathrm{\Omega })=\left\{f:\frac{{\partial }^{|\alpha |}f}{\partial x^{\alpha }}\in L^{\infty }(\mathrm{\Omega })\text{,}|\alpha |⩽m\right\}\text{,}\Vert f{\Vert }_{W_{\infty }^m}={\displaystyle \underset{|\alpha |⩽m}{\max }}{‖\frac{{\partial }^{|\alpha |}f}{\partial x^{\alpha }}‖}_{L^{\infty }}\text{,}m⩾0\text{.}\hfill \end{array}$$

In particular, $H^0(\mathrm{\Omega })=L^2(\mathrm{\Omega })$, $W_{\infty }^0(\mathrm{\Omega })=L^{\infty }(\mathrm{\Omega })$. Let $[a\text{,}b]\subset J$ and let X be any of the spaces just defined. If $f(x\text{,}t)$ represents functions on $\mathrm{\Omega }\times [a\text{,}b]$, we set

$$\begin{array}{cc}\hfill & H^m(a\text{,}b\text{;}X)=\left\{f:{\int }_a^b{‖\frac{{\partial }^{\alpha }f}{\partial t^{\alpha }}(·\text{,}t)‖}_X^2\mathit{dt}<\infty \text{,}\alpha ⩽m\right\}\text{,}\hfill \\ \hfill & \Vert f{\Vert }_{H^m(a\text{,}b\text{;}X)}={\left[{\displaystyle \sum _{\alpha =0}^m}{\int }_a^b{‖\frac{{\partial }^{\alpha }f}{\partial t^{\alpha }}(·\text{,}t)‖}_X^2\mathit{dt}\right]}^{1/2}\text{,}m⩾0\text{,}\hfill \\ \hfill & W_{\infty }^m(a\text{,}b\text{;}X)=\left\{f:\mathtt{ess}{\displaystyle \underset{[a\text{,}b]}{\sup }}{‖\frac{{\partial }^{\alpha }f}{\partial t^{\alpha }}(·\text{,}t)‖}_X<\infty \text{,}\alpha ⩽m\right\}\text{,}\hfill \\ \hfill & \Vert f{\Vert }_{W_{\infty }^m(a\text{,}b\text{;}X)}={\displaystyle \underset{0⩽\alpha ⩽m}{\max }}\mathtt{ess}{\displaystyle \underset{[a\text{,}b]}{\sup }}{‖\frac{{\partial }^{\alpha }f}{\partial t^{\alpha }}(·\text{,}t)‖}_X\text{,}m⩾0\text{,}\hfill \\ \hfill & L^2(a\text{,}b\text{;}X)=H^0(a\text{,}b\text{;}X)\text{,}{L}^{\infty }(a\text{,}b\text{;}X)=W_{\infty }^0(a\text{,}b\text{;}X)\text{.}\hfill \end{array}$$

If $[a\text{,}b]=J=[0\text{,}T]$, we drop it from the notation. We also drop $\mathrm{\Omega }$; thus, we write $L^{\infty }(W_{\infty }^1)$ for $L^{\infty }(0\text{,}T\text{;}{W}_{\infty }^1(\mathrm{\Omega }))$.

If $v=(v_1\text{,}{v}_2)$ is a vector function, we say that $v\in X$ if $v_1\in X$ and $v_2\in X$. We also define the special vector-function spaces and norms

$$H^m(\mathrm{div}\text{;}\mathrm{\Omega })=\{v:v_1\text{,}{v}_2\text{,}\nabla ·v\in H^m(\mathrm{\Omega })\}\text{,}$$

$$\Vert v{\Vert }_{H^m(\mathtt{div})}={[\Vert v_1{\Vert }_m^2+\Vert v_2{\Vert }_m^2+\Vert \nabla ·v{\Vert }_m^2]}^{1/2}\text{,}m⩾0\text{,}$$

$$H(\mathrm{div}\text{;}\mathrm{\Omega })=H^0(\mathrm{div}\text{;}\mathrm{\Omega })\text{.}$$

Our assumptions on the regularity of the solution of (1.1) are denoted collectively by

$$(\mathtt{R})\left\{\begin{array}{c}e\text{,}p\in L^{\infty }(H^2)\cap H^1(H^1)\cap L^{\infty }(W_{\infty }^1)\cap H^2(L^2)\text{,}\hfill \\ \psi \in L^{\infty }(H^1)\text{,}\hfill \\ u\in L^{\infty }(H^1(\mathtt{div}))\cap L^{\infty }(W_{\infty }^1)\cap W_{\infty }^1(L^{\infty })\cap H^2(L^2)\text{.}\hfill \end{array}\right.$$

The hypotheses in ($\mathtt{R}$) enforce tacit conditions on the coefficients in (1.1). We will make explicit use of :

$$(\mathtt{C})0⩽D_{\ast }⩽D_s(x)⩽D^{\ast }\text{,}|{\mu }_s(x)|+|\nabla {\mu }_s(x)|⩽{\mu }^{\ast }\text{,}s=e\text{,}p\text{,}$$

for constants $D_{\ast }\text{,}{D}^{\ast }$ and ${\mu }^{\ast }$. These assumptions are reasonable in physics [9]. We also assume that $R(e\text{,}p)$ is Lipschitz continuous in $\epsilon$-neighborhood of the solutions. Exactly, when $|{\epsilon }_i|⩽{\epsilon }_0$ ($1⩽i⩽2$), we have

$$|R(e(x\text{,}t)\text{,}p(x\text{,}t))-R(e(x\text{,}t)+{\epsilon }_1\text{,}p(x\text{,}t)+{\epsilon }_2)|⩽K\{|{\epsilon }_1|+|{\epsilon }_2|\}\text{,}(x\text{,}t)\in \mathrm{\Omega }\times J\text{.}$$

An outline of the paper follows. In Section 2, we formulate the characteristics-mixed finite element method for the concentrations and mixed finite element method for the electric field intensity and electrostatic potential. Then we present the approximation schemes that combine the characteristics-mixed finite element, mixed finite element method and a post-processing step. Properties of finite element spaces and projections are listed in Section 3. In Section 4, we derive the optimal-order error estimates. Finally in Section 5, we give a numerical experiment to validate the theoretical analysis.

Throughout, the symbols K and $\epsilon$ will denote, respectively, a generic constant and a generic small positive constant, which may have different values in different places.