AN HYBRID FINITE ELEMENT METHOD FOR A QUASI-VARIATIONAL INEQUALITY MODELING A SEMICONDUCTOR

. A problem of determining the characteristics of a semiconductor can be reduced to the study of the quasi-variational inequality, (J. Abouchabaka, R. Aboula¨ıch, A. Nachaoui and A. Souissi, COMPEL 18 (1999) 143–164.) where the obstacle ℳ ( 𝑢 ) is the solution of an elliptic problem depending on 𝑢. We present here an hybrid finite element method for the computation of obstacle ℳ ( 𝑢 ) and we discuss some numerical aspects appearing in its approximation.


Introduction
The stationary distribution of the carrier motions inside a semiconductor device is described by the following system of nonlinear elliptic partial differential equations: where (1) is the Poisson equation, (2) is the continuity equations for electrons and (3) is the continuity equations for holes.The domain Ω is an open of R 2 corresponding to the interior of the semiconductor device, of lipchitzian boundary Ω with , Ω = Γ  ∪ Γ  .These equations are subject to the following boundary conditions: =   ,  =   ,  =   ; on Γ  , where   ,   ,   are given functions on Γ  and  denotes the external normal vector of Ω.In the van Roosbroeck model ( 1)-( 3), the potential , the electron density  and the hole density  are the unknowns; the doping function , the charge constant  and the permittivity  are given.In the expression of the current density   the quantities   ,   represent the mobility and the diffusivity of the electrons.This problem has been studied by several authors [19,20,30,31], especially in a fairly general framework, [13,15,16,29] and for more realistic boundary conditions by the author [13,16].Uniqueness is not possible in general, it is obtained under restrictive hypotheses [14,17,23].
Several authors are interested in the approximation of this system.For example, a finite element approximation is done in [22,24,28,34], mixed finite elements is in [32], and several approximations based on the use of finite volumes can be consulted in [4,6,7] and the references therein.A finite-difference approximation is used in [2] for the equilibrium electrostatic potential.
Concerning the convergence of approximations, we refer the reader to [24] where he will find a general framework for the study of the convergence of approximations whatever their types.An application of this study to finite elements can be found in [22].
Under physical assumptions, we show that the system (1)-(2) leads as a limit case to a quasi-variational equation [1,12,25], given in the case of a unipolar field-effect semiconductor, by with () = { ∈  1 (Ω),  =  on Γ  and  ≤ ℳ() a.e.dans Ω} where ℳ() is the solution of problem with   =   = 0 on the source,   =   =  + ≥ 0 on the drain,   =  − ≤ 0 and   =  + on the grid (see Fig. 1).This IQV model is a formulation of a free boundary problem where we have to find the unknown boundary separating two sets Ω  and Ω  , corresponding to the charge neutrality region where the potential  and the Fermi quasi-potential  = ℳ() are identical and the depletion region where  <  [1,21,25].
The existence and uniqueness of the solution to the problem (6)-( 7) are discussed in [1] where the authors propose a technique based on topological degree theory which extends the results of [20] to cases where the solution is not regular.The solution of the quasi-variational inequality being obtained as the limit of a series of variational inequalities where the two problems ( 6) and ( 7) are decoupled The approximation of (6) poses no problem, any reasonable approximation scheme works.The difficulty appears in the discretization of (7) where it is necessary to take into account the electric field which can be very important in certain regions of the semiconductor.We are interested here in a hybrid discretization, [5], whose main property is the preservation of electric current, which is a particularly important property in practice, and an exact computation of the integral of the exponential.For recent use of this type of approximation in different fields see [3,9,18,26,27,33] and the references therein.

