Boundary layer analysis in the semiclassical limit of a quantum drift diffusion model

Quantum Drift Diffusion Model

Shen Bian^∗, Li Chen^†, and Michael Dreher^‡ March 10, 2011

Abstract

We study a singularly perturbed elliptic second order system in one space variable as it appears in a stationary quantum drift diffusion model of a semiconductor. We prove the existence of solutions and their uniqueness as minimizers of a certain functional and determine rigorously the principal part of an asymptotic expansion of a boundary layer of those solutions. We prove analytical estimates of the remainder terms of this asymptotic expansion, and confirm by means of numerical simulations that these remainder estimates are sharp.

1 Introduction and Known Results

The distribution of charged particles in a semiconductor can be described by various systems of partial differential equations, for instance the drift–diffusion equations. Typically, the relevant classes of particles are the movable electrons in the conduction band (which is an energy band at a higher level) and the so-called holes (which are vacant positions in the next lower energy band, with positive charge). Additionally, charged ions may exist at fixed positions in the semiconductor crystal. For smaller devices, it may be necessary to consider also quantum mechanical effects. Then, in the stationary case, thequantum drift diffusion model reads (after scaling)

F =V +h_n(n)−ε²4√n

√n , (1.1)

G=−V +h_p(p)−ξε²4√p

√p , (1.2)

div (µ_nn∇F) =R₀(n, p)R₁(F, G), (1.3)

div (−µ_pp∇G) =−R₀(n, p)R₁(F, G), (1.4)

−λ²4V =n−p−C(x), (1.5)

∗Department of Mathematical Sciences, Tsinghua University, Beijing, 100084, P.R. China, bs09@mails.tsinghua.edu.cn

†Department of Mathematical Sciences, Tsinghua University, Beijing, 100084, P.R. China, lchen@math.tsinghua.edu.cn

‡Corresponding author. Fachbereich Mathematik und Statistik, Universit¨at Konstanz, 78457 Konstanz, Germany, michael.dreher@uni-konstanz.de

1

Konstanzer Online-Publikations-System (KOPS) URL: http://nbn-resolving.de/urn:nbn:de:bsz:352-137644

for the unknowns n, p, F, G, V. The functions n and p describe the densities of electrons and holes, and F, G are known as quantum quasi Fermi levels. Finally, V is the electric potential as generated from the charged particles via the Poisson equation (1.5). HereCcharacterizes the known density of positive ions. The functionsh_n, h_p are called the enthalpy functions of the electrons and holes; typically they are of the form h(s) =Tlnswith some positive temperature constant T. The positive parameter ε is proportional to the Planck constant ~ and describes quantum effects, and the positive constant ξ is related to the quotient of the effective masses of the electrons and the holes. Next, the functions R₀, R₁ are known expressions, which model generation–recombination effects, they satisfy certain monotonicity assumptions, and an example is

R₀(n, p)R₁(F, G) = 1

a₀+a₁|n|+a₂|p| exp(F +G)−δ² ,

with δ >0 chosen in such a way that R₀R₁ ≡0 in thermal equilibrium. The constants µ_n, µ_p are related to the mobilities of the particles, and λis known as Debye length.

This system has been studied extensively in [6], [8], [1]; and we also refer to [5], [4].

The above system can be complemented with certain boundary conditions, for instance Dirichlet boundary conditions of (n, p, V, F, G) on a boundary part Γ₊ (with n, p positive there), homoge- neous Neumann boundary conditions of (n, p, V, F, G) on a boundary part Γ_N, and (n, p, V) = (0,0, V_extern+V_equil) on a further boundary part Γ₀ (note that the elliptic equations (1.3) and (1.4) degenerate at points wheren=p= 0).

Several results were proved in [1] under appropriate assumptions: the full system has a solution (n, p, V, F, G) ∈ L^∞(Ω)∩H¹(Ω) ∩C(Ω) with non-negative n and p. And if F, G ∈ L^∞(Ω) are given, then a solution (n, p, V) to (1.1), (1.2), (1.5) exists which is uniquely determined by the condition that a certain functional shall be minimized. Finally, the semiclassical limit ε→ 0 has been performed, under the assumption that the above mentioned boundary part Γ₀ is empty.

It is one of the goals of this article to remove this restriction on Γ₀, for a sub-class of the systems studied in [1] which we describe now.

S

D G

Figure 1: A rough schematic sketch of a MOSFET, with fictitious boundary in the bulk of the material

Our studies are related to MOSFET¹ devices, whose structure is sketched in Figure 1. At one end of the device, there are contacts called source, gate, drain, of which the gate contact is insulated

1metal-oxide-semiconductor field-effect transistor

by means of an oxide, explaining the name of the device. Depending on an applied voltage V_GS between gate and source, the density of movable charge carriers changes, and this effect decides whether a current can flow from source to drain. If such movable particles are available between source and drain in abundance, we say that an inversion layer has formed. A good knowledge of this inversion layer, be it analytical or numerical, is of high importance to applications. In an opposite situation, when no movable particle are available in a certain region, we say that this region depleted of particles. Various asymptotic expansions for such layers have been proposed in [9], [10], [2]; we also refer to [3] for a model involving quantum effects. The modifications of [3] to the model (1.1)–(1.5) can be summarized as follows. The holes are supposed to be in thermal equilibrium, which implies G≡0. In many situations, the parameter ξ from (1.2) is small, which motivates us to neglect quantum effects for the holes, and then (1.2) turns into

0 =−V +h_p(p),

