International audienceWe establish uniform L ∞ bounds for approximate solutions of the drift-diffusion system for electrons and holes in semiconductor devices, computed with the Schar-fetter-Gummel finite-volume scheme. The proof is based on a Moser iteration technique adapted to the discrete case

Bessemoulin-Chatard, Marianne

Chainais-Hillairet, Claire

Jüngel, Ansgar

INRIA a CCSD electronic archive server

HAL Id: hal-01472643https://hal.archives-ouvertes.fr/hal-01472643Submitted on 21 Feb 2017HAL is a multi-disciplinary open accessarchive for the deposit and dissemination of sci-entific research documents, whether they are pub-lished or not. The documents may come fromteaching and research institutions in France orabroad, or from public or private research centers.L’archive ouverte pluridisciplinaire HAL, estdestinée au dépôt et à la diffusion de documentsscientifiques de niveau recherche, publiés ou non,émanant des établissements d’enseignement et derecherche français ou étrangers, des laboratoirespublics ou privés.Uniform L ∞ estimates for approximate solutions of thebipolar drift-diffusion systemMarianne Bessemoulin-Chatard, Claire Chainais-Hillairet, Ansgar JüngelTo cite this version:Marianne Bessemoulin-Chatard, Claire Chainais-Hillairet, Ansgar Jüngel. Uniform L ∞ estimatesfor approximate solutions of the bipolar drift-diffusion system. FVCA 8, Jun 2017, Lille, France.￿hal-01472643￿Uniform L∞ estimates for approximate solutionsof the bipolar drift-diffusion systemM. Bessemoulin-Chatard∗ C. Chainais-Hillairet† A. Jüngel‡February 21, 2017AbstractWe establish uniform L∞ bounds for approximate solutions of the drift-diffusionsystem for electrons and holes in semiconductor devices, computed with the Schar-fetter-Gummel finite-volume scheme. The proof is based on a Moser iterationtechnique adapted to the discrete case.Keywords: volume scheme, drift-diffusion, Moser iterationsMSC (2010): 65M08, 35B40.1 IntroductionWe consider the Van Roosbroeck’s bipolar drift-diffusion system on Ω× (0,T ), whereΩ is a domain of Rd (d = 2, 3):∂tN + div(−∇N +N∇Ψ) =−R(N,P), (1a)∂tP+ div(−∇P−P∇Ψ) =−R(N,P), (1b)−λ 2∆Ψ = P−N +C. (1c)The unknowns are the electron density N, the hole density P and the electrostatic po-tential Ψ. The doping profile C(x) is given and λ is the scaled Debye length. Thissystem is supplemented with initial densities N0, P0, Dirichlet boundary conditionson ΓD (ND, PD, ΨD) and homogeneous Neumann boundary conditions on ΓN (with∂Ω = ΓD ∪ΓN , ΓD ∩ΓN = /0, and m(ΓD) > 0). The Dirichlet boundary conditionsare describing Ohmic contacts, while homogeneous Neumann boundary conditions arefor the insulated boundary segments. Dirichlet boundary conditions for (1) may de-pend on time, but we assume time-independent data to simplify. The recombination-generation rate is written under the following form which includes Shockley–Read–Hall and Auger terms:R(N,P) = R0(N,P)(NP− 1). (2)In what follows, we consider the following (standard) assumptions:∗Univ. Nantes, CNRS, UMR 6629 – Laboratoire Jean Leray, F-44000 Nantes, France,Marianne.Bessemoulin@univ-nantes.fr†Univ. Lille, CNRS, UMR 8524 – Laboratoire Paul Painlevé, F-59000 Lille, France,Claire.Chainais@math.univ-lille1.fr‡Institute for Analysis and Scientific Computing, TU Wien, 1040 Wien, Austria, juengel@tuwien.ac.at1(H1) C ∈ L∞(Ω),(H2) ND, PD ∈ L∞ ∩H1(Ω), ΨD ∈ H1(Ω),(H3) ND PD = 1,(H4) ∃M > 0 such that 0 ≤ N0, P0, ND, PD ≤ M a.e. on Ω,(H5) ∃R̄ > 0 such that 0 ≤ R0(N,P)≤ R̄(1+ |N|+ |P|) ∀N, P ∈ R.Hypothesis (H3) means that the boundary data are in thermal equilibrium. Existenceand uniqueness of weak solutions to system (1) have been proved in [9]. Nonnegativityof the densities and uniform-in-time upper bounds have also been shown in [9]. Theproof is based on an approach proposed by Alikakos [1], closely related to the Moseriteration technique [12].Let ∆t > 0 be the time step and let consider an admissible mesh of Ω. It is given bya family T of control volumes, a family E of edges (or faces) and a family of points(xK)K∈T which satisfy Definition 9.1 in [7]. In the set of edges E , we distinguish theset of interior edges Eint from the set of boundary edges Eext . We split Eext into Eext =E Dext ∪E Next , where E Dext and E Next is the set of Dirichlet and Neumann boundary edges,respectively. For a control volume K ∈ T , we denote by EK = Eint,K ∪E Dext,K ∪E Next,K .For all σ ∈ E , we define τσ = |σ |/dσ , where dσ = d(xK ,xL) for σ = K|L ∈ Eint , anddσ = d(xK ,σ) for σ ∈ Eext . We also need the following assumptions on the mesh:∃ξ > 0 such that d(xK ,σ)≥ ξ dσ , ∀K ∈ T , ∀σ ∈ EK , (3a)∃c0 > 0 such that τσ ≥ c0, ∀σ ∈ E . (3b)A finite volume discretization for (1) provides an approximate solution unT= (unK)K∈Tfor all n ≥ 0 and approximate boundary values uE D = (uσ )σ∈E Dext for u = N, P, Ψ. Forany vector uM = (uT ,uE D), we defineDK,σ u = uK,σ − uK, Dσ u = |DK,σ u|, ∀K ∈ T ,∀σ ∈ EK ,where uK,σ is either uL (σ = K|L), uσ (σ ∈ E DK,ext ) or uK (σ ∈ E NK,ext ). We also definethe discrete H1-seminorm | · |1,M by|uM |21,M = ∑σ∈Eτσ (Dσ u)2, ∀uM = (uT ,uE D).We define the initial conditions N0K , P0K as the mean values of N0 and P0 over K ∈ T .The boundary conditions are also approximated by taking the mean values of ND, PDand ΨD over each Dirichlet boundary edge σ ∈ E Dext .We are now in the position to define the scheme for (1), based on a backward Eulerin time discretization. For all K ∈ T and n ≥ 0,|K|Nn+1K −NnK∆t+ ∑σ∈EKFn+1K,σ =−|K|R(Nn+1K ,Pn+1K ), (4a)|K|Pn+1K −PnK∆t+ ∑σ∈EKGn+1K,σ =−|K|R(Nn+1K ,Pn+1K ), (4b)−λ 2 ∑σ∈EKτσ DK,σ Ψn = |K|(PnK −NnK +CK), (4c)2where F n+1K,σ and Gn+1K,σ are the Scharfetter-Gummel fluxesFn+1K,σ = τσ[B(−DK,σ Ψn+1)Nn+1K −B(DK,σ Ψn+1)Nn+1K,σ], (5a)Gn+1K,σ = τσ[B(DK,σ Ψn+1)Pn+1K −B(−DK,σ Ψn+1)Pn+1K,σ], (5b)and B is the Bernoulli function B(x) = x/(ex − 1) for x 6= 0, B(0) = 1.In [2], the existence of a solution to scheme (4)-(5) and the boundedness of theapproximate densities are shown, but the bounds depend on time and blow up whentime goes to infinity. The only case where the result is uniform in time is that of zerodoping profile. The purpose of this paper is to adapt the ideas developed in [9, 10]to the discrete framework to obtain uniform-in-time L∞ estimates for the approximatedensities obtained with scheme (4)–(5) for general doping profiles. Our main resultreads as follows.Theorem 1. Let (H1)–(H5) hold and let M = (T ,E ,P) be an admissible mesh of Ωsatisfying (3). Any solution (NnT,PnT,ΨnT)n≥0 to (4)–(5) satisfies∃κ > 0, ∀n ≥ 0, ‖NnT ‖L∞(Ω) ≤ κ and ‖PnT ‖L∞(Ω) ≤ κ . (6)The constant κ depends only on the initial and boundary data, C, λ , R̄, Ω and d, andon the constants ξ and c0 given in (3), but not on n.This theorem establishes a part of the assumptions needed to prove the exponentialdecay of approximate solutions given by scheme (4)–(5) towards an approximation ofthe thermal equilibrium [2, Theorem 3.1]. However, a uniform positive lower boundfor the densities is also required, which is not easy to prove, and its proof is an openproblem.The proof of (6) applies a Nash-Moser type iteration method based on Lr bounds [1,12]. Let us mention that this method has already been applied to a discrete setting in[8]. As we deal here with equations on a bounded domain, we have to take care aboutthe boundary conditions. Therefore, as in [11], we establish (6) for NM =(N−M)+ andPM = (P−M)+, where M is given in (H4), instead of N and P. The proof is detailed inSection 3. The uniform-in-time L1 bounds for the densities necessary to initialize theMoser iteration method are obtained thanks to an entropy-entropy production estimate,recalled in Section 2.2 Discrete entropy-entropy production inequalityThe thermal equilibrium is a steady state for which the electron and hole current den-sities and the recombination-generation term vanish. If (H3) is satisfied, there existsα ∈R such that the thermal equilibrium is defined byN∗ = eα+Ψ∗, P∗ = e−α−Ψ∗, (7a)−λ 2∆Ψ∗ = e−α−Ψ∗ − eα+Ψ∗ +C, (7b)Ψ∗ = ΨD on ΓD, ∇Ψ∗ ·ν = 0 on ΓN . (7c)An approximation of the thermal equilibrium (N∗T,P∗T,Ψ∗T) is given by−λ 2 ∑σ∈EKτσ DK,σ Ψ∗ = |K|(e−α−Ψ∗K − eα+Ψ∗K +CK), ∀K ∈ T , (8)N∗K = eα+Ψ∗K , P∗K = e−α−Ψ∗K , ∀K ∈ T . (9)3Let H(x) = x logx− x+ 1. The discrete relative entropy is defined byEn =λ 22|ΨnM −Ψ∗M |21,M + ∑K∈T|K|(H(NnK)−H(N∗K)− logN∗K(NnK −N∗K)+H(PnK)−H(P∗K)− logP∗K(PnK −P∗K)). (10)We also define the discrete entropy production:In = ∑σ∈E ;K=Kστσ[min(NnK ,NnK,σ )(Dσ (log(Nn)−Ψn))2+ min(PnK ,PnK,σ )(Dσ (log(Pn)+Ψn))2]+ ∑K∈T|K|R0(NnK ,PnK)(NnKPnK − 1) log(NnKPnK), (11)We recall the discrete entropy-entropy production inequality proved in [5].Proposition 1. For all n ≥ 0,0 ≤ En+1 +∆tIn+1 ≤ En. (12)Summing (12) over n, we have En ≤ E0, which gives a uniform-in-time estimatefor En. Then, if M satisfies (3b), and since there exists C∗ such that |Ψ∗M|1,M ≤ C∗(see [4, Lemma 3.3]), we have Dσ Ψn+1 ≤ D, where D > 0 only depends on E0, λ , c0and C∗. The properties of the Bernoulli function ensure that∃γ ∈ (0,1], B(Dσ Ψn+1)≥ γ, ∀σ ∈ E ,∀n ≥ 0. (13)3 Proof of Theorem 1We setNnM,K = (NnK −M)+, PnM,K = (PnK −M)+, ∀K ∈ T , ∀n ≥ 0,and V nq = ∑K∈T|K|[(NnM,K)q +(PnM,K)q], ∀n ≥ 0, ∀q ≥ 1.We start by establishing the following result about the evolution of V nq+1.Proposition 2. Let q ≥ 1. There exist positive constants µ and ν only depending on‖C‖∞, λ , M, R̄ and γ ∈ (0,1] such that1∆t(V n+1q+1 −V nq+1)+4qq+ 1γ ∑σ∈E[(Dσ (Nn+1M )q+12 )2 +(Dσ (Pn+1M )q+12 )2]≤ µqV n+1q+1 +ν|Ω|. (14)Proof. Multiplying (4a) (resp. (4b)) by (Nn+1M,K )q (resp. (Pn+1M,K )q) , summing over K andadding the two equations, we obtain S1 + S2 = S3, where S1 contains the discrete timederivatives, S2 the numerical fluxes and S3 the recombination-generation term. Usingthe elementary identity (x−y)xq ≥ (xq+1 −yq+1)/(q+1) for all x,y ≥ 0 and q ≥ 1, wefind thatS1 ≥1q+ 11∆t(V n+1q+1 −V nq+1). (15)4By a discrete integration by parts on S2, combined with some properties of the Bernoullifunction, we haveS2 ≥4q(q+ 1)2 ∑σ∈Eτσ B(Dσ Ψn+1)[(Dσ (Nn+1M )q+12 )2 +(Dσ (Pn+1M )q+12 )2]− qq+ 1 ∑σ∈Eτσ DK,σ Ψn+1[DK,σ ((Nn+1M )q+1)−DK,σ ((Pn+1M )q+1))]−M ∑σ∈Eτσ DK,σ Ψn+1[DK,σ ((Nn+1M )q)−DK,σ ((Pn+1M )q))].We perform a discrete integration by parts of the two last sums on the right-hand side,use scheme (4c) and the monotonicity of the functions x 7→ ((x−M)+)q and x 7→ ((x−M)+)q+1. Combined with (13), this yieldsS2 ≥4q(q+ 1)2γ ∑σ∈Eτσ[(Dσ (Nn+1M )q+12 )2 +(Dσ (Pn+1M )q+12 )2]− qq+ 1‖C‖∞λ 2V n+1q+1 −M‖C‖∞λ 2V n+1q . (16)Thanks to (H5) and the nonnegativity of the approximate densities, we haveR0(Nn+1K ,Pn+1K )(1−Nn+1K Pn+1K )[(Nn+1M,K)q+(Pn+1M,K)q]≤ R̄(1+ 2M)[(Nn+1M,K )q +(Pn+1M,K )q]+ R̄[(Nn+1M,K )q+1 +(Pn+1M,K )q+1]+ R̄[(Nn+1M,K )qPn+1M,K +(Pn+1M,K )qNn+1M,K]. (17)Then, applying the Young’s inequality, we obtainV n+1q ≤1q+ 1(qV n+1q+1 +m(Ω)),∑K∈T|K|[(Nn+1M,K)qPn+1M,K +(Pn+1M,K )qNn+1M,K]≤V n+1q+1 .Combining this with (15), (16) and (17) finishes the proof.Now, our aim is to control the term V n+1q+1 appearing on the right-hand side of (14).The discrete Nash inequality [3, Corollary 4.5] reads for functions χT that vanish on apart of the boundary as(∑K∈T|K|χ2K)1+ 2d≤ C̃ξ(∑σ∈Eτσ (Dσ χ)2)(∑K∈T|K||χK |)4d,where ξ is given in (3a) and C̃ only depends on Ω and d. Thanks to Young’s inequality,it follows for ε > 0 that, up to a change of the value of C̃,∑K∈T|K|χ2K ≤C̃εd/2ξ d/2(∑K∈T|K||χK |)2+ ε(∑σ∈Eτσ (Dσ χ)2).5Applying this inequality to χ =(Nn+1M)q+12 and χ =(Pn+1M)q+12 , we haveV n+1q+1 ≤C̃(εξ )d/2(V n+1q+12)2+ ε ∑σ∈E[(Dσ (Nn+1M )q+12 )2 +(Dσ (Pn+1M )q+12 )2]. (18)Arguing similarly as in [6] and using the fact that γ ∈ (0,1], we can find A > 0 depend-ing only on µ and hence only on ‖C‖∞, λ , M and R̄ such thatγAq(µq+γAq)≤ 4γqq+ 1, ∀q ≥ 1.Therefore, multiplying (18) by µq+ ε(q) with ε(q) = γA/q and adding the resultingequation to (14), we infer thatV n+1q+1 −V nq+1∆t≤−ε(q)V n+1q+1 +ν|Ω|+C̃ε(q)d/2ξ d/2(µq+ ε(q))(V n+1q+12)2. (19)Let us now define W nk = Vn2kfor all n ≥ 0 and k ∈ N. The definitions of M andthe initial condition ensure that W 0k = 0 for all k ∈ N. Moreover, the discrete entropy-entropy production inequality (12) ensures that En ≤ E0 for all n ≥ 0 and applying theinequalities∀x,y > 0 x log xy− x+ y ≥ (√x−√y)2 ≥ x2− y,we deduce a uniform bound of W n0 for all n ≥ 0. With q = 2k −1 = ζk and εk = γA/ζk,we infer from (19) thatW n+1k −Wnk∆t≤−εkW n+1k +B(ζd/2k (ζk + εk)(Wn+1k−1 )2 + 1)(20)with B = γ−d/2 max{νm(Ω), C̃ξ d/2A−d/2, C̃ξ d/2A−d/2µ}. Therefore, if W nk−1 is boundedfor all n by E , we conclude from (20) thatW nk ≤Bεk(ζd/2k (ζk + εk)E2 + 1), ∀n ≥ 0.Set δk = Bζd/2k (ζk + εk)/εk. As ζd/2k (ζk + εk)≥ 1, it follows thatW nk ≤ δk(E2 + 1), ∀n ≥ 0. (21)We prove by induction (see [1, 11]) that for all k ≥ 0,W nk ≤ 2δk(2δk−1)2 · · · (2δ1)2k−1K2k , ∀n ≥ 0,where K = max(1,supn≥0Wn0 ). This is a direct consequence of (21), remarking thatwith E = 2δk(2δk−1)2 · · · (2δ1)2k−1K 2kfor all k ≥ 0, we have 1 ≤ E2 (thanks to thedefinition of K ).To conclude, we first remark that δk ≤ D2(2+d/2)k with D = B/A. Hencek−1∏j=0(2δk− j)2 j ≤ (2D)2k−1 ·2(2+d/2)∑k−1j=0(k− j)2 j ≤ (2D)2k ·2(2+d/2)·2k ∑∞ℓ=1 ℓ2−ℓ ,6and since ∑∞ℓ=1 ℓ2−ℓ = 2, we find that W nk ≤ (25+dDK )2k. Taking the power 1/2k ofW nk we obtain‖NnM‖L2k (Ω) ≤ 25+dDK , ‖PnM‖L2k (Ω) ≤ 25+dDK , ∀n ≥ 0, ∀k ∈N,and passing to the limit k → ∞ gives‖NnM‖L∞(Ω) ≤ 25+dDK , ‖PnM‖L∞(Ω) ≤ 25+dDK , ∀n ≥ 0.Acknowledgements: The authors have been partially supported by the bilateralFrench-Austrian Amadée-ÖAD project. M. B.-C. thanks the project ANR-14-CE25-0001 Achylles. C. C.-H. thanks the team Inria/Rapsodi, the ANR Moonrise and theLabex Cempi (ANR-11-LABX-0007-01) for their support. A.J. acknowledges partialsupport from the Austrian Science Fund (FWF), grants P22108, P24304, and W1245.References[1] N. D. Alikakos. Lp bounds of solutions of reaction-diffusion equations. Comm.Partial Differential Equations, 1979.[2] M. Bessemoulin-Chatard and C. Chainais-Hillairet. Exponential decay of a finitevolume scheme to the thermal equilibrium for drift–diffusion systems. JNUM,2016.[3] M. Bessemoulin-Chatard, C. Chainais-Hillairet, and F. Filbet. On discrete func-tional inequalities for some finite volume schemes. IMA Journal of NumericalAnalysis, 2015.[4] C. Chainais-Hillairet and F. Filbet. Asymptotic behavior of a finite volumescheme for the transient drift-diffusion model. IMA J. Numer. Anal., 2007.[5] M. Chatard. Asymptotic Behavior of the Scharfetter–Gummel Scheme for theDrift-Diffusion Model. In FVCA VI. Springer Berlin Heidelberg, 2011.[6] M. Di Francesco, K. Fellner, and P. A. Markowich. The entropy dissipationmethod for spatially inhomogeneous reaction–diffusion-type systems. Proc. R.Soc. Lond. Ser. A Math. Phys. Eng. Sci., 2008.[7] R. Eymard, T. Gallouët, and R. Herbin. Finite volume methods. In Handbook ofnumerical analysis. North-Holland, 2000.[8] A. Fiebach, A. Glitzky, and A. Linke. Uniform global bounds for solutions of animplicit Voronoi finite volume method for reaction-diffusion problems. Numer.Math., 2014.[9] H. Gajewski and K. Gröger. On the basic equations for carrier transport in semi-conductors. J. Math. Anal. Appl., 1986.[10] H. Gajewski and K. Gröger. Semiconductor equations for variable mobilitiesbased on Boltzmann statistics or Fermi-Dirac statistics. Math. Nachr., 1989.[11] R. Kowalczyk. Preventing blow-up in a chemotaxis model. J. Math. Anal. Appl.,2005.[12] J. Moser. A new proof of de Giorgi’s theorem concerning the regularity problemfor elliptic differential equations. Commun. Pur. Appl. Math., 1960.7

Uniform L ∞ estimates for approximate solutions of the bipolar drift-diffusion system

HAL-Univ-Nantes

https://hal.archives-ouvertes.fr/hal-01472643

Uniform L ∞ estimates for approximate solutions of the bipolar drift-diffusion system

Abstract

Similar works

Full text

Available Versions

INRIA a CCSD electronic archive server

HAL-Univ-Nantes