Hybrid method for the obstacle approximation
Using the change of variable  =  −ℳ() , we obtain new expressions of the equations verified by the electric current  : The multiplication of the first equation of ( 9) by  ∈ [ 2 (Ω)] 2 and the second equation by  ∈  1 (Ω),  = 0 sur Γ  and integration on Ω gives rise to the following problem: find The hybrid discretization will be based on the approximation of (10).For this, consider a regular triangulation  ℎ of Ω and define the following finite-dimensional spaces: where P  ( ) is the set of polynomials in ,  of degree less than or equal to 1, restricted to .The current  and the function  are then approximated by  ℎ ∈  ℎ and  ℎ ∈    ℎ such that The choice of  ℎ and  0 ℎ results in the uniqueness of the solution of ( 11), (see [5]).On each triangle ,  ℎ and  ℎ are constant, the system ( 11) is then written The choice in the first equation of ( 12) of an element of  ℎ not zero on the triangle  and vanishing elsewhere results in: The introduction of the relation (13) in the second equation of (12) gives a system where the unique unknown is  ℎ where   is defined by   = ( ) The following proposition allows us to define an approximation of ℳ() Proposition 2.1.Suppose that the triangulation  ℎ is such that for all  ∈  ℎ the angles of  are ≤  2 , then min Proof.The matrix  associated with the system ( 14) is unreduced, symmetric with strictly dominant diagonal.These coefficients verify the following relation As the second member of ( 14) is zero, we can deduce (cf.Ciarlet [8]) that the principle of discrete maximum can be applied to  ℎ , that is: Since   () > 0 ∀ ∈ Γ  , we then define an approximation  ℎ de  by ℳ ℎ () =  ℎ {− ln  ℎ ()}.

Approximation of Q. V. I.
Consider now the convex  ℎ () defined by the finite element approximation of Q.V.I. is then given by Proposition 3.1.The approximate problem (15) admits a solution in    ℎ .
Proof.Let  ℎ be the solution of and let   be the operator who with  ℎ ∈  ℎ associates  ′ ℎ solution of the variational equation where Any fixed point of  1 is a solution of (15).As for the problem (1)-( 5) (see [23]), the fact that   () ∈  ℎ and an argument of the principle of the discrete maximum [8] implies that any fixed point  ℎ of   is such that  − ≤  ℎ ≤  + .Let  be the open defined by: where  is a constant much greater than max(| − |,  + ), it is clear that no fixed point of   belongs to  therefore deg[ −   , , 0] is defined and independent of , where deg[•, •, •] is the topological degree of Brouwer [10,11] and  is the identity of  ℎ .
The questions concerning the convergence of  ℎ towards  are natural and will be published in a future work.We will treat the case where the operator ℳ is approached by a classical finite element method as well as the case of the new hybrid finite element approach.

Numerical aspects
To solve (15) which is also a Q.V.I., we use an iterative scheme based on the operator  1 , specifically we define the following 0 ℎ being chosen as the solution of the variational inequality where  0 is the Fermi quasi-potential called the Shokley approximation.

Treatment of exponential nonlinearity
The computation of the Fermi quasi-potential  ℎ , which serves as an obstacle in the computation of the potential  ℎ , solution of the variational inequation (17), is reduced to solving the system of linear equations where  ℎ = ( ) ,  ℎ is the contribution of the triangle  to the discretization of the Laplacian,   ℎ is the set of triangles having the node  as vertex,  ℎ is the set of nodes which are vertices of a triangle  ∈   ℎ ,   ℎ =   ℎ ∩   ℎ , and  ℎ the value of  ℎ at node .In the form (19) the computation of the coefficients   leads to capacity overruns (underflow-overflow), especially for the computation of the diagonal terms.
To eliminate this problem we will describe a technique which does not require any additional storage and which is very efficient for sufficiently small ℎ steps.
Multiply equation ( 18) by   ℎ we get that we write in the form ∑︁ Typically, the coefficients  ′  keep the properties of   , i.e., the property of the dominant diagonal.However, the system of linear equation ( 20) is no longer symmetric.But we can still manipulate this system, by storing only the upper part of the matrix ( ′  ).Indeed, the non-zero subdiagonal terms can be obtained as follows: The coefficient  ′  ,  > ; is obtained from the coefficient  ′  ,  > , by the relation This means that the solution of the system of equation ( 20) does not require any additional storage than that of the upper part of the matrix (  ) and the vector [ ℎ ]  .System (20) therefore preserves the main characteristics of ( 18), with the following properties: Proposition 4.1.For ℎ sufficiently small, the calculation of all the coefficients  ′  is done without any capacity overrun (neither underflow nor overflow).
Proof.The proof of the proposition is obtained from the relation Since  ℎ ∈ P 1 ( ), we can calculate in an exact way the integral  appearing in the previous relation.For simplicity we assume that  = ∑︁ ∈ ( ) where  ( ) is the set of vertices of the triangle .Since  ∈  ( ), ∀ ∈   ℎ , there is an  0 ∈  () such as  ( ℎ − ℎ 0 ) = 1, which eliminates any division by zero in the calculation of  ℎ .Furthermore,  ℎ ∈ P 1 ( ) leads to It is clear that the sign of ( ℎ −  ℎ ) can be negative or positive.However, for ℎ small enough we can ensure that − < ( ℎ −  ℎ ) <  où  is the lower limit value acceptable by the machine for   , this guarantees that the capacity is not exceeded in the calculation  ( ℎ − ℎ ) .