or equivalentlyp= exp(V /T_p). Moreover, generation–recombination events are ignored: R₀R₁ ≡0.

Next we assume that the domain Ω of the device is a square (0,1)×(0,1)⊂R², and we write the spatial variable as (x, y), withx running from the contacts into the bulk of the crystal. Concerning the electron densityn, we assume thermal equilibrium in direction of thexvariable, which makesF a function ofy∈(0,1) only, and this function is supposed to be known. Thequasi 1D approximation is supposed, which says that all functions depend only weakly on the variable y. Hence we will neglect the derivatives with respect to y from now on, and the system becomes

F =V(x) +T_nlnn(x)−ε²(p

n(x))_xx

pn(x) , 0< x≤1, (1.6)

−λ²V_xx(x) =n(x)−exp(V(x)/T_p)−C(x), 0≤x≤1. (1.7) In the sequel, the dependence of the functions on the y–variable is no longer mentioned.

The boundary conditions at the fictitious boundaryx= 1 in the bulk material are

n(1) =nB, V(1) =VB, (1.8)

and the constants nB∈R₊ and VB∈Rsatisfy the compatibility condition

F =VB+TnlnnB, (1.9)

which expresses that quantum effects do not appear there.

And the boundary conditions at the location of the gate (x= 0) are

n(0) = 0, V_x(0) =β(V(0)−V_GS), (1.10)

with given constantsβ ≥0 and V_GS∈R.

The vanishing of the particle density atx= 0 makes two terms in (1.6) singular, and the formation of an additional layer inside the inversion layer is to be expected: the quantum layer.

The objective of this paper is to study this quantum layer with the rigor of analysis. Our modelling follows [3], where an asymptotic expansion was conjectured, but the key improvement presented in

0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 0

0.5 1 1.5 2 2.5 3 3.5 4

n p

0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1

−0.2

−0.15

−0.1

−0.05 0 0.05

Figure 2: Electron and hole densities, and electric potential in case of an inversion layer this article is a precise discussion of the error terms. This enables us to perform the limit ε→ 0 rigorously, even in the presence of a zero boundary value of the electron density.

We sketch the results presented in this paper. First, we prove the existence of a solution to the system (1.6), (1.7), (1.8), (1.10), which is unique as a minimizer to a certain functional. Our approach builds upon [1]; however, the additional nonlinearity in our equation (1.7) which is not present in [1] requires several changes.

Second, we prove rigorously that the electric potential V_ε converges to the corresponding solution V_∗ of the classical model, with an analytic estimate of the error V_ε−V_∗. We also determine the profile of the quantum layer of the electron density, again with several analytical error estimates.

And our third result are numerical simulations which confirm that our analytical error estimates are sharp.

2 Main Results

Our first result is about the existence of solutions (n, V), which are constructed in such a way that

%:=√

nis the unique non-negative minimizer of a certain functional, see Remark 3.7.

Theorem 2.1. Let us be given positive constants T_n,T_p, ε, λ, and real constants F,n_B, V_GS, and V_B, and a non-negative constant β. Suppose that C =C(x) is continuous and real-valued.

Then the boundary-value problem (1.6), (1.7), (1.8), (1.10) possesses a solution (n, V)∈C²([0,1]) with n(x) > 0 for 0 < x ≤ 1, such that % := √

n is the unique non-negative minimizer to F from (3.8).

There is a constant C0, independent of ε, with ε²

Z 1 0

(√

n)⁰(x)2

dx+ Z 1

0

V⁰(x)2

dx≤C₀, (2.1)

for this solution (n, V).

Our second result will describe the shape ofnandV, in particular for small values ofε. The precise formulation requires some preparations.

By variational methods, we will show in Proposition 3.6 that unique functions V_∗, n_∗ ∈ C²([0,1]) exist which solve







−λ²V_∗⁰⁰(x) =n_∗(x)−exp(V_∗(x)/T_p)−C(x), 0≤x≤1, n_∗(x) := exp((F −V_∗(x))/T_n), 0≤x≤1,

V_∗⁰(0) =β(V_∗(0)−V_GS), V_∗(1) =V_B.

(2.2)

Let Z=Z(y) be the unique function from C²([0,∞)) that solves Z⁰⁰(y) =Z(y) lnZ(y), 0< y <∞, Z(0) = 0, lim

y→+∞Z(y) = 1. (2.3)

This function approaches the value 1 exponentially for y → +∞, and all its derivatives decay exponentially, as will be shown in Lemma A.1.

Theorem 2.2. Assume (1.9). Then there is a positive ε₀, such that, for 0< ε≤ε₀, we have n−n₍₀₎

L²((0,1))+

V −V₍₀₎

W2¹((0,1)) ≤Cε,

where the zero-th order approximations (n₍₀₎, V₍₀₎) of (n, V) are defined as follows:

n₍₀₎(x) :=n_∗(x)Z²p 2T_n· x

ε

, V₍₀₎(x) :=V_∗(x).

Theorem 2.3. The approximations n₍₀₎ satisfy the uniform estimates n−n₍₀₎

L^∞((0,1)) ≤Cε^3/4, 0< ε≤ε₀, and in particular, we have

n(x)−n₍₀₎(x)

≤Cε^3/4· x

ε, 0≤x≤4ε. (2.4)

We follow the standard notation conventions. In particular, C denotes a generic positive constant that is independent of the unknown functions and may change its value from line to line.

3 Existence of Solutions

In this section, we demonstrate Theorem 2.1, whose proof is split into several parts.

First we consider the boundary value problem







