Juliette Bouhours
email: juliette.bouhours@ljll.math.upmc.fr

Bistable traveling wave passing an obstacle: perturbation results

Keywords: 2010 Mathematics Subject Classification: 35B40, 35K57, 35C07, 35B51, 35B53 parabolic equation, generalized transition front, obstacle, maximum principle

We study the existence of generalized transition fronts for a bistable reaction diffusion equation, u t -∆u = f (u), in a heterogeneous medium, Ω = R N \K where K is a compact set of R N , with Neumann boundary condition and t ∈ R. In the paper, Bistable traveling waves around an obstacle (2009), H. Berestycki, F. Hamel and H. Matano prove the existence of a generalized transition front when K is smooth enough and satisfies some geometric properties. We are interested in an extension of this result when Ω ε = R N \K ε and K ε is a small perturbation of K. We prove that as soon as K ε is close to K in the C 2,α topology generalized transition front still exist while it does not if the perturbation is not smooth enough.

1 Introduction and main results

Problem and motivations

This work is concerned with the following parabolic semilinear problem:

   ∂u ∂t = ∆u + f (u) in Ω ε = R N \K ε , ν • ∇u = 0 on ∂Ω ε = ∂K ε , (1.1)
considering a family of obstacles (K ε ) 0<ε≤1 which are C 2,α -compact sets in R N , for 0 < α < 1, such that for all 0 < ε ≤ 1, K ε ⊂ B R 0 , for some given R 0 > 0, and K ε → K for the C 2,α topology, where K is a C 2,α -compact subset of R N and is either star-shaped or directionally convex (see the definitions at the end of this section).

Remark 1.1 When we write K ε → K for the C 2,α topology we mean that for each x 0 in ∂K, and for some r > 0 such that ∂K ε ∩ B r (x 0 ) = ∅ there exists a couple of parametrization of K ε and K, ψ ε and ψ, C 2,α (B r (x 0 )) functions such that ψ εψ C 2,α (Br(x 0 )) → 0 as ε → 0.

For more details about the C 2,α topology one can look at [START_REF] Gilbarg | Elliptic partial differential equations of second order[END_REF], chapter 6.

We assume that f is of bistable type, meaning that it is a C 1,1 ([0, 1]) function such that, ∃ θ ∈ (0, 1)| f (0) = f (θ) = f (1) = 0, f ′ (0) < 0, f ′ (1) < 0, (1.2) f (s) < 0 ∀ s ∈ (0, θ), f (s) > 0 ∀ s ∈ (θ, 1), with positive mass:

∀ 0 ≤ s < 1, 1 s f (τ )dτ > 0.
(1.3)

Our goal is to study how does the shape of K influence the behavior of the solutions of (1.1) and as a consequence how a propagating planar traveling front interacts with our obstacles (K ε ) ε .

Existence of traveling fronts for the homogeneous equation u t = ∆u + f (u) in R N has been studied since the article of Kolmogorov, Petrovsky and Piskunov, [START_REF] Kolmogorov | Etude de l'équation de la diffusion avec croissance de la quantité de la matière et son application à un problème biologique[END_REF] in 1937. Recently more and more interest has been observed in the study of traveling front solutions in spatially dependent media because of its importance in several scientific fields. For more descriptions and references about those fronts one can look at the introduction of [START_REF] Berestycki | Bistable traveling waves around an obstacle[END_REF], and references therein. As in H. Berestycki, F. Hamel and H. Matano's paper we study here the existence of a particular type of generalized transition wave, a notion introduced by H. Berestycki and F. Hamel in [START_REF] Berestycki | Fronts and invasions in general domains[END_REF][START_REF] Berestycki | Generalized travelling waves for reaction-diffusion equations[END_REF][START_REF] Berestycki | Generalized transition waves and their properties[END_REF]] (see at the end of Section 1.2, Remark 1.7 for a precise definition).

Coming back to our framework, it follows from (1.2) and (1.3) that there exists a unique solution (up to translation) (c, φ) of

     φ ′′ (z) -cφ ′ (z) + f (φ(z)) = 0, z ∈ R, φ(-∞) = 0, φ(+∞) = 1, 0 < φ(z) < 1, z ∈ R, (1.4)
with c > 0. For further details see [START_REF] Paul | The approach of solutions of nonlinear diffusion equations to travelling front solutions[END_REF] for example. As stated in [START_REF] Berestycki | Bistable traveling waves around an obstacle[END_REF] (section 2 and 3), as soon as K ε ⊂ {x, x 1 ≤ L}, L ∈ R, there exists a unique time global solution u ε of (1.1) such that 0 < u ε < 1, (u ε ) t > 0 in R × Ω ε and that behaves like a planar traveling front for very small times, i.e

u ε (t, x) -φ(x 1 + ct) → 0 as t → -∞ uniformly in x ∈ Ω ε . (1.5)
To prove these results the authors constructed some sub and super solutions ω + and ω - depending only on x 1 and t. Then using (1.5) and some extension of the comparison principle they proved the uniqueness of the entire (time-global) solution u ε . Moreover, they proved that the behavior of u ε for large time is determined by the solution of the associated stationary problem. We will then start with the study of the following elliptic problem:

           ∆u ε,∞ + f (u ε,∞ ) = 0 in Ω ε = R N \K ε , ν • ∇u ε,∞ = 0 on ∂Ω ε = ∂K ε , 0 < u ε,∞ ≤ 1 in Ω ε , lim |x|→+∞ u ε,∞ (x) = 1. (1.6) Remark 1.2 We add the extra consition lim |x|→+∞ u ε,∞ (x) = 1, because from [3], section 5, as soon as K ε is a compact set, f ∈ C 1 ([0, 1]) satisfies (1.
2) and (1.3), (φ, c) solution of (1.4) with c > 0 and under some asymptotic and monotonicity conditions on u ε , which are satisfied in our framework, our entire solution u ε converges toward u ε,∞ as t → +∞ locally uniformly in x ∈ Ω ε , and u ε,∞ is a classical solution of

     ∆u + f (u) = 0 in Ω ε , ν • ∇u = 0 on ∂Ω ε , 0 < u eq1 in Ω ε , (1.7
)

such that u ε,∞ (x) → 1 as |x| → +∞.
Before stating the main results, let explain what we mean by star-shaped or directionally convex obstacles.

Definition 1.3 K is called star-shaped, if either K = ∅, or there is x ∈

• K such that, for all y ∈ ∂K and t ∈ [0, 1), the point x + t(yx) lies in • K and ν K (y) • (yx) ≥ 0, where ν K (y) denotes the outward unit normal to K at y. Definition 1.4 K is called directionally convex with respect to a hyperplane P if there exists a hyperplane P = {x ∈ R N , x • e = a} where e is a unit vector and a is some real number, such that • for every line Σ parallel to e the set K ∩ Σ is either a single line or empty,

• K ∩ P = π(K) where π(K) is the orthogonal projection of K onto P .

Main results

Our main result is the following theorem

Theorem 1.5 Assume f ∈ C 1,1 ([0, 1]) satisfies (1.2) and (1.3). Consider (K ε ) 0<ε≤1 a family of C 2,α -compact sets of R N (i.e for all ε ∈ (0, 1], K ε ⊂ B R 0 for some R 0 > 0) and let Ω ε = R N \K ε a smooth, open, connected subset of R N (with N ≥ 2). Assume that K ε → K
for the C 2,α topology as ε → 0, where K is a C 2,α -compact subset of R N that is either starshaped or directionally convex with respect to some hyperplane P. Then there exists ε 0 > 0 such that for all 0 < ε < ε 0 , the unique solution of (1.6) is u ε,∞ ≡ 1

This means that for obstacles that are compact sets in R N and close enough (in the C 2,α sense) to some star-shaped or directionally convex domains then the unique solution of (1.6) is the constant 1. This Theorem yields some properties about the solution u ε of the parabolic problem (1.1).

Corollary 1.6 Assume that f satisfies (1.2) and that there exists a solution φ to (1.4) with c > 0 (if f satisfies (1.2) and (1.3), φ exists and c > 0). Let (K ε ) 0<ε≤1 be a family of compact domains in R N such that for all 0 < ε ≤ 1, K ε ⊂ B R 0 , for some given R 0 > 0, and K ε → K for the C 2,α topology, with 0 < α < 1 and where K ⊂ R N is either star-shaped or directionally convex with respect to some hyperplane P. Then for all 0 < ε ≤ 1, there exists an entire solution u ε (t, x) of (1.1) such that 0 < u ε < 1 and ∂ t u ε > 0 over R × Ω ε and there exists ε 0 > 0 such that for all 0 < ε < ε 0 ,

u ε (t, x) -φ(x 1 + ct) → 0
as t → ±∞ uniformly in x ∈ Ω ε , and as |x| → +∞ uniformly in t ∈ R.

We will prove Theorem 1.5 in section 2 below, and Corollary 1.6 in section 3.

Remark 1.7 For all 0 < ε < ε 0 the solution u ε (t, x) given in Corollary 1.6 is a generalized, almost planar, invasion front between 0 and 1 with global mean speed c, in the sense that

sup (t,x)∈R×Ωε,x 1 +ct≥A |u ε (t, x) -1| ----→ A→+∞ 0 sup (t,x)∈R×Ωε,x 1 +ct≤-A |u ε (t, x)| ----→ A→+∞ 0
Before proving the previous statements, let give some examples of domains (K ε ) ε and K to illustrate our results.

Examples

We assume that N = 2 and we construct two families of obstacles; one which converges to a star shaped domain and the other which converges to a directionally convex domain. The black plain line represents the limit K and the dashed parts represent the small perturbations (of order ε).

Star-shaped domain: The long-dashed line is used during the construction of K and it is on this line that we could find the center(s) of the domain (i.e the point x in Definition 1.3). We can clearly see that for all ε > 0, K ε is not star-shaped, because the points just behind the dashed area cannot be linked to the center x by a straight line.

Directionally convex domain: In this second figure, the hyperplanes P , for which the plain black domain could be directionally convex are necessarily horizontal (i.e {(x, y) ∈ R 2 |(x, y) • (0, 1) = a}) and a has to be 0 (assuming the center of the ellipses are on the x-axis), else the second property in Definition 1.4 is not satisfied. Indeed the first property of Definition 1.4 eliminates every vertical hyperplanes P and the second property eliminates the diagonal ones. Adding small perturbations (of order ε) on each side of each ellipse, one gets that for all ε > 0 K ε does not satisfy the second property of Definition 1.4.

One need to be careful on the shape of the perturbations. Indeed considering an ellipse, which is star-shaped and directionally convex and adding on each side of the vertical axis some well chosen perturbations (keeping the domain smooth), the obstacle is not star-shaped, neither directionally convex anymore (see figure below), Figure 3: Obstacles converging only in C 0 and K ε → K, an ellipse, as ε → 0, but the convergence of K ε cannot be C 2,α (see section 4 for more details) but only C 0 which is not enough to get Theorem 1.5, as we will see in section 4.

Proof of the main Theorem

In this section, to simplify the notation we will write u ε in stead of u ε,∞ . As we are only working on the stationary problem, there is no confusion with the solution of the parabolic system. To prove Theorem 1.5, we will use the following Proposition:

Proposition 2.1 For all 0 < δ < 1, if u ε is a solution of (1.6), then there exists R = R(δ) > R 0 , such that u ε (x) ≥ 1 -δ for all |x| ≥ R and for all 0 < ε < 1.
This proposition means that u ε converges toward 1 as |x| → +∞ uniformly in ε.

Let first admit this result and prove Theorem 1.5.

Proof of Theorem 1.5

Using Schauder Estimates (see [START_REF] Gilbarg | Elliptic partial differential equations of second order[END_REF], chapter 6.7), we obtain that

u ε C 2,α (Ω) ≤ C.
As K ε → K for the C 2,α topology when ε → 0, C is independent of ε. So there exists a sequence (ε n ) n in ]0, η ∧ 1[ such that ε n → 0 and u εn → u in C 2 loc as n → +∞ and u satisfies:

∆u + f (u) = 0 in Ω = R N \K, ν • ∇u = 0 on ∂Ω = ∂K, (2.1)
because of the compact injection C 2,α loc ֒→ C 2 loc . Using Proposition 2.1 we get lim

|x|→+∞ u(x) = 1.
And K is either star-shaped or directionally convex. We now recall the following results from [START_REF] Berestycki | Bistable traveling waves around an obstacle[END_REF]:

Theorem 2.2 Let f be a Lipschitz-continuous function in [0, 1] such that f (0) = f (1) = 0 and f is nonincreasing in [1 -δ, 1] for some δ > 0. Assume that ∀ 0 ≤ s < 1, 1 s f (τ )dτ > 0. (2.2)
Let Ω be a smooth, open, connected subset of R N (with N ≥ 2) with outward unit normal ν, and assume that

K = R N \Ω is compact. Let 0 ≤ u ≤ 1 be a classical solution of      ∆u = f (u) in Ω, ν • ∇u = 0 on ∂Ω, u(x) → 1 as |x| → +∞. (2.3)
If K is star shaped or directionally convex, then u ≡ 1 in Ω.

(2.4)

Remark 2.3 The main difficulty here was that, for every ε ∈ (0, 1), K ε is not necessarily star-shaped nor directionally convex and thus we could not use directly Theorem 2.2. The purpose of the article is to find an equivalent to Theorem 2.2 but in the case of a family of obstacles (K ε ) ε that converges toward a star-shaped or directionally convex domain.

It follows from Theorem 2.2 that u ≡ 1. It also proves that the limit u is unique and thus u ε → u as ε → 0 in C 2 (and not only along a subsequence).

Now we need to prove that there exists ε 0 > 0 such that u ε ≡ 1 for all 0 < ε < ε 0 . Let assume that for all ε > 0, u ε ≡ 1. Then there exists

x 0 ∈ Ω ε such that u ε (x 0 ) = min x∈Ωε u ε (x) < 1.
As u ε is a solution of (1.6), the Hopf lemma yields that,

if x 0 ∈ ∂K ε then δu ε δν (x 0 ) < 0,
which is impossible due to Neuman boundary conditions. Hence

x 0 ∈ Ω ε . If u ε (x 0 ) > θ, -∆u ε (x 0 ) = f (u ε (x 0 )) > 0,
which is impossible since x 0 is a minimizer. So, for all 0 < ε < 1,

0 ≤ min x∈Ωε u ε (x) ≤ θ,
which contradicts u ε → u ≡ 1. Thus there exists ε 0 such that for all ε < ε 0 , u ε ≡ 1.

Proof of Proposition 2.1

We will now prove Proposition 2.1, starting with the following lemma:

Lemma 2.4 There exists ω = ω(r) with r ∈ R + such that

           -ω ′′ (r) = f (ω(r)), ∀r ∈ R + * , ω(0) = 0, ω ′ (0) > 0, ω ′ > 0, 0 < ω < 1 in R * + , lim r→+∞ ω(r) = 1. (2.5)
We can prove this lemma using some results about traveling fronts in the multistable case (a) or using an ODE approach (b).

Proof of Lemma 2.4

• (a) Proof using traveling waves We know from [START_REF] Paul | The approach of solutions of nonlinear diffusion equations to travelling front solutions[END_REF] that there exists a unique solution (ω

1 , c 1 ) ∈ C 2 (R) × R of          -ω ′′ 1 (x) + c 1 ω ′ 1 (x) = f (ω 1 (x)) in R, ω 1 (-∞) = 0, ω 1 (+∞) = 1, ω ′ 1 > 0, in R, 0 < ω 1 < 1, in R.
(2.6)

Because of (1.3), c 1 > 0. One can prove the existence and uniqueness of v

∈ C 2 (R + )          v ′′ (z) -cv ′ (z) + f (v(z)) = 0, ∀z ∈ (0, +∞), v(0) = 0, v(+∞) = 1, v ′ (z) > 0, ∀z ∈ (0, +∞), 0 < v(z) < 1, ∀z ∈ (0, +∞), (2.7)
for c ≤ c 1 (see [START_REF] Paul | The approach of solutions of nonlinear diffusion equations to travelling front solutions[END_REF][START_REF] Hamel | Existence of multidimensional travelling fronts with a multistable nonlinearity[END_REF]). Then taking c = 0 ≤ c 1 there exists ω ∈ C 2 (R + ) such that

         ω ′′ (z) + f (ω(z)) = 0, ∀z ∈ (0, +∞), ω(0) = 0, ω(+∞) = 1, ω ′ (z) > 0, ∀z ∈ (0, +∞), 0 < ω(z) < 1, ∀z ∈ (0, +∞).
The Hopf Lemma yields that ω ′ (0) > 0.

• (b) Proof using an ODE approach.

We want to prove the global existence and uniqueness of the following ODE:

     -ω ′′ = f (ω) in (0, +∞), (1) ω(0) = 0, ω ′ (0) = 2F (1), (2.8)
where

F (z) = z 0 f (s)ds. Using (1.3), F (1) > 0.
From Cauchy-Lipschitz theorem we know that there exists a unique maximal solution ω of (2.8) in I ⊂ (0, +∞). To prove the global existence, i.e I = (0, +∞), let prove that ω ′ > 0 in I and 0 < ω < 1 in I.

We start by proving that ω ′ > 0. We know that ω ′ (0) > 0. Suppose that there exists r 0 ∈ I such that ω ′ (r 0 ) = 0. Then multiplying (2.8) (1) by ω ′ and integrating between 0 and r 0 , one gets:

-F (1) + F (w(r 0 )) = 0 = ω(r 0 ) 1 f (z)dz.
(2.9)

Without loss of generality, we extend f linearly (as a C 1 function) outside [0, 1]. Note that by the Maximum Principle any solution with such a f will take values in [0, 1], hence is a solution of the original problem. The last equation (2.9) is impossible, which implies that ω > 0 in I.

Next assume by contradiction that there exists r 1 ∈ I such that ω(r 1 ) = 1. Using the same method as above (multiplying by ω ′ and integrating between 0 and r 1 ) one gets:

(ω ′ (r 1 )) 2 2 -F (1) + F (1) = 0,
which is impossible. Hence 0 < w < 1 in I.

If we assume that I (0, +∞), i.e there exists r ∞ ∈ (0, +∞) such that I = (0, r ∞ ), it means that lim ω(r) = +∞ as r → r ∞ . This is impossible because 0 < ω < 1 in I. Thus I = (0, +∞).

We have proved that there exists a unique global solution ω of (2.8) and that ω ′ (r) > 0 and 0 < ω(r) < 1 for all r ∈ (0, +∞). As ω is increasing and bounded from above, it has a limit when r → +∞ such that 0 < ω(+∞) ≤ 1. Moreover ω(+∞) > θ. Indeed if we assume that ω ≤ θ in R + then one gets that ω is convex (since ω ′′ = -f (ω) ≥ 0) and increasing, it then goes to +∞ when r → +∞, which is impossible. It immediately follows from elliptic regularity estimates that f (ω(+∞)) = 0. Hence ω(+∞) = 1. One has proved Lemma 2.4.

Proof of Proposition 2.1. Now we introduce a function f δ with the same hypothesis as f but such that

f δ ≤ f , f δ = f in [0, 1 -δ] and f δ (1 -δ 2 ) = 0 (see figure below). Notice that 1-δ 2 0 f δ (z)dz > 0 for δ small. 1 1 -δ 2 1 -δ x f f δ Figure 4: f δ (x)
For the same reason than in Lemma 2.4 there exists ω = ω δ such that

         -ω ′′ δ (x) = f δ (ω δ (x)) in (0, +∞), ω δ (0) = 0, ω δ (+∞) = 1 -δ 2 , 0 < ω δ < 1 -δ 2 in (0, +∞), ω ′ δ > 0 in (0, +∞).
(2.10)

Next, for any R > R 0 consider ω δ (|x| -R) for every |x| ≥ R and let z(x) = ω δ (|x| -R), one gets:

∆z + f (z) = ω ′′ δ + N -1 |x| ω ′ δ + f (ω δ ) = N -1 |x| ω ′ δ + f (ω δ ) -f δ (ω δ ) > 0 in {|x| > R}. So -∆z < f (z) in {|x| > R}.
(2.11)

We want to prove that

ω δ (|x| -R 0 ) < u ε (x), ∀x ∈ R N , |x| ≥ R 0 .
We know from (1.6) that u ε (x) → 1 as |x| → +∞. Hence there exists

A = A(ε) > 0 such that u ε (x) ≥ 1 -δ 3 , for all |x| ≥ A. One gets u ε (x) ≥ ω(|x| -A), for all |x| ≥ A. A 1 1 -δ 2 1 -δ 3 u ω δ (|x| -A) Figure 5: u ε and ω δ (|•| -A) Consider R = inf R ≥ R 0 ; u ε (x) > ω δ (|x| -R), for all {|x| ≥ R} . (2.12) As R ≥ R 0 and K ε ⊂ B R 0 , u ε is always defined on {|x| > R}. One will prove that R = R 0 .
As ω δ is increasing, we know that

∀R ≥ A u ε (x) ≥ ω(|x| -R), ∀|x| ≥ R. Hence R ≤ A.
Assume that R > R 0 . Then there are two cases to study:

• either inf u ε (x) -ω δ (|x| -R), ∀|x| > R > 0, ( 1
)
• or inf u ε (x) -ω δ (|x| -R), ∀|x| > R = 0. ( 2
)
In the first case (1), one gets u ε (x) > ω δ (|x| -R) for all |x| > R. As ∇u ε and ω ′ δ are bounded one can translate ω δ to the left such that both graphs touch at one point, i.e there exists R * < R such that u ε (x) ≥ ω δ (|x| -R * ) for all |x| > R * , and u ε (x 0 ) = ω δ (|x 0 | -R * ) for some |x 0 | > R * . This contradicts the optimality of R.

In the second case (2), there necessarily exists

x 0 with |x 0 | > R such that u ε (x 0 ) = ω δ (|x 0 | - R). Let v(x) = u ε (x) -ω δ (|x| -R), for all |x| > R.
As u ε is a solution of (1.6) and using (2.11), v satisfies:

-∆v > f (v) in {|x| > R}, v > 0 on {|x| = R}. (2.13)
From the maximum principle v(x) ≥ 0, for all |x| ≥ R. But there exists x 0 such that |x 0 | > R (x 0 is an interior point) and v(x 0 ) = 0, i.e v(•) reaches its minimum 0 inside the domain. This implies that v(•) ≡ 0, which is impossible because v(•) > 0, for all |x| = R.

Then R = R 0 which does not depend on ε and

∀|x| ≥ R 0 u ε (x) ≥ ω(|x| -R 0 ). But ω δ (x) → 1 -δ 2 as |x| → +∞ implies that there exists R, independent of ε, such that for all |x| > R + R 0 , u ε (x) > ω δ (|x| -R 0 ) ≥ 1 -δ. One has proved Proposition 2.1.

The associated parabolic problem and its properties

In this section we will use Theorem 1.5, to apply some results of [START_REF] Berestycki | Bistable traveling waves around an obstacle[END_REF] and derive Corollary 1.6. We investigate the following semilinear parabolic problem:

   ∂u ε ∂t = ∆u ε + f (u ε ) in Ω ε = R N \K ε , ν • ∇u ε = 0 on ∂Ω ε = ∂K ε , (3.1)
where K ε is a compact set in R N . Notice that in this section, u ε is the solution of the parabolic problem. Proof of Corollary 1.6 We know from H. Berestycki, F. Hamel and H. Matano's paper [START_REF] Berestycki | Bistable traveling waves around an obstacle[END_REF] that their exists an entire solution

u ε of (3.1) in Ω ε such that 0 < u ε (t, x) < 1, ∂ t u ε (t, x) > 0 for all (t, x) ∈ R × Ω ε and u ε (t, x) -φ(x 1 + ct) → 0 as t → -∞ uniformly in x ∈ Ω ε ,
and as |x| → +∞ uniformly in t ∈ R. Furthermore, there exists a classical solution u ε,∞ of

           ∆u ε,∞ + f (u ε,∞ ) = 0 in Ω ε = R N \K ε , ν • ∇u ε,∞ = 0 on ∂Ω ε = ∂K ε , 0 < u ε,∞ ≤ 1 in Ω ε , lim |x|→+∞ u ε,∞ (x) = 1, (3.2) such that u ε (t, x) -φ(x 1 + ct)u ε,∞ (x) → 0 as t → +∞ uniformly in x ∈ Ω ε .
Then using Theorem 1.5, there exists ε 0 < 1 such that for all 0 < ε < ε 0 , the only solution u ε,∞ of (3.2) is identically equal to 1. We have proved Corollary 1.6.

4 Discussion about the convergence of (K ε ) 0<ε<1

In this section we discuss the hypothesis of convergence of (K ε ) 0<ε<1 . Until now we assumed that K ε → K in C 2,α , with 0 < α < 1 in order to use the Schauder estimates and ensure the convergence of u ε as ε → 0 . One can wonder if we can weaken this hypothesis: would the C 0 or C 1 convergence be enough?

We will prove that the C 0 convergence is not enough.

Example of a family of obstacles that converges only in C 0

In this subsection we construct a family of obstacles that are neither star-shaped nor directionally convex but converges uniformly to K which is convex. We want to prove that for all ε ∈]0, 1] there exists a solution of (3.2) which is not identically equal to 1. To do so we will use the counterexample of section 6.3 in [START_REF] Berestycki | Bistable traveling waves around an obstacle[END_REF].