Computation of Drain current
The following proposition expresses the weak continuity of the electric current.
Proposition 4.2.Let ( ℎ ,  ℎ ) be the solution of the system (11), let  be a node of  ℎ ,  / ∈ Γ  and let   ℎ be the set of sides having  as vertex.Then ∑︁ Proof.Integration by parts, the second equation of the system (12) leads to the equation For the basis element   associated with the node  we obtain ∑︁ In semiconductor device technology the most important feature is the computation of the current flowing along a given contact.For FETs, it is the current of Drain   defined by where Γ  is the Drain contact.To computation of   by this formula requires the computation of the electric field and a numerical integration.In the following we describe how we can get around this problem by using a method that does not require any additional calculation.
Let  ℎ be the function of  ℎ defined by where   is the base function associated with node number  and   is the set of nodes belonging to the Drain.By construction  ℎ = 1 on Γ  and it is identically zero on the remainder of Ω, hence Since  has zero divergence, we get  We approach the current  by  ℎ , the current of Drain   is then approached by where  = (  ) is the matrix associated with the system (14).
Remark 4.3.The calculation of  ℎ by the formula ( 21) is reduced to the product of the matrix  by the vector [ ℎ ], whose elements are already computed and stored.

Numerical results
Numerical simulation is made using the following value of the main parameters:    = 1,  = 8 microns and  = 2 microns.
We set the gate voltage to 0 V and we apply a positive voltage to the drain contact thus the current flow through the channel and we can calculate   , the current of the Drain contact.
We increased the drain voltage, the current   increases up to the value 2.8 mA, then it becomes stationary.The same remark when   = −1 V the limit value of   approaches 1.2 mA.Then, the limit value only reaches 0.5 mA when   = −1.5 V.The current therefore becomes increasingly weak when the absolute value of   becomes increasingly large.This is consistent with the Figures 2-4 which show that the depth of the depletion region increases, when the value of   becomes smaller and smaller, until this zone extends through the active channel which stops the current movement in the channel.This gate-source potential phenomenon is known in physics as "Pinch-off voltage".
We set the drain voltage to +1 V and we vary the Gate voltage from 0 V to −1.8 V.The depth of the depletion region becomes higher when we increase the gate voltage.This continues until the boundary of the depletion  We observe from Figure 5 that the negative equipotentials do not cross the free boundary.They all remain confined inside line 7, equipotential of value equal to 0 V.The positive equipotentials go around line 7, cross the channel and therefore allow current to flow.
In Figure 6, where only the positive equipotentials have been represented, it can be seen that there are no more lines around the line number 1, where the potential is of value 0 V, which is no longer a closed curve.This shows that the channel is closed and there is no more current flow.

Conclusion
From these numerical results we have observed, (see the following figures), that -The depletion region is localized around the gate contact.
-As is known from physical experience, when applying a positive potential difference between source and drain contacts, an electric current flows from one end to an other.The gate contact allows us to control the drain current by modifying the thickness of the depletion zone: a sufficiently negative gate-source voltage cuts the conducting channel and thereby switches off the drain current.We can conclude that for the semiconductor device, used in the numerical experiments, the simplified model described in section introduction does not change the behavior of the full system.Also, we can conclude that the proposed hybrid finite element method provides a good alternative for solving free boundary problems arising in the depletion approximation model.This method combined with the scaling technique makes it possible to avoid capacity overruns and leads to a correct calculation of the current.