−λ²U⁰⁰(x) =q(x)−p(x, U(x)), 0≤x≤1, U⁰(0) =βU(0), β ≥0,

U(1) = 0,

(3.1)

with a given function q∈C([0,1]), under the following assumptions on the nonlinearity p:











p∈C¹([0,1]×R),

p_U(x, U)>0 ∀(x, U)∈[0,1]×R, P(x, U) :=

Z _U

u=0

p(x, u) du≥ −C ∀(x, U)∈[0,1]×R.

(3.2)

We start our considerations with a simple observation.

Lemma 3.1. Solutions U ∈C²([0,1]) to (3.1) are unique, and they satisfy an estimate kUk_W2¹_((0,1)) ≤C(1 +kqk_L²_((0,1))),

with a constant C depending only onkp(·,0)kL²((0,1)) and λ >0.

Proof. The uniqueness is proved by usual monotonicity arguments. And if U ∈ C²([0,1]) is a solution, then

λ² Z 1

0

(U⁰(x))²dx+λ²β(U(0))²

= Z 1

0

(q(x)−p(x,0))U(x) dx− Z 1

0

(p(x, U(x))−p(x,0))U(x) dx, and then the monotonicity of p andβ ≥0 imply

λ² U⁰

2

L²((0,1)) ≤ kq−p(·,0)kL²((0,1))· kUkL²((0,1)).

It remains to exploit Poincare’s inequality, which is possible since U(1) = 0.

To prove the existence of a solution to (3.1), we choose the variational ground space X_U :=

U ∈W₂¹((0,1)) : U(1) = 0 .

Lemma 3.2. For a fixed functionq ∈C([0,1]), define a functional F⁽⁰⁾(U) := λ²

2 Z 1

0

(U⁰(x))²dx+λ²

2 β(U(0))²+ Z 1

0 −q(x)U(x) +P(x, U(x)) dx. (3.3) This functional possesses a unique minimizer U₀ ∈X_U, and this minimizer is a classical solution to the boundary value problem (3.1).

Proof. By Poincare’s inequality and β≥0, the functional F⁽⁰⁾ is coercive:

F⁽⁰⁾(U)≥ λ²

4 kUk²W2¹ −C, ∀U ∈X_U,

and then the existence of a minimizer U₀ ∈ X_U follows by standard arguments, and U₀ ∈ X_U ⊂ C([0,1]) is bounded. Take ϕ∈C^∞([0,1]) withϕ(1) = 0, and considerF⁽⁰⁾(U₀+δϕ):

F⁽⁰⁾(U₀+δϕ)

=F⁽⁰⁾(U₀)−δ Z 1

0

λ²U₀⁰⁰+q−p(x, U₀)

ϕdx+δλ²

−U₀⁰(0) +βU₀(0)

ϕ(0) +O(δ²), and therefore U₀ solves the differential equation −λ²U₀⁰⁰ =q−p(x, U₀) in the distributional sense, and it follows that U₀∈C²([0,1]), as well as−U₀⁰(0) +βU₀(0) = 0.

Next we discuss how this minimizer U₀ of F⁽⁰⁾ depends on q.

Lemma 3.3. Define a mapping Φ :C([0,1])→C_B²([0,1]), with C_B²([0,1]) as the vector space of all functions U fromC²([0,1]) with U⁰(0) =βU(0) andU(1) = 0, via the relation Φ{q}:=U0, and U0

is defined as the unique minimizer of the functional F⁽⁰⁾ from (3.3).

Then Φis a homeomorphism.

Proof. Clearly, Φ is bijective, and Φ⁻¹ is continuous. It remains to establish the continuity of Φ. To this end, let q₁ ∈C([0,1]) be given, and q₂ ∈C([0,1]) with kq₁−q₂k_L^∞ ≤1. Define U_k := Φ{q_k}. Then −λ²U_k⁰⁰+p(x, U_k) =q_k, and we quickly find

λ² Z 1

0

U₁⁰ −U₂⁰2

dx+λ²β

U1(0)−U2(0)2

+ Z 1

0

(p(x, U1)−p(x, U2))(U1−U2) dx

= Z ₁

0

(q₁−q₂)(U₁−U₂) dx,

which brings us, together with β ≥ 0, the monotonicity of p, and Poincare’s inequality, that kU₁−U₂kW2¹ ≤ Ckq₁−q₂kL², hence also kU₁−U₂kL^∞ ≤ Ckq₁−q₂kL². Next it follows that kp(·, U₁)−p(·, U₂)kL^∞ ≤Ckq₁−q₂kL², and then also kU₁⁰⁰−U₂⁰⁰kL^∞ ≤Ckq₁−q₂kL^∞.

For more information, we determine the Fr´echet derivative of Φ.

Lemma 3.4. Take q₀ ∈ C([0,1]). Then there is a constant C such that for all q ∈ C([0,1]) with kq−q₀k_C([0,1])≤1 we have the expansion

Φ{q}= Φ{q₀}+W +R,

with W ∈C²([0,1]) as the unique solution to

−λ²W⁰⁰+p_U(x,Φ{q₀})W =q−q₀, W⁰(0) =βW(0), W(1) = 0, satisfying kWkC²([0,1])≤Ckq−q₀kC([0,1]), and kRkC²([0,1])≤Ckq−q₀k²C([0,1]).

Proof. From Lemma 3.3 we know already that kΦ{q} −Φ{q₀}kC²([0,1]) ≤ Ckq−q₀kC([0,1]). We defineR := Φ{q} −Φ{q₀} −W and have