K η ε K = B R0 R 0
C ε B R0 R 1 R 2 A β ε x 0 2η 4η
x Figure 7: Zoom on the perturbation V ε

We consider an obstacle K ε = K η ε (see figure 6 and7), for ε > 0, such that:

       A ∩ {x; x 1 ≤ x 0 1 } ∪ B R 0 ∪ C ε ⊂ K η ε , A ∩ x; x 1 > x 0 1 , |x ′ | > 2η ⊂ K η ε , K η ε ⊂ A ∩ x; x 1 > x 0 1 , |x ′ | > η ∪ B R 0 ∪ A ∩ x; x 1 ≤ x 0 1 ∪ C ε . (4.1)
where x ′ = (x 2 , ..., x N ) and

A = {x : R 1 ≤ |x -x 0 | ≤ R 2 }, R 0 , R 1 < R 2
, are three positive constants, x 0 = (x 0 1 , 0, 0, ..., 0) is the center of the annular region A with x 0 1 = R 0 + R 2 + β ε , C ε is some corridor that links smoothly A and B R 0 which length is β ε and η > 0, small enough. Now let explain why the convergence of K ε is only true for the C 0 topology. We want that for all ε > 0, the annular region A stays an annular region (see arguments in the next section, Corollary 4.1 and 4.2). To do so we need to reduce vertically and horizontally the perturbation in order to stay with an annular region. To simplify the proof, assume that N = 2. Let g be the parametrization of K = B R 0 , i.e

K = (x, y) ∈ R 2 |(x, y) = g(t) = R 0 cos(t), R 0 sin(t) ∀t ∈ [0, 2π[ . (4.2)
Let f ε be the parametrization of K ε for all 0 < ε ≤ 1. To define f ε , we start with the case when ε = 1:

K 1 = (x, y) ∈ R 2 |(x, y) = f 1 (t) = g(t) ∀t ∈]θ, 2π -θ[, h(t) ∀t ∈ [0, θ] ∪ [2π -θ, 2π[ , (4.3)
where θ is some small positive number and h is such that f 1 is a C 2,α function. Now one can define f ε and K ε :

K ε = (x, y) ∈ R 2 |(x, y) = f ε (t) = g(t) ∀t ∈]εθ, 2π -εθ[, h ε (t) ∀t ∈ [0, εθ] ∪ [2π -εθ, 2π[ , (4.4)
where h ε is such that f ε is a C 2,α function and such that for every (x, y) ∈ A ∩ K ε , (x, y) = h ε (t) = εh(t). This last condition ensures that A stays an annular region.

Then one can easily see that K ε → K as ε → 0 for the C 0 topology, i.e ff ε C 0 ([0,2π[) → 0 as ε → 0 Now assume that K ε → K as ε → 0 for the C 1 topology. One can notice (see figure 8 below) that for ε = 1 there exists a t m ∈ [0, θ] such that f ′ 1 (t m ) = e 1 = (1, 0) or ν(x m , y m ) = e 2 = (0, 1), where ν is the outward unit normal and (x m , y

m ) = f 1 (t m ) ∈ A ∩ K ε . K = B R0 ν(xm , ym) (xm, ym) (R0, 0) ν(R0, 0) x y θ (R0, 0) Figure 8: Outward unit normal Let consider t ε ∈ [0, εθ] such that ν(x ε , y ε ) = e 2
, where (x ε , y ε ) = f ε (t ε ) the point at the top of the big sphere in A. This point always exists because we parametrize K ε such that A stays an annular region. As the convergence is C 1 one should have that ν(x ε , y ε ) → ν(R 0 , 0) as ε → 0, because (x ε , y ε ) → (R 0 , 0) as ε → 0. But ν(x ε , y ε ) = e 2 for every ε ∈ (0, 1] and ν(R 0 , 0) = e 1 , which is impossible. The convergence can not be C 1 .

Existence of a non constant solution u ε of (1.6)

We want to prove that for all 0 < ε < 1 their exists a solution 0

< u ε < 1 of      -∆u ε = f (u ε ) in R N \K η ε = Ω η ε , ν • ∇u ε = 0 on ∂K η ε = ∂Ω η ε , u ε (x) → 1 as |x| → +∞. (4.5)
Here as in section 2, u ε denotes the solution of the elliptic problem. We will follow the same steps as in [START_REF] Berestycki | Bistable traveling waves around an obstacle[END_REF], section 6. First, let notice that it is enough to find ω ≡ 1 solution of

     -∆ω = f (ω) in B R \K η ε , ν • ∇ω = 0 on ∂K η ε , ω = 1 on ∂B R , (4.6)
And then using Proposition 6.8 of [START_REF] Berestycki | Bistable traveling waves around an obstacle[END_REF] and Corollary 4.1 one gets the following Corollary Corollary 4.2 There exist γ > 0 and η 0 > 0 (which depend on ε) such that for all 0 < η < η 0 and v ∈ H 1 0 such that vv 0 2 = δ, then

J B R \K η ε (v 0 ) < J B R \K η ε (v) -γ. The functional J B R \K η
ε admits a local minimum in the ball of radius δ around v 0 in H 1 (B R \K η ε ) ∩ {v = 0 on ∂B R }. This yields a (stable) solution v of (4.8) for small enough η > 0. Furthermore, provided that δ is chosen small enough, this solution does not coincide either with 1 or with 0 in B R \K η ε . We have proved that for all ε ∈]0, 1], u ε ≡ 1.

One has proved that C 0 convergence of the domain is not sufficient. One can now wonder whether this C 1 convergence would be sufficient. The main difficulties here would be the following • We cannot apply Theorem 1.5 directly to prove the sufficiency of the C 1 convergence.

In deed without the C 2,α convergence of the obstacles we cannot use Schauder theory.

One can try to relax some assumptions on the regularity of the domain in the Schauder theory, as it has already been done in the literature (L ∞ assumption for coefficient in an elliptic or parabolic equation is enough to get the maximum principle). This remark refers to some technical arguments that will not be further explored in this article.

• One can try to construct a counterexample, as for the C 0 convergence. The problem here would be that we could not use the energy arguments anymore because if we want the convergence to be C 1 , the perturbations cannot draw any holes (most important argument in the construction of v 0 ).

Thus the optimal space of convergence for the obstacles K ε is still an open problem.

Figure 1 :

 1 Figure 1: Obstacles converging toward a star-shaped domain

Figure 2 :

 2 Figure 2: Obstacles converging toward a directionally convex domain

Figure 6 :

 6 Figure 6: Liouville counterexample Zooming on the dashed part:

Acknowledgments

The author thanks Henri Berestycki, Grégoire Nadin and François Hamel for their fruitful discussions and encouragements, Sylvain Arguillère for his geometric insights. This study was supported by the French "Agence Nationale de la Recherche" through the project PREF-ERED (ANR 08-BLAN-0313) and by the French "Région Ile de France" through the fellowship "Bouse Hors DIM".

for some R > 0 large enough such that K η ε ⊂ B R . Indeed then ω extended by 1 outside B R is a supersolution of (4.5) and one can define:

where U : R + → (0, 1) satisfies

is a subsolution, ψ is a subsolution of (4.5). Hence there exists a solution ψ < u ε < ω of (4.5). If we prove that ω ≡ 1 then 0 < u ε < 1 (with the maximum principle). Now let consider our problem (4.6) and replace ω by v = 1ω. The problem becomes

Using exactly the same arguments as in [START_REF] Berestycki | Bistable traveling waves around an obstacle[END_REF] one proves that, if we considere:

, (4.9) then there exists v ∈ H 1 (B R \K η ε ) ∩ {v = 0 on ∂B R } = H 1 0 such that vv 0 H 1 0 < δ for some δ > 0 small enough, that is a local minimizer of the associated energy functional in H 1 0 when the width η of the channel is small enough. For more clarity we will give the main step of the proof but for details see [START_REF] Berestycki | Bistable traveling waves around an obstacle[END_REF], section 6.3.

We introduce the energy functional in a domain D:

More precisely, their exist α > 0 and δ > 0 for which

for all v ∈ H 1 (B R 1 (x 0 )) such that vv 0 2 ≤ δ.


P Pietra
email: paola.pietra@imati.cnr.it

N Vauchelet
email: vauchelet@ann.jussieu.fr.

Numerical simulations of an energy-transport model for partially quantized particles

Keywords: Schrödinger equation, energy-transport system, subband model, nanotransistor, Gummel iterations, mixed finite elements AMS Subject Classification : 65M60, 65Z05, 82D37, 82D80, 35J10, 76P05

A coupled quantum-classical model describing the transport of electrons confined in nanoscale semiconductor devices is considered. Using the subband decomposition approach allows to separate the transport directions from the confinement direction. The motion of the gas in the transport direction is assumed to be classical. Then a hierarchy of adiabatic quantum-classical model is obtained, leading to subband SHE and energy-transport models, with explicit expression of the diffusion coefficients. The energy-transport-Schrödinger-Poisson model is then used for the numerical simulation of the transport of the electron gas in an ultra-scaled Double-Gate-MOSFET.

Introduction

In nanoscale semiconductor devices, electrons might be extremely confined in one or several directions refered to as the confining direction. This leads to a partial quantization of the energy which can be modelled by the subband decomposition method [START_REF] Polizzi | Subband decomposition approach for the simulation of quantum electron transport in nanostructures[END_REF]. This subband decomposition approach allows to separate the confinement direction from the transport direction. Thanks to the reduction of the dimension of the transport problem, the computational gain is significant. In the confined direction electrons behave like waves; the system is at thermodynamical equilibrium and it is described by a statistical mixture of eigenstates of a Schrödinger-Poisson system. In the transport direction the transport can be of classical [START_REF] Ben Abdallah | Diffusive transport of partially quantized particles : Existence, uniqueness and long time behaviour[END_REF][START_REF] Ben Abdallah | Adiabatic quantum-fluid transport models[END_REF] or quantum nature [START_REF] Polizzi | Subband decomposition approach for the simulation of quantum electron transport in nanostructures[END_REF].

Here, we are interested in deriving adiabatic quantum-classical models accounting for thermal effects, aiming at accurate and efficient numerical simulation of confined devices. In [START_REF] Bourgade | On quantum extensions to classical spherical harmonics expansion/Fokker-Planck models[END_REF] several spherical harmonic expension (SHE) models incorporating quantum effects are proposed. However, with their strategy the obtained models have a complicated nonlocal structure which is not suitable for numerical purposes. Quantum energy-transport description of the gas, whereas in the third direction z ∈ [0, ℓ] quantum effects occur. The quantum confinement of the electron gas is described thanks to the eigen-elements of the 1D Schrödinger operator. They are denoted (ǫ n , χ n ) n∈N * and solve the eigenvalue problem :

       - 2 2 d dz 1 m * (z) d dz χ n + Uχ n = ǫ n χ n , χ n (x, •) ∈ H 1 0 (0, ℓ), ℓ 0 χ n χ n ′ dz = δ nn ′ .
(2.1)

In this equation is the reduced Planck constant, m * the effective mass. It is known that the eigenvalues ǫ n in (2.1) form an increasing sequence tending to +∞. These functions depend on the potential energy defined by U = -eV , where e is the elementary charge and V denotes the self-consistent electrostatic potential, solution of the Poisson equation

div x,z (ε R (x, z)∇ x,z V ) = e ε 0 (N e -N D ). (2.2)
Here ε R (x, z) denotes the relative permittivity, ε 0 the permittivity constant in vacuum, N D (x, z) is the prescribed doping density and N e (t, x, z) is the electron density. This density is described by a sequence of distribution functions (f n ) n∈N * describing the repartition on each subband for the classical direction x ∈ ω and the corresponding momentum variable k ∈ R 2 . It is written as

N e (t, x, z) = +∞ n=1 R 2 f n (t, x, k) dk |χ n | 2 (t, x, z).
The evolution of distribution functions is governed by classical transport model in the x direction parallel to the gas. The total energy of the nth subband is defined by

ε n (t, x, k) = |k| 2 2
2m * + ǫ n (t, x).

(2.3)

Therefore the energy-band diagram of the semiconductor crystal is spherically symmetric and strictly monotone with respect to |k|. Then the Brillouin zone (which represents the elementary cell of the dual lattice L * ) is equal to R 2 . Moreover, we point out that, in contrast with the classical counterpart, the energy-band depends on space and time. In a kinetic collisional framework, the distribution function f n of the nth subband satisfies the rescaled Boltzmann transport equation [START_REF] Ben Abdallah | An energy transport model for the semiconductors derived from the Boltzmann equation[END_REF][START_REF] Pietra | Modeling and simulation of the diffusive transport in a nanoscale Double-Gate MOSFET[END_REF] :

α 2 ∂ t f α n + α(∇ k ε n • ∇ x f α n -∇ x ε n • ∇ k f α n ) = Q ld (f α ) n + α 2 β Q e (f α ) n , (2.4)
where Q ld is the collision operator for the lattice defect collisions and Q e is the collision operator for the elastic, non linear electron-electron collisions, accounting for intra-band scattering as well as for transitions between subbands. α and β are dimensionless parameters that satisfy α ≪ β ≪ 1.

The main classes of lattice-defects that we shall consider are impurities and phonons [START_REF] Ashcroft | Solid State Physics[END_REF] :

Q ld (f ) = Q imp (f ) + Q ph (f ).
The elastic character of the impurity scattering leads to

Q imp (f ) n (k) := n ′ ∈N * R 2 Φ imp n,n ′ (k, k ′ )δ(ε n (k) -ε n ′ (k ′ ))(f n ′ (t, x, k ′ ) -f n (t, x, k)) dk ′ ,
where δ is the Dirac measure and the dependence on t, x of ε and Φ imp has been omitted. The cross-section is assumed to be symmetric:

Φ imp n,n ′ (k, k ′ ) = Φ imp n ′ ,n (k ′ , k
). The electron-phonon collision operator is considered as

Q ph (f ) n (k) = n ′ ∈N * R 2 Φ ph n,n ′ (k, k ′ ) ([(N ph + 1)δ(ε n (k) -ε n ′ (k ′ ) + α 2 ε ph ) + N ph δ(ε n (k) -ε n ′ (k ′ ) -α 2 ε ph )]f n ′ (k ′ )(1 -ηf n (k))- [(N ph + 1)δ(ε n ′ (k ′ ) -ε n (k) + α 2 ε ph ) + N ph δ(ε n ′ (k ′ ) -ε n (k) -α 2 ε ph )]f n (k)(1 -ηf n ′ (k ′ )))dk ′ ,
where again,

Φ ph n,n ′ (k, k ′ ) = Φ ph n ′ ,n (k ′ , k
), ε ph is the phonon energy, η ≥ 0 is a dimensionless distribution function scaling factor and the terms 0 ≤ 1-ηf n ≤ 1 express the Pauli exclusion principle. N ph is the phonon occupation number, given by the Bose-Einstein statistics

N ph = e α 2 ε ph /(k B T L ) -1 -1 , (2.5)
with T L the lattice temperature and k B the Boltzmann constant. Formally expanding the phonon collision operator in power of α 2 , we get

Q α ph (f ) = Q ph,0 (f ) + α 2 Q α ph,1 (f ) where Q α
ph,1 is of order 1 when α goes to 0. The electron-electron collision operator is given by [START_REF]Hot electon transport in semiconductors[END_REF]

Q e (f ) n (k) = n ′ ,r,s (R 2 ) 3 Φ e n,n ′ ,r,s (k, k ′ , k 1 , k ′ 1 )δ(ε n + ε n ′ ,1 -ε ′ r -ε ′ s,1 )δ(k + k 1 -k ′ -k ′ 1 ) [f ′ r f ′ s,1 (1 -ηf n )(1 -ηf n ′ ,1 ) -f n f n ′ ,1 (1 -ηf ′ r )(1 -ηf ′ s,1 )] dk ′ dk 1 dk ′ 1 .
(2.6)

The notation

f n ′ ,1 , f ′ r and f ′ s,1 stands for f n ′ (k 1 ), f r (k ′ ) and f s (k ′ 1 )
, respectively. We define then the elastic collision operator

Q 0 (f ) n = Q imp (f ) n + Q ph,0 (f ) n = n ′ ∈N * R 2 Φ 0 n,n ′ (k, k ′ )δ(ε n (k) -ε n ′ (k ′ ))(f n ′ (k ′ ) -f n (k)) dk ′ , (2.7)
where

Φ 0 n,n ′ = Φ imp n,n ′ + (2N ph + 1)Φ ph n,n ′ . We set Q α 1 (f ) = Q α ph,1 (f ) + 1 β Q e (f ). (2.8)
Then, the kinetic equation, starting point for the diffusive limits, is written in the following form

α 2 ∂ t f α n + α(∇ k ε n • ∇ x f α n -∇ x ε n • ∇ k f α n ) = Q 0 (f α ) n + α 2 Q α 1 (f α ) n .
(2.9)

Definitions and notations

We first recall the coarea formula : for any C 1 function g : B → R, and any test function ψ ∈ C 0 (B), we have :

B ψ(k) dk = R g -1 (ε) ψ(k) dS ε (k) |∇g(k)| dε,
where dS ε (k) denotes the Euclidian surface element on the manifold g -1 (ε). We denote

dN ε (k) = dS ε (k)/|∇g(k)|. Taking g(k) = |k| 2 2
/(2m * ), the set of possible wave vectors of electrons belonging to the nth subband and having total energy ε is given by S

ε-ǫn = {k ∈ R 2 s. t. |k| 2 = 2m * -2 (ε -ǫ n )} and dN ε-ǫn = dS ε-ǫn
|k| 2 /m * where dS ε-ǫn is the surface measure of the ball S ε-ǫn . The coarea formula leads to :

n∈N * R 2 ψ n (k) dk = n∈N * +∞ ǫn S ε-ǫn ψ n (k)dN ε-ǫn (k) dε, (2.10)
and

n∈N * R 2 ψ n (k)δ(ǫ n + |k| 2 2 2m * -ε) dk = n∈N * S ε-ǫn ψ n (k)dN ε-ǫn (k).
Definition 2.1 We will use the following notations :

• The density of states is defined by :

N(t, x, ε) := n∈N * S ε-ǫn dN ε-ǫn (k) = 2π m * 2 N (t, x, ε),
where N (t, x, ε) = max{n ∈ N * / ǫ n (t, x) ≤ ε}, with the convention of N (t, x, ε) = 0 if ε < ǫ 1 (t, x).

• The Fermi-Dirac function is given by

F µ,T (t, x, ε) = η + exp ε -µ k B T -1
.

• We introduce the two Hilbert spaces

L 2 := {f = (f n ) n∈N * , +∞ n=1 R 2 |f n (k)| 2 dk < +∞},
endowed with the natural scalar product

f, g = n∈N * R 2 f n (k)g n (k) dk, and
L 2 F = {f ∈ L 2 (R) s. t. R f 2 (ε) dε F (ε)(1 -ηF (ε))
< +∞}, endowed with the weighted scalar product defined by

f, g F = R f (ε)g(ε) dε F (ε)(1 -ηF (ε))
.

We will make the following assumption on the cross-section :

Assumption 2.2 The coefficient Φ 0 n,n ′ satisfies, for λ 0 and λ 1 two positive constants,'

0 < λ 0 < Φ 0 n,n ′ N(t, x, ε n ) < λ 1 , Φ 0 n,n ′ (k, k ′ ) = Φ 0 n ′ ,n (k ′ , k),
where N is the density of state defined above.

First macroscopic scaling : the spherical harmonic expansion model

For the sake of completness of this work, we present in this section the limit α → 0 of the kinetic equation (2.9). All calculations will be done formally and we refer the reader to [START_REF] Ben Abdallah | Adiabatic quantum-fluid transport models[END_REF] where the rigorous derivation is studied. We consider the Hilbert expansion

f α = f 0 + αf 1 + α 2 f 2 + • • •
By linearity of the operator Q 0 and by identifying the term of equal powers of α in (2.9), we obtain

Q 0 (f 0 ) = 0, (2.11)
Q 0 (f 1 ) n = 2 k m * • ∇ x f 0 n -∇ x ǫ n • ∇ k f 0 n , (2.12)
Q 0 (f 2 ) n = ∂ t f 0 n + 2 k m * • ∇ x f 1 n -∇ x ǫ n • ∇ k f 1 n -Q 0 1 (f 0 ) n , (2.13)
where Q 0 1 is obtained by taking α = 0 in the expression (2.8). We will then make use of the following properties of the collision operator (see [START_REF] Ben Abdallah | Adiabatic quantum-fluid transport models[END_REF]). Proposition 2.3 Under Assumption 2.2, the elastic collision operator Q 0 defined in (2.7) satisfies the following properties :

1. The linear operator Q 0 : L 2 → L 2 is a bounded, symmetric, non-positive operator.

For any bounded function

ψ : R → R, we denote ψ(ε) n (k) = ψ( |k| 2 2 2m * + ǫ n ). Then, ∀ f ∈ L 2 , Q 0 (ψ(ε)f ) = ψ(ε)Q 0 (f ). 3. The Kernel of Q 0 is the set Ker Q 0 = {f ∈ L 2 , s.t. ∃ ψ : R → R, f = ψ(ε)}.
4. The range R(Q 0 ) is closed and coincide with the orthogonal of the kernel of Q 0 given by :

(Ker Q 0 ) ⊥ = {f ∈ L 2 , s.t. n∈N * S ε-ǫn f n (k) dN ε-ǫn (k) = 0, for a.e. ε ≥ ǫ 1 }.
From Proposition 2.3 and (2.11), we deduce that f 0 is an energy dependent function :

f 0 n (t, x, k) = F (t, x, ε n ).
Choosing ψ : R → R such that kψ(ε) ∈ L 2 , we deduce from Proposition 2.3 that there exists a unique solution ξ in (Ker Q 0 ) ⊥ such that

-Q 0 (ξ) = 2 k m * ψ(ε).
We can write ξ = Θ • ψ(ε) and from the second item of Proposition 2.3, Θ is independent of the choice of the function ψ. Then, equation (2.12) leads to

f 1 n (t, x, k) = -Θ n (t, x, k) • ∇ x F (t, x, ε n ). (2.14)
Finally, the solvability condition of equation (2.13) is that the right hand side belongs to (Ker Q 0 ) ⊥ . This leads to

n∈N * S ε-ǫn (∂ t f 0 n + ∇ k ε n • ∇ x f 1 n -∇ x ǫ n • ∇ k f 1 n -Q 0 1 (f 0 ) n ) dN ε-ǫn (k) = 0, for a.e. ε ≥ ǫ 1 .
(2.15) Let us denote

S e (F ) = n∈N * S ε-ǫn Q e (F ) n dN ε-ǫn (k), and S 1 (F ) = n∈N * S ε-ǫn Q 0 ph,1 (F ) n dN ε-ǫn (k).
(2.16) Multiplying (2.15) by an energy-dependent test function φ(ε) and integrating with respect to the variable ε, we obtain for the first term :

∞ ǫ 1 n∈N * S ε-ǫn ∂ t f 0 n dN ε-ǫn (k)φ(ε) dε = ∞ ǫ 1 n∈N * S ε-ǫn (∂ t F + ∂ ε F ∂ t ǫ n ) dN ε-ǫn (k)φ(ε) dε = ∞ ǫ 1 N∂ t F φ(ε) dε + ∞ ǫ 1 ∂ ε F n∈N * ∂ t ǫ n S ε-ǫn dN ε-ǫn (k) φ(ε) dε.
Using the coarea formula (2.10), we deduce that

∞ ǫ 1 n∈N * S ε-ǫn (∇ k ε n • ∇ x f 1 n -∇ x ǫ n • ∇ k f 1 n ) dN ε-ǫn (k) φ(ε) dε = = n∈N * R 2 (∇ x • ( 2 k m * f 1 n ) -∇ k • (f 1 n ∇ x ǫ n ))φ(ε n ) dk = ∇ x • n∈N * R 2 2 k m * f 1 n φ(ε n ) dk = = - ∞ ǫ 1 ∇ x • n∈N * S ε-ǫn 2 k m * ⊗ Θ n dN ε-ǫn (k) • ∇ x F φ(ε) dε,
where the last identity is a consequence of (2.14). We define the diffusion matrix by

D(t, x, ε) := n∈N * S ε-ǫn 2 k m * ⊗ Θ n dN ε-ǫn (k), (2.17)
and the current density by

J(t, x, ε) = -D(t, x, ε) • ∇ x F (t, x, ε). (2.18)
With these notations, we get that in the distributional sense, equation (2.15) is equivalent to the spherical harmonic expansion (SHE) model

N∂ t F + ∇ x • J -κ∂ ε F = 1 β S e (F ) + S 1 (F ), (2.19)
where κ is given by

κ(t, x, ε) = -2π m * 2 ∂ t n∈N * (ε -ǫ n ) + .
(2.20)

The notation u + = max{0, u} denotes the positive part of u. We recall moreover a property of the diffusion matrix D stated in Lemma 2.8 of [START_REF] Ben Abdallah | Adiabatic quantum-fluid transport models[END_REF]. We point out that the effect of the confinement is reflected in the special form of the coefficients of (2.19), which involve the subband energies.

Lemma 2.4

The diffusion matrix D(t, x, ε) define in (2.17) is a symmetric and nonnegative 2 × 2 matrix.

Remark 2.5 We end this section with a particular choice of the cross-section Φ 0 n,n ′ , which allows to compute explicitly the diffusion matrix. If the cross-section is an energy-dependent function of the form

Φ 0 n,n ′ (t, x, k, k ′ ) = Φ 0 (t, x, ε n )
, then, after a straightforward computation, we have that

Θ n (t, x, k) = 1 Φ 0 (t, x, ε n )N(t, x, ε n ) ∇ k ε n .
Therefore, the diffusion matrix defined in (2.17) has the expression

D(t, x, ε) = π 2 /m * Φ 0 (t, x, ε)N(t, x, ε) n∈N * (ε -ǫ n ) + Id.
(2.21)

Second macroscopic scaling : the energy-transport model

We start from the SHE model (2.19) and we assume that the electron-electron collision operator is dominant with respect to the second order correction of the phonon collision operator and therefore β ≪ 1, in order to obtain an ET model. Passing through the SHE model, instead of starting directly from the Boltzmann equation allows to get an explicit expression of the coefficients, which is needed for numerical purpose. Moreover, the considered dominant scattering mechanisms provide an energy relaxation term in the macroscopic limiting model.

The formal limit β → 0 in (2. [START_REF] Gadau | A three-dimensional mixed finite-element approximation of the semiconductor energy-transport equations[END_REF]) is again performed by means of a Hilbert expansion

F = F 0 + βF 1 + • • • Identifying equal powers of β implies S e (F 0 ) = 0, (2.22)
N∂ t F 0 + ∇ x J 0 -κ∂ ε F 0 -S 1 (F 0 ) = D F 0 S e (F 1 ), (2.23)
where D F 0 S e denotes the Fréchet derivative of S e at F 0 . We summarize below some useful properties of the collision operator S e defined in (2.16) and of its Fréchet derivative.

Proposition 2.6 Under micro-reversibility assumptions on the cross-section Φ e , the operator S e satisfies the following properties :

(i) For all f, g ∈ L 2 (R), we have :

R S e (f )(ε)g(ε)dε = = - 1 4 n,n ′ ,r,s (R 2 ) 4 Φ e n,n ′ ,r,s δ ε δ k [f (ε r (k ′ ))f (ε s (k ′ 1 ))(1 -ηf (ε n (k)))(1 -ηf (ε n ′ (k 1 ))) -f (ε n (k))f (ε n ′ (k 1 ))(1 -ηf (ε r (k ′ )))(1 -ηf (ε s (k ′ 1 )))[g ′ r + g ′ s,1 -g n -g n ′ ,1 ] dkdk 1 dk ′ dk ′ 1 .
(ii) kernel :

Ker S e = {f ∈ L 2 (R); ∃ µ(t, x), T (t, x) s. t. f (t, x, ε) = F µ,T (t, x, ε)},
where F µ,T is the so-called Fermi-Dirac distribution function (see Definition 2.1).

Proposition 2.7

The linear operator D F S e satisfies (i) D F S e is bounded, symmetric, non-positive on L 2 F . (ii) The kernel of D F S e is given by

Ker (D F S e ) = Span{F (1 -ηF ), F (1 -ηF )ε} (iii) The range of D F S e is closed and we have R(D F S e ) = Ker (D F S e ) ⊥ = f ∈ L 2 F / R f (ε) 1 ε dε = 0 .
These properties are an easy consequence of Proposition 3.16, Proposition 3.17 and Proposition 3.19 of [START_REF] Ben Abdallah | Adiabatic quantum-fluid transport models[END_REF], using the fact that :

R S e (f )(ε)g(ε) dε = +∞ ǫ 1 n∈N * S ε-ǫn Q e (f ) n (ε)g(ε) dN ε-ǫn (k)dε = n∈N * +∞ ǫn S ε-ǫn Q e (f ) n (ε)g(ε) dN ε-ǫn (k)dε = n∈N * R 2 Q e (f ) n (ε n )g(ε n ) dk,
where we use the coarea formula (2.10) for the last identity.

Formal derivation of energy-transport model. Let us come back to the formal limit β → 0 in (2.19). Thanks to Proposition 2.6, equation (2.22) implies that there exist µ(t, x) and T (t, x) such that

F 0 (t, x, ε) = F µ,T (t, x, ε) (2.24)
From Proposition 2.7 we deduce that equation (2.23) admits a solution iff

R (N∂ t F 0 + ∇ x • J 0 -κ∂ ε F 0 -S 1 (F 0 )) 1 ε dε = 0. (2.25)
For the first term, the definition of the density of states N (see Definition 2.1) implies that

N(t, x, ε) = 2π m * 2 n if ε ∈ [ǫ n , ǫ n+1 ) and vanishes for ε < ǫ 1 . Then R N∂ t F 0 1 ε dε = n∈N * ǫ n+1 ǫn 2π m * 2 n∂ t F 0 1 ε dε = ∂ t R NF 0 1 ε dε -2π m * 2 n∈N * ∂ t ǫ n F 0 (ǫ n ) 1 ǫ n .
Using the expression of the current (2.18), we can rewrite the second term of (2.25) :

R ∇ x • J 0 1 ε dε = -∇ x • R D(t, x, ε) • ∇ x F 0 (t, x, ε) 1 ε dε .
From (2.20), we deduce

R κ∂ ε F 0 1 ε dε = n∈N * ∞ ǫn 2π m * 2 ∂ t ǫ n ∂ ε F 0 1 ε dε = - n∈N * 2π m * 2 ∂ t ǫ n F 0 (ǫ n ) 1 ǫ n - n∈N * 2π m * 2 ∂ t ǫ n ∞ ǫn F 0 (ε) 0 1 dε,
where we use an integration by part for the last identity. Finally, the solvability condition (2.25) writes in the following form :

∂ t R NF 0 1 ε dε + n∈N * 2π m * 2 ∂ t ǫ n ∞ ǫn F 0 (ε) 0 1 dε -∇ x • R D(t, x, ε) • ∇ x F 0 (t, x, ε) 1 ε dε = R S 1 (F 0 ) 1 ε dε.
(2.26)

Let us denote by ρ and ρE the charge density and the energy density, respectively, associated to the Fermi-Dirac distribution function F µ,T :

ρ µ,T (t, x) := R NF µ,T (t, x, ε) dε = n∈N * R 2 F µ,T (t, x, ε n ) dk, (2.27) ρE µ,T (t, x) := R NF µ,T (t, x, ε)ε dε = n∈N * R 2 ε n F µ,T (t, x, ε n ) dk. (2.28)
We can state easily that for a Fermi-Dirac function, we have

∇ x F µ,T (t, x, ε) = -F µ,T (1 -ηF µ,T ) ε∇ x ( 1 k B T ) -∇ x ( µ k B T ) .
Then, equation (2.26) reads, using R S 1 (F ) dε = 0,

∂ t ρ µ,T -∇ x • J 1 = 0, (2.29)
∂ t (ρE µ,T ) + n∈N * 2π m * 2 ∂ t ǫ n ∞ ǫn F µ,T (t, x, ε) dε -∇ x • J 2 = W, (2.30)
where we denote

J 1 := R D(t, x, ε) • ∇ x F µ,T (t, x, ε) dε = D 00 ∇ x ( µ k B T ) -D 01 ∇ x ( 1 k B T ), (2.31)
J 2 := R ε D(t, x, ε) • ∇ x F µ,T (t, x, ε) dε = D 10 ∇ x ( µ k B T ) -D 11 ∇ x ( 1 k B T ), (2.32)
and where the diffusion coefficients are defined by

D ij (t, x) = R D(t, x, ε)ε i+j F µ,T (1 -ηF µ,T ) dε, for i, j = 0, 1, (2.33)
with D being defined in (2.17). The system (2.29)-(2.33) forms the energy-transport model in the transport direction for a partially confined electron gas. We recover the general form for energy-transport models for semiconductors (see [START_REF] Ben Abdallah | On a hierarchy of macroscopic models for semiconductors[END_REF][START_REF] Ben Abdallah | An energy transport model for the semiconductors derived from the Boltzmann equation[END_REF][START_REF] Jüngel | Transport Equations for Semiconductors[END_REF] and references therein). The right hand side of the energy equation (2.30) is the so-called relaxation term W and it is defined by

W := R εS 1 (F µ,T ) dε = n∈N * +∞ ǫn S ε-ǫn Q 0 ph,1 (f ) n dN ε-ǫn (k) ε dε. (2.34)
As for the SHE model (2.19), the presence of the subband energies in the diffusion matrix and in the relaxation term reflects the effect of the confinement in the transport equation.

Then we have the following important property for the diffusion matrix which is an easy consequence of expression (2.33) and Lemma 2.4.

Lemma 2.8 The diffusion matrix

D := D 00 D 01 D 10 D 11
is a symmetric and positive definite matrix.

Remark 2.9 After straightforward calculations, we can have an explicit expression of ρ and ρE. In fact,

ρ µ,T (t, x) = n∈N * ǫ n+1 ǫn 2π m * 2 nF µ,T (t, x, ε) dε = 2πm * k B T η 2 n∈N * log(1 + η exp( µ -ǫ n k B T )), ρE µ,T (t, x) = n∈N * 2π m * 2 n ǫ n+1 ǫn ε F µ,T (t, x, ε) dε if η > 0.
For Boltzmann statistics η = 0, we have

ρ µ,T (t, x) = n∈N * 2π m * 2 k B T (t, x) exp( µ(t, x) -ǫ n (t, x) k B T (t, x) ), (2.35
)

ρE µ,T (t, x) = n∈N * 2π m * 2 k B T (t, x)(ǫ n (t, x) + T (t, x)) exp( µ(t, x) -ǫ n (t, x) k B T (t, x)
).

(2.36)

Relaxation term W

The relaxation term defined in (2.34) measures the influence of the interaction of phonons with the charge carriers. In [START_REF] Pietra | Modeling and simulation of the diffusive transport in a nanoscale Double-Gate MOSFET[END_REF] we have formally established that a diffusion limit of the kinetic Boltzmann transport equation coupled to subband model in the scaling of dominant phonon-electron interaction leads to a drift-diffusion system in the transport direction coupled to the subband model. Using the coarea formula we can rewrite (2.34) as

W = n∈N * R 2 ε n Q 0 ph,1 (F µ,T ) dk.
Moreover, we have

n∈N * R 2 ε n Q α ph (F µ,T ) n dk = n,n ′ R 4 Φ ph n,n ′ (k, k ′ ) ([(N ph + 1)δ(ε n -ε ′ n ′ + α 2 ε ph ) + N ph δ(ε n -ε ′ n ′ -α 2 ε ph )]ε n F ′ n ′ (1 -ηF n )- [(N ph + 1)δ(ε ′ n ′ -ε n + α 2 ε ph ) + N ph δ(ε ′ n ′ -ε n -α 2 ε ph )]ε n F n (1 -ηF ′ n ′ )) dkdk ′ ,
where the notation ε ′ n ′ , F n and F ′ n ′ stands, respectively, for

ε n ′ (k ′ ), F µ,T (ε n ) and F µ,T (ε ′ n ′ ). Since the Fermi-Dirac distribution function is energy-dependent, we have n∈N * R 2 ε n Q α ph (F µ,T ) n dk = n,n ′ R 4 Φ ph n,n ′ (k, k ′ )(1 -ηF ′ n ′ )(1 -ηF n ) ([(N ph + 1)δ(ε n -ε ′ n ′ + α 2 ε ph ) + N ph δ(ε n -ε ′ n ′ -α 2 ε ph )]ε n M ′ n ′ - [(N ph + 1)δ(ε ′ n ′ -ε n + α 2 ε ph ) + N ph δ(ε ′ n ′ -ε n -α 2 ε ph )]ε n M n ) dkdk ′ ,
where M n = e (µ-εn)/(k B T ) is the Maxwellian and it satisfies

F n = M n (1 -ηF n ). Moreover, n∈N * R 2 ε n Q α ph (F µ,T ) n dk = n,n ′ R 4 Φ ph n,n ′ (k, k ′ )(1 -ηF ′ n ′ )(1 -ηF n ) (δ(ε n -ε ′ n ′ + α 2 ε ph )[(N ph + 1)M ′ n ′ -N ph M n ]ε n - δ(ε ′ n ′ -ε n + α 2 ε ph )[(N ph + 1)M n -N ph M ′ n ′ ](ε ′ n ′ + α 2 ε ph )) dkdk ′ .
By changing the variable ε n with ε ′ n ′ in the first term of the sum, we notice that only the α 2 -factor term does not vanish in this last identity, which can be rewritten as

n,n ′ R 4 Φ ph n,n ′ (1 -ηF ′ n ′ )(1 -ηF n )ε ph δ(ε n -ε ′ n ′ -α 2 ε ph )[N ph M ′ n ′ -(N ph + 1)M n ] dkdk ′ = n,n ′ R 4 Φ ph n,n ′ (1 -ηF ′ n ′ )(1 -ηF n )ε ph δ(ε n -ε ′ n ′ -α 2 ε ph )N ph M n (e α 2 ε ph k B T -e α 2 ε ph k B T L ) dkdk ′ ,
where the phonon occupation number N ph is defined in (2.5). Letting α → 0, we have that

N ph (e α 2 ε ph /k B T -e α 2 ε ph /k B T L ) → T L ( 1 T -1 T L ). Thus W = T L ( 1 T - 1 T L ) n,m R 4 Φ ph n,n ′ (k, k ′ )F n (1 -ηF n )ε ph δ(ε n -ε ′ n ′ ) dkdk ′ .
(2.37)

The following lemma proves that W is a temperature relaxation term which relaxes T to the lattice temperature T L . Lemma 2.10 Let W be defined in (2.34). Then, we have

W • (T -T L ) ≤ 0.
Proof. The proof of this result is an immediate consequence of (2.37).

Formal derivation of drift-diffusion equation.

In the case where the electron-phonon scattering is dominant, after a rescaling we have W = 1 γ W with a parameter γ ≪ 1. Then equation (2.30) with expression (2.37) implies that formally in the limit γ → 0, we have T = T L . Then (2.29) leads to the well-known drift-diffusion model

∂ t ρ µ -∇ x D 00 k B T L ∇ x µ = 0. (2.38)
Moreover, assuming η = 0, the equilibrium is then given by a Boltzmann statistics. We deduce therefore from (2.35) that for

T = T L ρ µ = 2π m * 2 k B T L e µ/k B T L n∈N * e -ǫn/k B T L . (2.39) Then ∇ x ρ µ = 1 k B T L ρ µ ∇ x µ + ρ µ ∇ x ( n e -ǫn/k B T L ) n e -ǫn/k B T L .
We can introduce as in [START_REF] Ben Abdallah | Diffusive transport of partially quantized particles : Existence, uniqueness and long time behaviour[END_REF] the effective potential energy defined by

V s = -k B T L log n∈N * e -ǫn/k B T L . (2.40)
Then, we have

ρ µ ∇ x µ = (k B T L ∇ x ρ µ + ρ µ ∇ x V s ) ,
such that, denoting D = D 00 /ρ µ , we recover from (2.38) the standard formulation of the drift-diffusion system [START_REF] Ben Abdallah | Diffusive transport of partially quantized particles : Existence, uniqueness and long time behaviour[END_REF] :

∂ t ρ µ -∇ x • (D(k B T L ∇ x ρ µ + ρ µ ∇ x V s )) = 0.
(2.41) Remark 2.11 We conclude the section noticing the similarities and the differences between the classical and the partially confined energy-transport model. Due to partial confinement in the z direction, the electron density in the transport direction contains T as factor rather than T 3/2 . Moreover, we point out that the system (2.29)-(2.33) in the variables µ/k b T , -1/k b T is in symmetric form, with the electric forces appearing in the diffusion coefficients through the eigenenergies ǫ n (which in turns depend on the electrostatic potential).

In the partially confined framework, the diffusive limiting process brought directly to a set of variables which can be interpreted as dual entropy variables (using the denomination of nonequilibrium thermodynamics [START_REF] De Groot | Nonequilibrium thermodynamics[END_REF][START_REF] Kreuzer | Nonequilibrium thermodynamics and its statistical fundation[END_REF]). In this framework, T can be interpreted as electron temperature, and the variable µ as quasi-Fermi potential energy. This fact is clear considering the drift-diffusion model (2.41). Indeed, using the effective potential energy V s defined in (2.40), we can write the electron density (2.39) as

ρ µ = 2π m * 2 k B T L e (µ-Vs)/k B T L ,
which gives the classical relation between electron density and quasi-Fermi energy. A chemical potential can then be defined as µ chem = µ -V s .

Numerical simulation

The device we are modelling in this work is a nanoscale Double-Gate MOSFET (Metal Oxide Semiconductor Field Effect Transistor) such as in [START_REF] Pietra | Modeling and simulation of the diffusive transport in a nanoscale Double-Gate MOSFET[END_REF]. This device consists of a silicon film, characterized by two highly doped regions near the Ohmic contacts (denoted by source and drain) and an active region, called channel, with lower doping. The silicon film is sandwiched between two thin layers of silicon dioxide SiO 2 , each of them with a gate contact. We assume invariance in the x 2 direction (infinite boundary conditions), so that the problem is studied in a (x 1 , z)-domain. The device occupies a region of a 2D domain denoted by Ω=[0, L] × [0, ℓ]. A schematic representation of the device is shown in Figure 1.

Energy-transport-Schrödinger-Poisson system

In the following, we describe the collisional transport in the Double-Gate MOSFET, schematized in Figure 1, by means of the energy-transport model (2.29)-(2.33). The confinement is described by the subband decomposition approach, which involves the resolution of the eigenvalue problem (2.1), taking also into account the presence of the oxide. Moreover, in order to provide explicitly computable diffusion and relaxation terms, the following physical assumptions will be used. • (H1) The cross-sections Φ 0 and Φ ph are assumed to be energy dependent functions and to have the following expression (see [START_REF] Ben Abdallah | On a hierarchy of macroscopic models for semiconductors[END_REF][START_REF] Ben Abdallah | High field approximation of the spherical harmonic expansion model for semiconductors[END_REF][START_REF] Degond | Numerical discretization of energy-transport models for semiconductors with non-parabolic band structure[END_REF])

Φ 0 n,n ′ (t, x, k, k ′ ) = φ 0 (t, x)ε s n , Φ ph n,n ′ (t, x, k, k ′ ) = φ ph (t, x)ε s n , (3.1)
with -2 < s < 2. In the physical literature, the values s = 0 and s = 1/2 have been used [START_REF] Chen | An improved energy transport model including nonparabolicity and non-maxwellian distribution effects[END_REF][START_REF] Jüngel | Transport Equations for Semiconductors[END_REF][START_REF] Lyumkis | Transient semiconductor device simulation including energy balance equation[END_REF].

• (H2) The electron density and the energy are assumed to be given by non-degenerate Boltzmann statistics, i.e. η = 0, as in (2.35), (2.36).

Using Assumption (H1), we deduce that the diffusion matrix (2.21) has the following expression

D(t, x, ε) = 1 φ 0 (t, x)ε s N (ε) n∈N * (ε -ǫ n ) + Id, (3.2)
where N (ε) = max{n ∈ N * / ǫ n ≤ ε} is the number of non-zero terms in the sum (see Definition 2.1). Moreover, thanks to the coarea formula in (2.37), we have

W = -4π 2 ε ph (1 - T L T ) +∞ ǫ 1 φ ph ε s N 2 (ε)F (ε)(1 -ηF (ε)) dε,
which under Assumption (H2), reads

W = -4π 2 ε ph (1 - T L T )e µ/(k B T ) +∞ ǫ 1 φ ph ε s N 2 (ε)e -ε/(k B T ) dε.
By defining

W 0 = 4π 2 ε ph +∞ ǫ 1 φ ph ε s N 2 (ε)e -ε/(k B T ) dε, (3.3)
we have the compact expression

W = -W 0 1 - T L T e µ/(k B T ) . (3.4)
Assumption (H2) implies also that the density of charge carriers N e (t, x, z) is given by

N e = n∈N * R 2 F µ,T (t, x, |k| 2 2 2m * + ǫ n ) dk|χ n | 2 = 2πm * k B T 2 n∈N * e (µ-ǫn)/k B T |χ n | 2 .
Finally, the coupled subband energy-transport model under Assumption 3.1 is given by : Find µ(t, x), T (t, x), (ǫ n (t, x), χ n (t, x)) for n ≥ 1, and V (t, x, z) such that

∂ t ρ µ,T -∇ x • J 1 = 0, in (0, L) (3.5) ∂ t (ρE µ,T ) + n∈N * 2π m * 2 ∂ t ǫ n k B T (t, x)e (µ-ǫn)/k B T -∇ x • J 2 = W, in (0, L) (3.6)        - 2 2 d dz 1 m * (z) d dz χ n -e(V + V c )χ n = ǫ n χ n , in (0, ℓ) χ n (t, x, •) ∈ H 1 0 (0, ℓ), ℓ 0 χ n χ n ′ dz = δ nn ′ , (3.7) div x,z (ε R ∇ x,z V ) = e ε 0 2πm * k B T 2 n∈N * e (µ-ǫn)/k B T |χ n | 2 -N D , in Ω, (3.8)
where the expressions of ρ µ,T and E µ,T with respect to the unknows are given in (2.35)-(2.36). In (3.7) the effective mass m * takes different values in the Si and in the SiO 2 domain. Moreover, V c represents a given potential barrier between the silicon and the oxide. The currents J 1 and J 2 are given by the expressions

J 1 = D 00 ∇ x µ k B T -D 01 ∇ x 1 k B T , (3.9
)

J 2 = D 10 ∇ x µ k B T -D 11 ∇ x 1 k B T , (3.10)
where, under Assumption 3.1, the diffusion coefficients are given by

D ij = 1 φ 0 n∈N * +∞ ǫn ε i+j-s (ε -ǫ n ) N (ε) e (µ-ε)/k B T dε.
The relaxation term W is given by (3.4). This system is complemented with initial and boundary conditions. In particular, at the ohmic contacts and at the gate (see Figure 1), we will impose Dirichlet boundary conditions for the potential, otherwise we fix homogeneous Neumann boundary conditions, which model isolating conditions.

V (x, z) = V Gate , for z ∈ {0, ℓ}, x ∈ Gate ; (3.11) V (x, z) = V D , for x ∈ {0, L}, z ∈ (0, ℓ); (3.12) ∂V ∂ν = 0, elsewhere, (3.13)
where ν is the outward unit normal. Since the transport occurs only in the longitudinal direction, we just have to impose boundary conditions in x = 0 and x = L for µ and T . The temperature is assumed to be at the lattice temperature T L , thus

T (x) = T L , for x ∈ {0, L}. (3.14)
Then, we consider that the surface density of the charge carriers is almost constant near the frontiers x = 0 and x = L and given by N b s . The surface density being the integral over z of the total density (N b s = N + × ℓ Si ), we deduce

µ(x) = µ b := k B T L log N b s 2 2πmk B T L n e -ǫn/k B T L , for x ∈ {0, L}. (3.15)

Stationary system

Let us introduce the notations

u = µ k B T , v = - 1 k B T . (3.16)
Then we can rewrite the expressions of the current (3.9)-(3.10) as

J 1 = D 00 (u, v)∇ x u + D 01 (u, v)∇ x v, (3.17)
J 2 = D 10 (u, v)∇ x u + D 11 (u, v)∇ x v, (3.18)
where the diffusion coefficients are given by

D ij (u, v) = 1 φ 0 n∈N * +∞ ǫn ε i+j-s (ε -ǫ n ) N (ε) e u+εv dε. (3.19)
We define the relaxation coefficient in the same way :

W 0 (u, v) = 4π 2 ε ph +∞ ǫ 1 φ ph ε s N 2 (ε)e (ε-µ)v dε. (3.20)
Then, the stationary version of the energy-transport subband system (3.5)-(3.6) in variable u and v writes

-∇ x (D 00 (u, v)∇ x u + D 01 (u, v)∇ x v) = 0, (3.21)
-∇ x (D 10 (u, v)∇ x u + D 11 (u, v)∇ x v) = -W 0 (u, v) (1 + k B T L v) , (3.22)
Boundary conditions (3.14)-(3.15) become

v(x) = v b := - 1 k B T L , for x ∈ {0, L}, and
u(x) = u b := log N b s 2 2πmk B T L n e -ǫn/k B T L , for x ∈ {0, L}.

Numerical approach for the energy-transport system

We introduce a partition of [0, L] with nodes x i , i = 0, • • • , N x , and a partition of [0, ℓ] with nodes z j , j = 0, • • • , N z . We assume that the partitions are uniform and denote h = x ix i-1 . Then, we mesh the domain [0, L] × [0, ℓ] with rectangular triangles using the nodes (x i , z j ) previously defined. The Schrödinger equations and the Poisson equation are discretized with conforming P 1 finite elements.

We consider here in details the discretization scheme for the equations governing u and v, assuming first that the eigenenergies ǫ n are known. Using the following notations

U = (u, v) ⊤ , J = (J 1 , J 2 ) ⊤ , W (U) = (0, W 0 (u, v)(1 + k B T L v)) ⊤
equations (3.17

J = D(U)∇ x U, ∇ x • J = W (U). (3.23)
Denoting by U i an approximation of U(x i ), we take the piecewise constant approximation of U given, in the interval

I i := (x i-1 , x i ), by U i = U i-1 + U i 2 .
and define the piecewise constant diffusion coefficients and relaxation terms as

D kl = D kl (U ), for k, l = 0, 1; W = W (U). (3.24)
We are going to use a mixed finite element discretization of lowest order in hybridized form [START_REF] Brezzi | Mixed and hybrid finite element methods[END_REF] (see also [START_REF] Gadau | A three-dimensional mixed finite-element approximation of the semiconductor energy-transport equations[END_REF][START_REF] Gadau | A mixed finite-element scheme of a semiconductor energy-transport model using dual entropy variables, Hyperbolic problems: theory, numerics and applications[END_REF] for applications to ET). Let us introduce the following finite dimensional spaces :

X h = {σ ∈ L 2 (Ω) : σ(x) = a i + b i (x -x B i ) in I i , i = 1, . . . , N x }, Y h = {ξ ∈ L 2 (Ω) : ξ is constant in I i , i = 1, . . . , N x }, Z h,χ =
{q is defined at the nodes x 0 , . . . , x Nx , q(x 0 ) = χ(0), q(x Nx ) = χ(1)}, where x B i denotes the central point of the interval I i , and χ is prescribed.

Then, the mixed-hybrid formulation of (3.23) reads as follows: Find Thanks to the discontinuity of the spaces X h and Y h , we can apply static condensation in order to reduce the size of the discrete system and obtain an algebraic system for the variable U h only. More precisely, choosing first the local basis

J h ∈ X 2 h , P h ∈ Y 2 h , and U h ∈ Z h,u b × Z h,v b such that Nx i=1 I i φ h • D(U i ) -1 J h dx + I i P h • ∇ x φ h dx -[φ h • U h ] x i x i-1 = 0, (3.25)
Nx i=1 I i Ψ h • ∇ x J h dx - I i W (U i ) • Ψ h dx = 0, (3.26)
- Nx i=1 [µ h • J h ] x i x i-1 = 0 (3.27) for all φ h ∈ X 2 h , Ψ h ∈ Y 2 h ,
φ h =
(1, 0) ⊤ in I i (0, 0) ⊤ elsewhere , φ h = (0, 1) ⊤ in I i (0, 0) ⊤ elsewhere in (3.25) and then, analogously choosing Ψ h in (3.26), we obtain the piecewise linear (discrete) current

J h|I i = D(U i ) U i -U i-1 h + W (U i )(x -x B i ). (3.28)
Imposing continuity at the nodes (through (3.27)) we obtain the final system

-D(U i )U i-1 + (D(U i ) + D(U i+1 ))U i -D(U i+1 )U i+1 = - h 2 2 (W (U i ) + W (U i+1 )), (3.29)
for i = 1, . . . , N x -1. We point out explicitly that, since the first component of W (U) is null, the approximation of the current J 1 is piecewise constant (see (3.28)) and that, thanks to (3.27), it is indeed globally constant. System (3.29) forms a non-linear system in the unknown (u, v) that can be solved using a Newton algorithm. We point out that the Jacobian corresponding to this non-linear system can be easily computed noticing that, from the expressions (3.24) and (3.19), we have

∂D i kℓ ∂u i = ∂D i kℓ ∂u i-1 = 1 2 D i kℓ , ∂D i kℓ ∂v i = ∂D i kℓ ∂v i-1 = 1 2 D i kℓ+1 ,
and that similar relations hold for the partial derivatives of W i .

Remark 3.2 The expression (3.19) is not practical for numerical purpose. However, with the definition N (ε) = max{n ∈ N * : ǫ n ≤ ε} , we have

D ij (u, v) = 1 φ 0 +∞ n=1 +∞ m=n ǫ m+1 ǫm (ε -ǫ n )ε i+j-s m e u+εv dε = 1 φ 0 +∞ m=1 ǫ m+1 ǫm (ε - m n=1 ǫ n m )ε i+j-s e u+εv dε,
by interchanging the sums over m and n. We can rewrite

D ij (u, v) = 1 φ 0 +∞ ǫ 1 ε i+j+1-s e u+εv dε + 1 φ 0 +∞ m=1 m-1 n=1 ǫ n m -1 - m n=1 ǫ n m +∞ ǫm ε i+j-s e u+εv dε,
with the convention that m-1 n=1 ǫ n /(m -1) = 0 for m = 1. Then, in the actual numerical computation, we can get an accurate approximation of D ij by truncating the infinite sum to a finite number of eigenmodes. In fact, (ǫ n ) n forms an increasing sequence going to +∞, thus, since v < 0, we have that (e ǫnv ) n fast decreases to 0.

Algorithm

We are now ready to describe the algorithm used for the numerical resolution of the stationary subband energy-transport model. The first step of the algorithm is the computation of the thermal equilibrium solution, with no applied drain-source bias. In this case the temperature and the Fermi level are constant along the device, therefore the problem reduces to solving the Poisson equation (3.8) for a given temperature and Fermi level computed thanks to the boundary conditions. The computed potential at thermal equilibrium is used as a starting data for the following Gummel [START_REF] Gummel | A self-consistent iterative scheme for one-dimensional steady state transistor calculations[END_REF] iteration process :

1. Lets V old be a given potential.

2. We solve the eigenvalue problem (3.7) on each vertical slice of the domain by diagonalization of the Hamiltonian. Therefore we obtain the set {χ n (x i , z j )} and {ǫ n (x i )}.

3. We implement the Newton procedure which has been described above for the computation of (u, v).

4. We compute the density of charge carriers corresponding to the right hand side of (3. 8)

N e = n∈N * 2π m * 2 e u v e ǫnv |χ n | 2 .
We are then able to solve the Poisson equation (3.8) with boundary condition (3.11). Indeed, due to the strong coupling of the entire system, the simple resolution of equation (3.8) does not provide a converging algorithm. Following [START_REF] Ph | Finite element approximation of electrostatic potential in one dimensional multilayer structures with quantized electronic charge[END_REF] a Gummel iteration algorithm is used, amounting to compute the new potential V new by solving the following modified Poisson equation

div(ε R ∇V new ) + e ε 0 N e V new V ref = e ε 0 N e (1 - V old V ref ) -N D , in Ω, (3.30)
with

V ref = k B T L /e.
5. We repeat the last three steps until the quantity V old -V new L ∞ becomes sufficiently small. Once the convergence is reached, we increment the applied drain-source bias V DS of 0.02 V and start a new iteration.

Numerical results

In this section we present and comment the performed numerical results. The modeled device is schematized in Figure 1 and the physical values are chosen as the ones in [START_REF] Pietra | Modeling and simulation of the diffusive transport in a nanoscale Double-Gate MOSFET[END_REF] which are recalled in Table 1. We take N x = 50 points in the transport direction and N z = 50 in the confined direction for all the tests. The results are presented for s = 1 2 in (3.1), which corresponds to the so-called Chen model [START_REF] Chen | An improved energy transport model including nonparabolicity and non-maxwellian distribution effects[END_REF].

Some other physical coefficients should be determined. The effective mass is m * = 0.19 m e with m e the electron mass, the lattice temperature is T L = 300 K and the scattering coefficient for the elastic collisions is the one used in [START_REF] Pietra | Modeling and simulation of the diffusive transport in a nanoscale Double-Gate MOSFET[END_REF] φ 0 = 1 µ 0 n i , where the low field

m -3 L S 10nm N - 10 21 m -3 L C 30nm U c 3 eV L D 10nm ε R [Si] 11.7 ℓ ox 3nm ε R [SiO 2 ]
3.9 ℓ Si 5nm mobility is taken as µ 0 = 0.12 m 2 V -1 s -1 and the intrinsic density is given by n i = 10 10 m -2 .

We have to fix the value of the scattering coefficient φ ph for the electron-phonon interaction.

As noticed in subsection 2.6, when φ ph → +∞, the model converges formally to the subband drift-diffusion system presented in [START_REF] Ben Abdallah | Diffusive transport of partially quantized particles : Existence, uniqueness and long time behaviour[END_REF][START_REF] Pietra | Modeling and simulation of the diffusive transport in a nanoscale Double-Gate MOSFET[END_REF]. It is then interesting to compare the numerical results for large and small values of φ ph . Figure 2 displays the I -V DS characteristics for φ ph = 10 -4 /φ 0 and φ ph = 10 5 /φ 0 . As expected and as noticed in [START_REF] Baccarani | A compact double-gate MOSFET model comprising quantum-mechanical and nonstatic effects[END_REF], the energytransport model gives higher currents compared with the drift-diffusion model, due to the electron velocity overshoot within the channel. Figure 2 (right) shows the temperature for φ ph = 10 5 /φ 0 , confirming that we are in the drift-diffusion regime, where the temperature is constant. Indeed we check numerically that max Tmin T = 2.6584 10 -7 K. In the rest of the section, we present the results only for φ ph = 10 -4 /φ 0 which corresponds to the energy-transport regime. Figure 3 displays the computed current vs drainsource applied bias characteristics with this chosen value. We present in Figure 3 (left) the characteristics for different numerical values of the Gate voltage V G and with ℓ Si = 5 nm. In Figure 3 (right), we display the characteristics for V G = 0 V and for different geometry of the devices : ℓ Si = 4, 5 or 7 nm. These characteristics are comparable to the one obtained in [START_REF] Baccarani | A compact double-gate MOSFET model comprising quantum-mechanical and nonstatic effects[END_REF][START_REF] Negulescu | Simulation schemes in 2D nanoscale MOSFETs: a WKB based method[END_REF][START_REF] Pietra | Modeling and simulation of the diffusive transport in a nanoscale Double-Gate MOSFET[END_REF]. We present in Figure 4 the evolution of the temperature in the device with respect to the drain-source voltage for two different values of the Gate voltage. Figure 5 displays the evolution of the mean velocity, defined by J 1 /(qρ µ,T ), where the one dimensional density is given in (2.35) and the current in (3.17). As expected, we notice an overshoot of the velocity at the frontier between the channel and the drain for high value of V DS . We plot in Figures 6 the shape of the density in the device for two differents drain-source voltage. For V DS = 0V , we are at equilibrium and the density is symmetric in the device. For V DS = 0.5V , we notice transport of the charge carriers in the device.

Conclusion

A coupled quantum-classical model has been obtained for describing the transport of a partially confined electron gas. In a subband decomposition framework, the transport model is obtained by means of diffusive approximation from adiabatic quantum-kinetic models. The final system in the transport direction is obtained through two steps. First, under the assumption of dominant elastic scattering, a SHE system is derived (referring to [START_REF] Ben Abdallah | Adiabatic quantum-fluid transport models[END_REF]). Then, under dominant electron-electron collisions, an energy transport model is given, obtaining diffusion coefficients well suited for numerical purposes and with a relaxation term taking into account the electron-phonon interactions. In particular, with energy dependent crosssection of the collision operator explicit expression of the diffusion coefficients and of the relaxation term is derived and used for the numerical simulation of transport in a Double-Gate MOSFET. In the limit of large electrons-phonons collisions, we recover the model of [START_REF] Pietra | Modeling and simulation of the diffusive transport in a nanoscale Double-Gate MOSFET[END_REF]. We point out that a saturation of the current is observed without need of resorting at mobility modeling as done in [START_REF] Pietra | Modeling and simulation of the diffusive transport in a nanoscale Double-Gate MOSFET[END_REF].

Figure 1 :

 1 Figure 1: Schematic representation of the modeled device.

  )-(3.18) and (3.21)-(3.22) can be written in compact form as

  and µ h ∈ Z 2 h,0 . Equation (3.25) is derived from the weak formulation of the first equations in (3.23); (3.26) comes from the weak form of the second equations in (3.23); and finally, (3.27) imposes the continuity of the currents at the nodes.

Figure 2 :

 2 Figure 2: Left : I -V DS characteristics for V G = 0V and V G = 0.2V . The dashed line corresponds to φ ph = 10 -4 /φ 0 , the solid line corresponds to φ ph = 10 5 /φ 0 which is a good approximation to the drift-diffusion model. Right : Temperature in the device for φ ph = 10 5 /φ 0 . We see that the temperature is almost constant.

Figure 3 :

 3 Figure 3: I -V DS characteristics for different Gate voltages V G with ℓ Si = 5nm (left) and for different width of the Silicon in the DG-MOSFET with V G = 0 V (right).

Figure 4 :

 4 Figure 4: Evolution of the temperature in the device for a Gate voltage V G = 0V (left) and V G = 0.2V (right).

Figure 5 :

 5 Figure 5: Mean velocity for different drain-source potentials V DS and for V G = 0V (left) and V G = 0.2V (right).

Figure 6 :

 6 Figure 6: Density of electrons in the device for V DS = 0V (left) and V DS = 0.5V (right); in this simulation we take V G = 0.2V .

Table 1 :

 1 Table of the main physical values

	Parameter	Value	Length	Value
	N +	10 26

Acknowledgements. The authors acknowledge partial support from the Galilée Project no 25992ND of the Hubert Curien program : "Modèles numériques du transport collisionnel dans des dispositifs nano-électroniques".