−λ²R⁰⁰=−λ²(Φ{q}⁰⁰−Φ{q₀}⁰⁰−W⁰⁰)

= (q−p(x,Φ{q}))−(q₀−p(x,Φ{q₀})) +p_U(x,Φ{q₀})W −(q−q₀)

=−p_U(x,Φ{q₀})·(Φ{q} −Φ{q₀} −W) +O(kΦ{q} −Φ{q₀}k²_C([0,1]))

=−p_U(x,Φ{q₀})R+O(kq−q₀k²_C([0,1])).

Note that p_U(x,Φ{q₀}) is positive, hence we can conclude thatkRk_C²_([0,1])≤Ckq−q₀k²_C([0,1]). Lemma 3.5. Define a function K=K(x, U) :=U p(x, U)−P(x, U), and set

F⁽¹⁾(q) := λ² 2

Z 1 0

((Φ{q})⁰(x))²dx+λ²

2 β(Φ{q}(0))²+ Z 1

0

K(x,Φ{q}(x)) dx.

Then K(x, U)≥0 for all (x, U)∈[0,1]×R, and the Fr´echet derivative of F⁽¹⁾ is given via F⁽¹⁾(q) =F⁽¹⁾(q0) +

Z 1 0

Φ{q0} ·(q−q0) dx+O(kq−q0k²C([0,1])).

Here the remainder term is positive for q6=q₀, and F⁽¹⁾ is strictly convex.

Proof. Concerning the bound for K, we remark that P(x,0) = 0, p_U >0 and K(x, U) =

Z _U

u=0

up_U(x, u) du.

For the proof of the second claim, we putU₀= Φ{q₀}and U := Φ{q}. ThenU =U₀+W +R with W and R as given in Lemma 3.4, and it follows that

F⁽¹⁾(q) = λ² 2

Z 1 0

(U⁰(x))²dx+λ²

2 β(U(0))²+ Z 1

0

K(x, U(x)) dx

=F⁽¹⁾(q₀) +λ² Z 1

0

U₀⁰ ·(U −U₀)⁰dx+λ²βU₀(0)·(U(0)−U₀(0)) +λ²

2 Z ₁

0

(U⁰−U₀⁰)²dx+λ²

2 β(U(0)−U₀(0))²+ Z ₁

0

K(x, U)−K(x, U₀) dx

≥ F⁽¹⁾(q₀)−λ² Z ₁

0

U₀·(U −U₀)⁰⁰dx+ Z ₁

0

K(x, U)−K(x, U₀) dx

=F⁽¹⁾(q₀) + Z ₁

0

U₀·(q−q₀) dx +

Z ₁

0 −U₀·(p(x, U)−p(x, U₀)) +K(x, U)−K(x, U₀) dx.

Now we have, with some ˜U between U and U₀,

−U₀·(p(x, U)−p(x, U₀)) +K(x, U)−K(x, U₀)

=−U₀p(x, U) +U p(x, U)−P(x, U) +P(x, U₀)

=−

P(x, U) +p(x, U)·(U₀−U)

+P(x, U₀)

= 1

2P_{U U}(x,U˜)·(U₀−U)²≥0, because of P_{U U} =p_U >0.

Proposition 3.6. Suppose thatλand T_p are positive, V_GS,V_B are real, and n−C is a continuous given real-valued function on [0,1].

Then there is exactly one solution V ∈C²([0,1]) to (1.7) with the boundary conditions V(1) =V_B and V⁰(0) =β(V(0)−V_GS).

Moreover, the nonlinear boundary value problem (2.2) possesses exactly one solutionV_∗∈C²([0,1]).

Proof. To solve (1.7) with the mentioned boundary conditions, we put V =V_inh+U withV_inh as the unique solution to

(−λ²V_inh⁰⁰ (x) =−C(x), 0≤x≤1,

V_inh⁰ (0) =β(V_inh(0)−V_GS), V_inh(1) =V_B. (3.4) We then have U = Φ{n}, with Φ as defined in Lemma 3.3, and

p(x, U) := exp(V_inh(x)/T_p) exp(U/T_p). (3.5)

The uniqueness of V is proved as in Lemma 3.1. The result concerning (2.2) is proved likewise.

For n≥0 we may write n=%², and then we transform (1.6), (1.7), (1.8), (1.10) into ( %F =%(V_inh+ Φ{%²}) + 2T_n%ln%−ε²%⁰⁰, on [0,1],

%(0) = 0, %(1) =%_B :=√n_B, (3.6)

withV_inh given by (3.4).

Seth(s) =T_nlnsfors >0 and H(s) :=

Z s σ=1

h(σ) dσ =T_n(slns−s+ 1). (3.7)

Our intuition is to look for% as a non-negative minimizer to the functional F(%) :=

Z ₁

0

ε²|%⁰|²+H(%²) +%²V_inh−%²F

dx+F⁽¹⁾(%²), (3.8)

over the set X_% :={%∈W₂¹((0,1)) :%(0) = 0, %(1) =%_B}.

Remark 3.7. For the convenience of the reader, we collect all information about how to construct F in one place: for functions % from X_% we shall discuss the functional

F(%) :=

Z ₁

0

ε²|%⁰|²+H(%²) +%²V_inh−%²F dx +λ²

2 Z 1

0

((Φ{%²})⁰)²dx+λ²

2 β(Φ{%²}(0))² +

Z 1 0

Φ{%²}(x)·p(x,Φ{%²}(x))−P(x,Φ{%²}(x)) dx,

with V_inh as the unique solution to (3.4),p=p(x, U) defined in (3.5), and P =P(x, U) from (3.2).

Finally, U = Φ{q} is defined as the unique solution to (3.1). For the definition of F, we take q :=%².

However, some problem occurs here. The pole ofh=h(s) ats= 0 and the boundary values of%at x= 0 make the functional F irregular, and then the Euler–Lagrange equations can not be derived.

To overcome this difficulty, we perform a regularization step: for a parameterγ ∈(0,1), we set

h_γ(s) :=







h(γ) : 0≤s≤γ, h(s) : γ ≤s≤γ⁻¹, h(γ⁻¹) : γ⁻¹ ≤s, and the we put H_γ(s) := R_s

σ=1h_γ(σ) dσ. The functional for which we wish to find a non-negative minimizer is

Fγ(%) :=

Z ₁

0

ε²|%⁰|²+H_γ(%²) +%²V_inh−%²F

dx+F⁽¹⁾(%²), (3.9)

where we restrict γ to the interval (0, γ₀), and γ₀ with 0 < γ₀ 1 is selected by the condition H_γ(s)−skV_inh−FkL^∞((0,1))≥sfors≥γ₀⁻¹.

Lemma 3.8. For functions %taking only non-negative values, the functionalsF andF^γ from (3.8) and (3.9) are strictly convex functionals of %², in the following sense: if %₁, %₂ ∈X_% with %_1,2 ≥0 and %₁6=%₂, and if0< t <1, then

F q

t%²₁+ (1−t)%²₂

< tF(%₁) + (1−t)F(%₂), Fγ

q

t%²₁+ (1−t)%²₂

< tFγ(%₁) + (1−t)Fγ(%₂).

Non-negative minimizers of F or F^γ are unique.

Proof. First, the strict convexity of the functional %² 7→ R₁

0 |∇%|²dxwas shown in [6]. Second, the scalar functions H and H_γ are weakly convex functions of their arguments. And third, we recall the strict convexity ofF⁽¹⁾ from Lemma 3.5.

Lemma 3.9. For0< γ < γ₀, the functionalFγ possesses a non-negative minimizer%_γ ∈C²([0,1])∩ X_%, and each such minimizer satisfies the uniform in γ estimate

εk%γkW2¹((0,1))+

Φ{%²_γ}

W2¹((0,1)) ≤C. (3.10)

Proof. We clearly have Fγ(%)≥ ε²

2 k%k²W2¹((0,1))−C, %∈X_%,

with someC ∈R₊depending only onkV_inh−FkL^∞((0,1)), but not onγ. Thenm_γ := inf{Fγ(%) :%∈ X_%} exists. Let (%₁, %₂, . . .) be a sequence in X_% with lim_j→∞Fγ(%_j) =m. Then this sequence is bounded inW₂¹((0,1)), hence we can assume strong convergence inC([0,1]), and weak convergence in W₂¹((0,1)) of this sequence to some limit %_γ ∈ X_%. Then F⁽¹⁾(%²_j) has the limit F⁽¹⁾(%²_γ) for j→ ∞, because Φ : C([0,1])→C_B²([0,1]) is continuous, see Lemma 3.3. By classical arguments [7], we conclude that Fγ is weakly sequentially lower semi-continuous on X_%, and consequently Fγ has a minimizer on X%. If %γ ∈ X% is such a minimizer, then |%γ| also belongs to X%, and it has the same value of Fγ. Independent of γ estimates of such non-negative minimizers %_γ can be found by choosing a function%^∗∈X_% that is non-negative and bounded from above. ThenFγ(%_γ)≤ Fγ(%^∗), and the right-hand side is bounded from above independently ofγ becauseH_γ is uniformly bounded from below. This gives us the uniform estimates (3.10).

Lemma 3.10. Let %_γ ∈C²([0,1])∩X_% be a non-negative miminizer of Fγ, and V_γ :=V_inh+ Φ{%²_γ}. Then %_γ satisfies

−ε²%⁰⁰_γ + (h_γ(%²_γ) +V_γ−F)%_γ = 0. (3.11)

Proof. By Lemma 3.5, (3.11) is just the Euler–Langrange equation of the functional Fγ.

Proof of Theorem 2.1. By the uniform bound ofk%_γk_W2¹_((0,1))and the compactness of the embedding W₂¹((0,1)) ⊂ C([0,1]), we can assume that we have a sequence (%_γ)_γ→0 of non-negative solutions to (3.11) that converges in C([0,1]) to a non-negative limit %. Since the function s 7→ sh(s²) is continuous on [0,∞), we deduce that

(h_γ(%²_γ)−V_γ+F)%_γ →(h(%²)−V +F)%, V :=V_inh+ Φ{%²},

with convergence in C([0,1]), for γ →0. Then the limit% is a non-negative distributional solution to (3.6), and by elliptic regularity, %∈C²([0,1]), andε²

%⁰⁰_γ−%⁰⁰

L^∞((0,1)) →0.

To demonstrate the positivity of %(x) at all x∈(0,1), we recall that lim_s→+0h(s) =−∞. Assume that %(x^∗) = 0 at an interior point x^∗ ∈(0,1). By the continuity of s7→ √

sh(s), the function x7→%(x)(h(%²(x)) +V(x)−F)

takes only non-positive values near x^∗. Therefore %⁰⁰ ≤0 in an open neighbourhood Ω ⊂(0,1) of x^∗. On the other hand, %(x) ≥ 0 in Ω, which is only possible if % ≡ 0 in Ω. Therefore, the set of zeroes of % is an open subset of (0,1). On the other hand, it is a closed subset of (0,1) in the

relative topology because % is continuous. Hence the set of zeroes of % in (0,1) is either empty or all of (0,1). The latter is impossible because % > 0 at x = 1. Therefore % > 0 in (0,1), and from F(σ) = lim_γ→0Fγ(σ), with F from (3.8) and any functionσ ∈X_%, we learn that % is indeed a minimizer of F. It remains to show (2.1). Choose %^∗ ∈ X_%, for instance %^∗(x) = %_Bx. Then F(%)≤ F(%^∗), which is bounded independent ofε∈(0,1]. Observe that there is a constantC∈R₊ such that

Z ₁

0

H(ψ²) +ψ²(V_inh−F) dx≥ −C, F⁽¹⁾(ψ²)≥ −C, for all ψ∈C([0,1]), from which (2.1) quickly follows.

4 Asymptotics of the Solutions

Now we begin to demonstrate Theorem 2.2, whose proof will span the whole Section 4. From now on, let (n_ε, V_ε) be the solution to (1.6), (1.7), (1.8), (1.10), as constructed in Theorem 2.1. We introduce the notation %_ε:=√n_ε. Then the pair (%_ε, V_ε) solves

( ε²%⁰⁰_ε =g(%_ε) +%_ε·(V_ε−F),

−λ²V_ε⁰⁰=%²_ε−p(V_ε)−C(x), (4.1)

withg(s) := 2T_nslnsand p(s) := exp(s/T_p). Moreover, we have the boundary conditions (%_ε(0) = 0, %_ε(1) =%_B :=√n_B,

V_ε⁰(0) =β(Vε(0)−VGS), Vε(1) =VB.

4.1 Properties of the Solutions

Lemma 4.1. There is a constantC1, independent ofε, such that kV_εkC²([0,1])+kn_εkC([0,1])≤C₁.

Proof. By the embedding W₂¹((0,1)) ⊂ C([0,1]), (2.1), and V_ε(1) = V_B, we find kV_εkC([0,1]) ≤ C.

Suppose that %_ε takes a local maximum at an inner point x^∗ ∈ (0,1). Then %⁰⁰_ε(x^∗) ≤ 0, hence h(%²_ε(x^∗)) +V_ε(x^∗)−F ≤0. We then have

%²_ε(x^∗)≤exp((F −V_ε(x^∗))/T_n),

and then also kn_εkC([0,1])≤C. Then (1.7) gives us the remaining bound for kV_εkC²([0,1]). The next result gives us a first information on the graph of n, depending on ε.

Lemma 4.2. There is a positive constant g_∗ such that

%_ε(x)≥g_∗, 2ε≤x≤1,

%⁰_ε(x)≥ g_∗

2ε, 0≤x≤4ε.

Proof. The function %_ε solves (4.1). Defineg_∗ as g_∗ = 1

3min %B

2 ,exp 1

2T_n(−1−C1)

,

withC₁ from Lemma 4.1 as an upper estimate of kV_ε−FkC([0,1]). Then it follows that whenever 0< %_ε(˜x)≤3g_∗, then

2T_nln%_ε(˜x) +V_ε(˜x)−F ≤ −1,

and consequently %⁰⁰_ε(˜x) < 0. Then %⁰_ε(˜x) ≤ 0 is impossible because this would imply %⁰_ε(x) < 0 for all x ∈ (˜x,1], contradicting %_ε(1) = %_B > 3g_∗. Hence there is a uniquely determined number x3∈(0,1) with

%_ε(x)

(<3g_∗ : 0≤x < x3,

>3g_∗ : x₃ < x≤1.

Moreover, on the interval [0, x₃], %_ε is strictly increasing. Define x₁, x₂ ∈ (0, x₃) uniquely by

%_ε(x₁) =g_∗ and%_ε(x₂) = 2g_∗. Then

%⁰⁰_ε(x)≤ −g_∗

ε², x∈[x₁, x₃].

Now we make use of a simple fact: if|ψ⁰⁰(x)| ≥d₀ on [a, b] for some ψ∈C²([a, b]), then max_[a,b]ψ− min_[a,b]ψ ≥d0(b−a)²/4. Hence we conclude that x3 −x2 ≤ 2ε and x2−x1 ≤2ε. Since %ε(x2,3) are independent of ε, and because of %⁰⁰_ε <0 on (0, x₃], we find that

%⁰_ε(x)≥ g_∗

2ε, x∈[0, x₂].

From the definition ofx₁ and %_ε(0) = 0 we then getx₁ ≤2ε, hencex₂ ≤4ε as desired.

The next result is a first step in showing that the quantum term ε²%⁰⁰_ε is of less relevance for x≥ O(√

ε).

Lemma 4.3. Assume (1.9). Then there is a constant C, independent ofε∈(0,1/8], such that k%ε(2Tnln%ε+Vε−F)kL^∞((2√

ε,1)) ≤Cε^1/4, (4.2)

k%_ε(2T_nln%_ε+V_ε−F)kL²((2√

ε,1)) ≤Cε^3/4, (4.3)

k%_ε(2T_nln%_ε+V_ε−F)kL²((4ε,1)) ≤Cε^1/2, (4.4)

%²_ε−exp((F −V_ε)/T_n)

L²((4ε,1)) ≤Cε^1/2, (4.5)

%⁰_ε L^∞((2√

ε,1)) ≤Cε^−3/4. (4.6)

Proof. Take a functionχ∈C^∞([0,1]) with 0≤χ≤1, χ⁰≥0,

χ(x) =

(3x : 0≤x≤1/4, 1 : 1/3≤x≤1,

and setχ_ε(x) =χ(x−2ε) forx∈[2ε,1]. Then we conclude, using (2T_nln%_ε+V_ε−F)|x=1= 0, that Z ₁

2ε

ε²χ_ε

%_ε (%⁰⁰_ε)²dx= Z ₁

2ε

χ_ε%⁰⁰_ε(2T_nln%_ε+V_ε−F) dx

=−2T_n Z ₁

2ε

χ_ε

%_ε(%⁰_ε)²dx− Z ₁

2ε

χ_ε%⁰_εV_ε⁰dx− Z ₁

2ε

χ⁰_ε%⁰_ε(2T_nln%_ε+V_ε−F) dx.

Now we estimate

|χ_ε%⁰_εV_ε⁰| ≤χ_εT_n(%⁰_ε)²

%_ε +C_Tnχ_ε%_ε(V_ε⁰)² ≤T_nχ_ε

%_ε(%⁰_ε)²+C(V_ε⁰)², Z ₁

2ε

χ⁰_ε%⁰_εln%_εdx= Z ₁

2ε

χ⁰_ε(%_εln%_ε−%_ε)⁰dx

=−

χ⁰_ε(%_εln%_ε−%_ε) x=2ε−

Z ₁

2ε

χ⁰⁰_ε(%_εln%_ε−%_ε) dx, Z ₁

2ε

χ⁰_ε%⁰_εV_εdx=−(χ⁰_ε%_εV_ε) x=2ε−

Z ₁

2ε

%_ε(χ⁰⁰_εV_ε+χ⁰_εV_ε⁰) dx, and consequently we have

Z 1 2ε

χ_ε

%_ε (ε%⁰⁰_ε)²+T_n(%⁰_ε)²

dx≤C.

Choose a positive σ(ε). Then Lemma 4.2 brings us to Z 1

2ε+σ(ε)

(ε%⁰⁰_ε)²+Tn(%⁰_ε)²dx≤ C σ(ε).

We clearly have ε²%⁰⁰⁰_ε = (%_ε(2T_nln%_ε+V_ε−F))⁰, hence ε²

%⁰⁰⁰_ε

L²((2ε+σ(ε),1)) ≤ C σ^1/2(ε).

Interpolating the estimates k%_εkW2³ ≤Cε⁻²σ^−1/2(ε) and k%_εkW2¹ ≤Cσ^−1/2(ε), we then derive k%εkW2²((2ε+σ(ε),1)) ≤ C

εσ^1/2(ε). By Nirenberg–Gagliardo interpolation,

u⁰

L^∞ ≤Ckuk^1/2_W²

2 kuk^1/2_W¹

2 ≤Ckuk^3/4_W²

2 kuk^1/4_L² , (4.7)

we then have %⁰_ε

L^∞((2ε+σ(ε),1))≤ C ε^1/2σ^1/2(ε). Now the differential equation implies

k%_ε(2T_nln%_ε+V_ε−F)kL²((2ε+σ(ε),1)) ≤ Cε σ^1/2(ε) Then (4.2) follows from choosing σ(ε) = √

ε and interpolating this estimate with the inequality k%_ε(2T_nln%_ε+V_ε−F)kW2¹((2√

ε,1))≤Cε^−1/4; and this choice ofσ(ε) also gives (4.6). Finally, (4.4) is found when we take σ(ε) = 2ε.

4.2 First Remainder Estimates

We continue our preparations for the proof of Theorem 2.2.

Lemma 4.4. The sequence(V_ε)_ε→0 converges to a limit V_∗ ∈C²([0,1]),

kV_ε−V_∗k_W2¹_((0,1)) ≤Cε^1/2, (4.8)

and V_∗ solves (2.2). Moreover, the sequence (%_ε)_ε→0 converges uniformly on compact subsets of (0,1) to the limit %_∗:=√n_∗ in the sense of

k%_ε−%_∗kL^∞((2√

ε,1))=k%_ε−exp((F −V_∗)/(2T_n))kL^∞((2√

ε,1))≤Cε^1/4. (4.9)

See (2.2) for the definition of n_∗.

Proof. For parameters 0< ε₂ < ε₁ <1/8, we conclude that

−λ²(V_ε1−V_ε2)⁰⁰=%²_ε1 −%²_ε2−exp(V_ε1/T_p) + exp(V_ε2/T_p), λ²

Z 1 0

(V_ε1 −V_ε2)⁰2

dx+λ²β(V_ε1(0)−V_ε2(0))²

=− Z 1

0

(exp(V_ε1/T_p)−exp(V_ε2/T_p)) (V_ε1−V_ε2) dx +

Z 4ε1

0

(%²_ε1 −%²_ε2)(V_ε1 −V_ε2) dx+ Z 1

4ε¹

(%²_ε1 −%²_ε2)(V_ε1−V_ε2) dx, and then the representation (4.5) and monotonicity arguments bring us to

λ² Z 1

0

(V_ε1 −V_ε2)⁰2

dx≤Cε₁,

hence there is a limit V_∗ ∈ W₂¹((0,1)), and (4.8) holds. Combined with the uniform bound kV_εkC²([0,1])≤C from Lemma 4.1, then also lim_ε→0kV_ε−V_∗kC¹([0,1])= 0.

Now Lemma 4.1, Lemma 4.2, and (4.2) give us k%_ε−exp((F −V_ε)/(2T_n))kL^∞((2√

ε,1))≤Cε^1/4,

and combining this estimate with (4.8) then yields (4.9). It follows that V_∗ is a distributional solution to the differential equation (2.2), and then V_∗ ∈C²([0,1]) by elliptic regularity.

Next we discuss the behaviour of %_ε near the boundaryx= 0, and our result is:

Proposition 4.5. On the interval [0,4ε], we have the uniform expansion

%_ε(x) =%_∗(x)Zp 2T_n· x

ε

+O ε^1/2· x

ε

, (4.10)

%_ε−%_∗(0)Z

√2Tn

ε ·

C¹([0,4ε]) ≤Cε^1/2, (4.11)

with Z =Z(y) as in (2.3). And on the intermediate interval [4ε,3√

ε], we uniformly have

%_ε(x) =%_∗(x)Zp 2T_n· x

ε

+O(ε^1/4). (4.12)

Proof. Observe that, for 0< x≤1,

ε²%⁰⁰_ε(x) =g(%_ε) +%_ε(V_ε−F), ε²

2 (%⁰_ε)²₀

=g(%_ε)%⁰_ε+1

2(%²_ε)⁰(V_ε−F)

= 2Tn

1

2%²_εln%ε−1 4%²_ε

₀ +1

2(%²_ε)⁰(Vε−F), ε² %⁰_ε(x)2

−ε² %⁰_ε(0)2

= 2T_n

%²_ε(x) ln%_ε(x)−1 2%²_ε(x)

+%²_ε(x)(V_ε(x)−F)

− Z x

0

%²_ε(t)V_ε⁰(t) dt

= 2T_n%²_ε(x)

ln%_ε(x)−1

2 +Vε(x)−F 2T_n

−V_ε⁰(ξ) Z x

0

%²_ε(t) dt, (4.13) withξ ∈(0, x). If 0≤x≤3√

ε, then (4.8) andV_∗ ∈C²([0,1]) imply V_ε(x) =V_∗(0) + (V_ε(x)−V_∗(x)) + (V_∗(x)−V_∗(0))

=V_∗(0) +O(ε^1/2) =F −2T_nln%_∗(0) +O(ε^1/2).

In the following, let R_ε be a generic function of x withkR_εkL²((2√ ε,3√

ε))≤Cε^3/4. For such R_ε, we then also obtain kRεkL¹((2√

ε,3√

ε))≤Cε.

On the interval [2√ ε,3√

ε], we then have, by (4.6) and (4.3), ε² %⁰_ε(x)2

≤Cε^1/2,

%²_ε(x) = exp

F −V_ε(x) T_n

+R_ε(x) =%²_∗(0) +R_ε(x), ln%_ε(x) = F −V_ε(x)

2T_n +R_ε(x),

V_ε⁰(ξ) Z x

0

%²_ε(t) dx

≤Cε^1/2.

Therefore, for 2√

ε≤x≤3√

ε, we conclude that

−ε² %⁰_ε(0)2

= 2Tn%²_ε(x)

−1

2+Rε(x)

+O(ε^1/2), Z 3√

ε 2√

ε

ε² %⁰_ε(0)2

−T_n%²_ε(x) dx≤

Z 3√ ε 2√

ε

Cε^1/2dx+kR_εk_L¹₍₍₂^√_ε,3^√_ε))≤Cε, Z 3√

ε 2√

ε

%²_∗(0)−%²_ε(x) dx=

Z 3√ ε 2√

ε |R_ε(x)|dx≤Cε,

and then we ultimately deduce that ε²(%⁰_ε(0))² = Tn%²_∗(0) +O(ε^1/2) and (under the assumption x∈[2√

ε,3√

ε]) that ε²(%⁰_ε(0))² =T_n%²_∗(x) +O(ε^1/2). Now the differential equation becomes ε² %⁰_ε(x)2

= 2Tn

%²_ε(x) ln%ε(x) +1

2%²_∗(0)− 1

2%²_ε(x)−%²_ε(x) ln%_∗(0)

+O(ε^1/2)

=T_n%²_∗(0)

%²_ε(x)

%²_∗(0)ln

%²_ε(x)

%²_∗(0)

+

1− %²_ε(x)

%²_∗(0)

+O(ε^1/2), 0≤x≤3√ ε.

We introduce the scaling y= x

ε, %_ε(x) =:%_(I,ε)(y)·%_∗(0),

and then we are led to a discussion of a differential equation

%⁰_(I,ε)(y)2

=H(%²_(I,ε)(y)) +O(ε^1/2), 0≤y≤3ε^−1/2, (4.14)

compare (3.7) for the definition of H. Now we restrict our interest to the shorter interval [0,4]. By Lemma 4.2, we have

%⁰_(I,ε)(y)≥ g_∗

2%_∗(0), 0≤y ≤4, (4.15)

and this brings us to

%⁰_(I,ε)(y) =q

H(%²_I,ε(y)) +O(ε^1/2), 0≤y≤4.

Note that the radicand is bounded from below by a positive constant, by (4.15), which makes the functions7→p

H(s²) uniformly Lipschitz continuous on the relevant interval.

The mentioned functionZ =Z(y) from (2.3) solves Z⁰(y) = 1

√2T_n

pH(Z²(y)), 0≤y <∞, Z(0) = 0.

PutW(y) =%_(I,ε)(y)−Z(√

2Tny). Then we have (W²(y))⁰ = 2W(y)·

q

H(%²_(I,ε)(y))− q

H(Z²(p

2Tny)) +O(ε^1/2)

≤C0W²(y) +O(ε),