A quantum Drift-Diffusion model and its use into a hybrid strategy for strongly confined nanostructures

In this paper we derive by an entropy minimization technique a local Quantum Drift-Diffusion (QDD) model that allows to describe with accuracy the transport of electrons in confined nanostructures. The starting point is an effective mass model, obtained by considering the crystal lattice as periodic only in the one dimensional longitudinal direction and keeping an atomistic description of the entire two dimensional cross-section. It consists of a sequence of one dimensional device dependent Schrodinger equations, one for each energy band, in which quantities retaining the effects of the confinement and of the transversal crystal structure are inserted. These quantities are incorporated into the definition of the entropy and consequently the QDD model that we obtain has a peculiar quantum correction that includes the contributions of the different energy bands. Next, in order to simulate the electron transport in a gate-all-around Carbon Nanotube Field Effect Transistor, we propose a spatial hybrid strategy coupling the QDD model in the Source/Drain regions and the Schrodinger equations in the channel. Self-consistent computations are performed coupling the hybrid transport equations with the resolution of a Poisson equation in the whole three dimensional domain.


Introduction
The extreme miniaturization reached in nanoelectronics brings the necessity of developing new models to describe the electron transport. Indeed, to a reduced channel length it corresponds also a strong reduced lateral dimension. When the cross-section diameter is below 3 nm, the strong confinement affects the energy band structure and bulk material quantities cannot be used in the simulations (see [15] e.g., and references therein). In particular, the potential generated by the crystal structure can be considered periodic only in the longitudinal direction and the assumption of infinite periodic structure, which allows to derive the commonly used effective mass approximation, is not valid anymore. Atomistic ab-initio computations give an accurate description of the electron transport but they are computationally too demanding and are not feasible in a device design framework. A challenge is to develop computationally efficient models that describe with accuracy the most relevant features of such confined devices.
This issue has been the subject of recent works. In [9], a quantum model has been proposed to describe the ballistic transport of electrons in ultra-scaled confined nanostructures. More precisely, an envelope function decomposition has been used to derive from a three dimensional (3D) Schrödinger equation a one dimensional (1D) longitudinal effective mass model, where device dependent effective quantities retain the effects of the confinement and of the transversal crystal structure. In the non degenerate case, adiabatic decoupling occurs and the model consists of a sequence of 1D effective Schrödinger equations, one for each band.
Interactions of charged particles with phonons have also been considered to describe these confined nanostructures. A drift-diffusion (DD) model has been derived and analyzed in [20]. It consists of a single macroscopic equation in which atomistic quantities, similarly to [9], are integrated. However, in promising miniaturized nanoelectronic components, the transport model must describe collisions as well as quantum effects. Since including collision terms in a quantum model is a quite complicate matter, a reasonable possibility consists in using a quantum macroscopic model.
In this work, we propose a formal derivation of a Quantum Drift-Diffusion (QDD) model in this context of strongly confined nanostructures. For that, we follow the theory developed in [13,12], based on an entropy minimization technique. It relies on an extension of Levermore's moment approach [25] to quantum systems. Notice also that, in the context of semiconductor modeling, the principles of extended thermodynamics were first used to derived macroscopic models of hydrodynamic nature (see [2] e.g.).
More precisely, we first identify a ballistic transport model. Starting from the Schrödinger system proposed in [9], we associate to each wave function a "density-matrix" function and we obtain a sequence of Wigner equations. We emphasize that, similarly to [9,20], it is a one dimensional model with device dependent effective quantities that retain the effects of the 2D transversal structure. Then, we heuristically include a description of collisions, adding a simple collisional term of Bhatnagar-Gross-Krook (BGK) type. The local equilibrium to which the system is relaxing is chosen, following the ideas of [13,12], as the extremum of an entropy functional subject to the constraint that the zeroth moment is prescribed. Next, we apply the Chapman-Enskog method and we obtain a non local QDD model. Finally, a semiclassical expansion of the quantum Maxwellian up to second order leads to a local QDD model. Along this paper, the term QDD is generally used to name the local model. Compared to the DD model derived in [20], the QDD model obtained in this paper contains a nonlinear term with a quantum correction to the potential that can be interpreted as a quantum potential similar to the so-called Bohm potential well known in quantum mechanics.
Generalized variants of QDD models have recently been proposed. For instance, the diffusive limit of a two band k.p model has been investigated in [4]. A multiband model with an arbitrary number of bands has been considered in [30] to describe the quantum transport in a strong force regime. In [5], a model has also been derived to describe a spin-polarized bidimensional electron gas. Two species of particles (labeled spin-up and spindown) are taken into account and a vectorial model is obtained. Our paper gets into this rich multiband framework. The novelty is the way to integrate into the QDD type model the atomistic information of the strongly confined transversal section. As consequences, the quantum correction that we obtain includes the contributions of the different energy bands and contains additional terms compared to the usual Bohm potential. Along this paper, we make the restrictive assumption that adiabatic decoupling occurs. In [9], a system of coupled Schrödinger equations is obtained in the degenerate case with a coupling between the energy bands through the potential. The derivation of a diffusive model in the degenerate case is far from the scope of this paper and is not discussed here.
In order to assess the capability of this QDD model to describe the transport of electrons in confined nanostructures, numerical experiments are performed for a zig-zag single-walled Carbon Nanotube Field-Effect Transistor (CNTFET) with a gate-all-around. The peculiar electronic properties of CNT, that strongly depend on the geometry of the tube, make them promising components in future FETs (see [29,16,26] e.g.). Self-consistent computations are performed coupling the one dimensional quantum transport model with the resolution of a Poisson equation in the whole three dimensional domain. Numerically, we observe that the quantum correction improves significatively the results compared to those obtained with the DD model derived in [20].
In order to complete the picture of possible models for CNTFET's, we also investigate in a second part the use of this QDD model in a hybrid strategy, spatially coupling it with the effective mass Schrödinger model proposed in [8]. More precisely, the macroscopic collisional QDD model is used in the Source/Drain regions and the Schrödinger equations are used in the active region. We emphasize that the computational cost of a hybrid strategy is generally cheaper than the one of a full Schrödinger model. In the literature, in devices where different effects are predominant in different part of the domain, spatial hybrid strategies have been designed, prescribing appropriate transmission conditions at interfaces. A hybrid kinetic-quantum model was first considered by N. Ben Abdallah in [7] where a Boltzmann equation is used to define the density in the classical zones instead a Schrödinger equation is chosen to describe the density in the quantum domain. At interfaces, reflection-transmission coefficients are defined to give the boundary conditions of the Boltzmann equation. Inversely, the distribution function is used as a statistical function to construct the quantum density. In [7], the author proves that the reflection-transmission conditions preserve the current.
Next, the diffusive limit has been considered in [11] where the Boltzmann equation and the reflection-transmission conditions have been replaced by a DD equation and appropriate interface conditions derived using a boundary layer analysis. In [6], the strategy to couple the drift-diffusion Schrödinger system is quite different since the coupling is direct and authors get an analytic expression of the connection conditions by writing the exact continuity of the current at interfaces. In [14], a QDD equation is coupled with a Schrödinger model. Due to the fourth order term appearing in the QDD equation, additional conditions are needed and the coupling is realized by assuming the continuity of both the electron density and the current at interfaces. In [11,6,14], numerical results are performed for a one-dimensional resonant tunneling diode.
Here, our hybrid strategy follows the one of [14]: we impose the continuity of the electron density and the current at interfaces. Again numerical results are performed for a CNTFET. We compare, analyzing the different current-voltage characteristics, this hybrid approach with the two non hybrid models (effective mass QDD or Schrödinger used in the entire domain) and also with the hybrid DD-Schrödinger strategy described in [19].
The paper is organized as follows. In Section 2, we present the effective mass Schrödinger model proposed in [9] to describe the electron transport in strongly confined nanostructures. We assume that adiabatic decoupling occurs and, associating to each wave function a "densitymatrix" function, we obtain a sequence of Wigner equations. This model is the starting point for the formal derivation of the macroscopic QDD model presented in Section 3. After a description of the obtained QDD model (Proposition 3.2), we detail the three main steps of the derivation: the entropy minimization (Section 3.1), the diffusive limit using a Chapman-Enskog method (Section 3.2) and the semiclassical expansion up to the second order (Section 3.3). Section 4 is dedicated to the description of the hybrid strategy. In particular, we explicit the interface conditions and present the equation discretizations. Finally, in Section 5, numerical simulations are performed for a gate-all-around CNTFET. We describe the physical device and compare the numerical results obtained with the different models.

Presentation of the multiband Wigner-BGK model
In [9], a novel quantum effective mass model has been derived by performing an asymptotic process which consists in using an envelope function decomposition to obtain a new effective mass approximation (see also [3] for a similar approach for 3D periodic crystals). We recall it briefly here. Let us consider an infinite wire defined in a physical domain R×ω , where is the typical spacing between lattice sites. As starting point, the transport is described by a scaled Schrödinger equation in R × ω containing a potential W generated by the crystal lattice, fast oscillating in the scale defined by the crystal spacing, and a slowly varying external potential V : where ψ in, is a given initial value. Here, is the reduced Planck constant and m e is the electron mass. Since the 2D cross-section ω comprises few ions, W is considered periodic only in the longitudinal x-direction (transport direction) and the variable z of the transverse section can be considered as fast variable, rescaled as z = z . Denoting by ω the scaled cross-section, we consider the following Bloch-type problem (with a quasi momentum equal to 0) in the 3D cell ∆χ n + W χ n = E n χ n , χ n (y, z ) = 0 on ∂ω, χ n 1-periodic in y, U |χ n | 2 dydz = 1. (2) Here y denotes the transport variable in the cell. The peculiarity of the strongly confined structure is reflected in the choice of the unit cell problem of Bloch type (2). We point out that this unit cell U comprises the entire cross-section of the nanostructure. Thus, the eigenvectors depend on the device under consideration, for instance on the device geometry, on the number of atoms, on the chirality, and so on. Moreover, the homogeneous Dirichlet condition imposes confinement in the transverse directions, while periodicity is considered only in the transport direction. Consequently, the eigenvectors are 3D quantities but the Brillouin zone and the associated energy bands are one dimensional.
The physically relevant potentials are nonnegative potentials given in L ∞ (U). Consequently, the eigenfunctions χ n , solutions of (2), form an orthonormal basis of L 2 (U), with real eigenvalues which satisfy Moreover, from the min-max principle, it is clear that we have E n ≥ Λ n for all n ∈ N, where Λ n are the eigenvalues of the Laplacian problem in the unit cell U with Dirichlet boundary conditions imposed on ∂ω and periodic conditions in the longitudinal direction. It is wellknown that for all λ > 0, n∈N e −λΛn < +∞ (eigenvalue properties for a Laplacian-Dirichlet problem [18]). Thus, we deduce that The asymptotic process developed in [9], by using the functions defined in (2) as basis for the envelope function decomposition, allows to average out not only the lattice potential, but also the lateral dimension. The limit problem consists of an infinite set of Schrödinger equations, whose structure depends on the multiplicity of the eigenvalues (3). Assumption 2.1. Along this paper, we assume that the eigenvalues E n are all simple.
Then, in the limit system adiabatic decoupling occurs and, for each energy band, the Schrödinger equation, which incorporates the relevant averaged quantities, based on the Bloch functions, has the form The n-th band effective mass m * n is defined by Also, the effective potential can be computed and it is given by Remark 2.2. In the degenerate case, to each multiple eigenvalue corresponds, instead of equation (5), a system of coupled Schrödinger equations with dimension equal to the multiplicity of the eigenvalue. The kinetic part of the effective mass Hamiltonian is diagonal and the coupling occurs through the potential.
The first step towards the definition of a macroscopic quantum model is the construction from (5) of the corresponding Wigner equation. For each band n, let us introduce the "densitymatrix" (function) ρ n associated to the wave function ψ n solution of (5). It is defined by ρ n (t, r, s) = ψ n (t, r)ψ n (t, s) (8) and it is solution of We also need the so-called Wigner transform of an operator ρ acting on functions of L 2 (R). It is a function in the phase space (x, p) ∈ R 2 first used by Wigner in 1932 [31] and defined by where ρ(r, s) is the integral kernel of the operator ρ : ρϕ(r) = R ρ(r, s)ϕ(s)ds, for any function ϕ ∈ L 2 (R). Inversely, the inverse Wigner Transform, also called Weyl quantization, transforms a function w(x, p) into an operator ρ acting on functions of L 2 (R) in the following way In our case, the Wigner function w n = W (ρ n ) constructed from the density matrix (8) is formally a solution of the following Wigner equation where v n is a velocity defined by v n = p m * n and Θ [V n ] is the operator given by We also recall that 1 2π R w n (t, x, p)dp represents the electron density in the n-th band.
We now need to include a collision mechanism into the ballistic model just described. Since the main information about the microscopic collision mechanism entering the macroscopic model that we will derive is the form of the local equilibrium, we will use along this paper a simple collision operator Q 0 in the BGK form. It writes where τ n is a relaxation time related to the n-th band and f (w n ) is a function representing an equilibrium state that will be specified in the next section. can be seen as a diffusive constant D. It will naturally appear in the derivation of the macroscopic model in the next section. As we will see, the obtained macroscopic model consists of a single equation and not of a sequence of equations, one for each band. Consequently, it is also interesting to rewrite this diffusive constant D as where τ is the relaxation time of the device and m * G a global effective mass given by Finally, we write the Wigner equation (12) in the diffusive scaled form, introducing a small parameter α proportional to the mean free path. Consequently, we formally obtain, for each energy band, the following Wigner-BGK scaled equation In the sequel, we will use this set of decoupled Wigner-BGK equation as starting point to derive a quantum macroscopic model in a multiband setting.

Derivation of the Quantum Drift-Diffusion model
In this section, we perform the formal derivation of a Quantum Drift-Diffusion type model, following an approach proposed by Degond et al. [13,12]. It is based on entropy minimization techniques and a definition of the so-called "quantum logarithm" LOG and "quantum exponential" EXP . They are phase space functions defined, using definitions (10) and (11), as We underline that the quantum logarithm and the quantum exponential involve heavily nonlocal operations. It is essential to recover the nonlocality of quantum mechanics. This approach, that is detailed afterwards in Sections 3.1 and 3.2, allows us to obtain the following result: Proposition 3.1. Let V be a given potential and w α = (w α n ) n∈N be a sequence of Wigner functions w α n solutions of (17) with initial data w 0 = (w n,0 ) n∈N . Then, as α goes to 0, w α n formally converges to with A (appearing through (N, J)) solution of the non local conservation equation where the electron density N , the current J and the initial density N 0 are respectively defined by and We emphasize that this macroscopic model is not made of a sequence of equations, one for each energy band. On the contrary, it consists in a single conservation equation similarly to the Drift-Diffusion model analyzed in [20]. As we will see in Subsection 3.1, the single macroscopic equation is the result of the use of an entropy in which, thanks to the adiabatic decoupling assumption, the contributions of each band are summed up.
The extremal solution of the entropy minimization problem will be given by the quantum exponential EXP A − E n − V n − p 2 2m * n that appears in (19). In the sequel, this extremal solution will be called quantum Maxwellian and denoted M n (w n ). In other words, Proposition 3.1 tells that w α n converges to the quantum Maxwellian as α goes to 0. The model described in Proposition 3.1 is nonlocal due to the fact that the quantum exponential involves non-local operators. However, a local model can be obtained thanks to a O( 4 ) expansion of the quantum Maxwellian. Following the computations that will be detailed in Section 3.3, we obtain the so-called Quantum Drift-Diffusion (QDD) model: Let A be the solution of the non local equation (20) with J and N defined by (22) and (23), respectively. Then, we formally have where N q satisfies the following conservation equation up to order 2 with V G is a global effective potential defined by and Q[N q ] is a quantum correction expressed as where the notation U denotes the operator We remind that D is the diffusive constant defined in (15), m * G the global effective mass (16) and n i the intrinsic carrier concentration that is expressed (see [1] e.g.) as All these physical quantities are device dependent since they depend on the effective masses m * n and the energy bands E n . So, this novel QDD model is peculiar to the device under consideration.
We also emphasize that the term −  (27) (through (29)) is similar to the so-called Bohm potential. Here, the quantum correction includes a sum over n to take into account the contributions of the different bands. Also, a notable difference is that Q[N q ] contains additional terms depending on V G − V n . They retain the difference between an effective potential V n (7) defined for each band and the global effective potential V G (28).
Remark 3.3. The current J q (27) can be rewritten in terms of a quantum quasi-Fermi energy A q : It is this form that will be used in the numerical part.
Remark 3.4. We can consider without restrictions that V is a nonnegative potential. Using (4), we can immediately say that n∈N e −En−Vn < +∞. So, all the quantities such as (28) are well defined.
Let us now detail in the following subsections the formal derivation that allows us to establish Propositions 3.1 and 3.2.

Entropy minimization technique
The starting point of the approach consists in defining an equilibrium state as the minimizer of a suitable quantum entropy functional under the constraint on the zero order moment. Since the moments are naturally defined in terms of Wigner functions, we have chosen here to express the entropy in the Wigner formalism (by invoking the Wigner and the inverse Wigner transforms) even if it is more natural to define it in terms of the density operator. So, the entropy is expressed as where s is the convex function defined by s(f ) = f (log(f ) − 1). Following [13] (Lemma 3.3) and using the definition of the quantum logarithm LOG (18), it can be proved that, for all functions g, Since we are interested in this paper in the derivation of an isothermal model, we have to take into account the total particle energy minimizing the relative entropy For any w = (w n ) n∈N in E = {w ∈ l 1 (L 1 (R 2 )) s.t. E(w) < +∞}, we associate a density N [w](x) = 1 2π n∈N R w n dp. Then, the minimization problem that we consider is the following: for a given electron density N > 0, we look for A rigorous proof of existence and uniqueness of a similar constrained minimization problem has been obtained by Méhats and Pinaud [28,27]. Introducing a Lagrange multiplier A that is independent of the energy bands, we consider the functional in which the constraint is weakly imposed. It can formally be proved that for each n, the extremal condition for w = (w n ) n∈N imposed for all functions g implies that The extremal solution w n is called the quantum Maxwellian and will be denoted M n (w n ) in the sequel. Moreover, A can be interpreted as a quantum quasi-Fermi energy.

Formal diffusive limit
Now, we can perform the diffusive limit when α goes to 0. The starting point is the Wigner-BGK equation (17) that, to simplify the notations, can be written introducing the operator Γ = v n ∂ x − Θ [V n ]. In order to specify the equilibrium state of the BGK collision operator defined in (14), we use the quantum Maxwellian defined in (37): The collision operator has the following properties (see [21] e.g.): We also remind (see e.g. Lemma 12.9 in [21]) that the operator Θ [V n ] defined in (13) is such that R Θ [V n ]w α n dp = 0 (41) and R pΘ [V n ]w α n dp = −∂ x V n R w α n dp.
Passing to the limit α → 0 in (38), we obtain that Q 0 (w n ) = 0, where w n = lim α→0 w α n . It means that w n = M n (w n ). Secondly, we insert a Chapman-Enskog expansion w α n = M α n (w α n )+ αg α n in (38). It gives Passing to the limit α → 0, we obtain τ n ΓM n (w n ) = M n (g n ) − g n , where g n = lim α→0 g α n . Thirdly, we take (38) and we integrate it over p. Using (41), we obtain α∂ t R w α n dp + ∂ x R v n w α n dp = 0.
We insert the Chapman-Enskog expansion in this expression and we use the fact that pM α n dp = 0 since it is an odd function. It gives ∂ t R M α n (w α n ) + αg α n dp + ∂ x R v n g α n dp = 0.
Passing to the limit α → 0 and summing on each band n, we obtain after calculations, thanks to the property (42), the non local conservation equation (20) with (22)-(23).

Expansion in power of 2
Furthermore, in order to obtain formally the local model stated in Proposition 3.2, we expand M n (w n ) in terms of 2 . Up to order O( 4 ), the following expansion (see Proposition 5.3 in [12]) holds for the quantum Maxwellian (for all x ∈ R and p ∈ R): Notice that the zeroth order term corresponds to the classical Maxwellian. Integrating over p, it leads to where U has been defined in (30). Using the expression of the intrinsic carrier concentration (31), we obtain, for the density N defined in (22), the following expansion where Z n has been defined in (28). Now, we would like to express A in terms of N . Using the global potential V G defined in (28), we can write Moreover, we have Differentiating (45) with respect to x, we find ∂ x (A − V G ) = ∂xN N + O( 2 ) and consequently Therefore, we have Putting l Z l in factor in the second logarithm, we finally obtain where Q[N ] is the quantum correction defined as in (29).
The last step consists in expanding the current defined in (23). Since Under Assumption 2.3 and using the definition of the diffusive constant (15), we obtain thanks to (47): Introducing the expansions obtained above into the conservation equation (20), the result formulated in Proposition 3.2 follows.

Self-consistent computations
Finally, in preparation for self-consistent computations, given N q solution of (25), we need to define a charge density N n q for each band. We do so, recalling that N = 1 2π n∈N R M n (w n )dp. Consequently, we define N n = 1 2π R M n (w n )dp.

From (44), we can write
Thus, a reasonable choice for N n q is to omit the remainder of order O( 4 ). Rewriting the dominant term in terms of N q , it gives where S n has been defined in (46). We obviously notice that the sum over all the bands gives the density N q . Next, as in [9] and [19], the transformation from the one dimensional transport direction to the entire nanostructure is done by means of the quantities g n 's (7). It leads to the definition of a three dimensional density Consequently, the 1D transport model presented in Proposition 3.2 can be self-consistently coupled with the following 3D Poisson equation for the electrostatic potential V P q is the elementary charge, 0 the permittivity in vacuum, r the relative permittivity and N dop the prescribed doping density. Then, V in (7) is given by V = −qV P .

Description of the hybrid strategy
In the following, we use the (stationary) QDD model derived in the previous section into a hybrid approach. Indeed, we propose to couple it, spatially in the transport direction, with the Schrödinger model derived in [9]. So, we consider a bounded domain Ω = (x L , x R ) × ω, where (x L , x R ) is the transport domain and ω is the two dimensional strongly confined cross-section. We assume that the device domain in the transport direction x is partitioned into a Schrödinger zone S = (x I 1 , x I 2 ), with x L < x I 1 < x I 2 < x R (where the Schrödinger system is used) and a quantum macroscopic zone Q = (x L , x R )\S (where the transport is described by the QDD model). For a Field-Effect Transistor (FET), which is the kind of device that we will consider in Section 5, the Schrödinger zone S corresponds to an active zone where quantum ballistic effects are predominant whereas the quantum macroscopic zone Q is made of two electron reservoirs called Source and Drain in which the transport is considered highly collisional.
First, we recall the Schrödinger model used in the bounded domain S. Afterwards, we describe the interface conditions which preserve the continuity of the total current and the electron density. We present the algorithm used to couple the two models and the discretization of equations for a given potential. Finally, we consider self-consistent computations and we detail the coupling between the transport equations and the Poisson equation.

Description of the Schrödinger system
The quantum transport in S is given by the sequence of Schrödinger equations (5) in the stationary framework. The system is considered as an open system in (x I 1 , x I 2 ), with Transparent Boundary Conditions (TBCs) supplementing the equations (see [8,24] e.g.). To fix ideas, we detail the case V n (x I 1 ) ≥ V n (x I 2 ). For each band and for each wave vector k, we consider the following stationary Schrödinger equation where Defining the coefficients the TBCs are written for k > 0 as and for k < 0, we have The reflection and transmission amplitudes r n (k) and t n (k) of the wave functions are determined by and t n (k) = ψ k n (x I 1 ) for k < 0.
Finally, we define the reflection coefficients as R n (k) = |r n (k)| 2 and the transmission coefficients T n (k), corresponding to the proportion of incident electrons which are transmitted, as Electrons are considered in mixed states and the one dimensional electron density in the S region is defined by where the density carried by the n-th band is given by superimposing the densities of states injected from the reservoirs, that is φ n (k) is a given statistical function which characterizes the electron injection from reservoirs. Finally, the current is defined by where the n-th band current is given by It can be easily seen that the current does not depend on x. Furthermore, using the TBCs of the Schrödinger equation (54) and (55), as well as the properties of the transmission coefficients (56), J n s can be expressed in the following form The choice of the function φ n is important in order to perform the coupling with the QDD model. Our choice is to use the classical Maxwellian expressed in terms of the quantum quasi-Fermi energy variable A q (32) (computed with the QDD equation) at the interfaces where f is the Boltzmann statistics defined by We emphasize that this function is the zeroth order term of the quantum Maxwellian (19).

Implementation of the hybrid model
We assume at this point that the electrostatic potential V is given (and consequently also the effective potentials V n (7) and V G (28) are given). The spatial coupling between the Schrödinger system and the QDD model is realized through interface conditions. Since the QDD equation contains fourth-order derivatives, two conditions are required at each interface. Following [14], we assume that, at the interface points x I 1 and x I 2 , the electron density and the current are continuous. It gives the following conditions For the transport direction x, we introduce a uniform mesh with N points (that includes the interface points x I 1 and x I 2 ). It is defined by the abscissa for l = 0, ..., N −1. We also denote N Q the number of mesh points for each macroscopic region [x L , x I 1 ] and [x I 2 , x R ] and N S = N − 2 N Q + 2 the number of mesh points for the Schrödinger zone [x I 1 , x I 2 ]. Thus, the quantum macroscopic region Q is described by the points ξ l defined by ξ l = x L + l∆x for l = 0, ..., N Q − 1 and ξ l = x I 2 + (l − N Q )∆x for l = N Q , ..., 2 N Q − 1.
We detail here the different steps of the hybrid coupling: • First, we solve on S the Schrödinger equations (52) with TBCs (54)-(55) for each wave vector k and each band n. Each equation is transformed to an initial value problem and discretized with a Crank-Nicolson scheme (see [9] for details).
• We compute the transmission coefficients T n (56) as well as the following quantities and The integrals are computed with a trapezoidal quadrature rule, using a constant momentum step ∆k enough refined to take into account the contribution to each significant energy.
• Next, the stationary QDD model (see Proposition 3.2 and Remark 3.3) is approximated by a central finite differences as in [14]. This scheme has been proved to be positivity preserving (see [22] for details). For that, we introduce the functions = N q and V n = V G − V n . We also denote by A l , l and V n,l the quantities A q (ξ l ), (ξ l ) and V n (ξ l ). In order to make easier the implementation of the interface conditions, we keep two unknowns: the square root of the electron density and the quantum quasi-Fermi energy A q . So, for all l ∈ {1, ..., N Q − 2} ∪ { N Q + 1, ..., 2 N Q − 2}, we obtain the discretized equations and and C n,l = V n,l+1 −2 V n,l + V n,l−1 2 n,l 4 . Since we assume that no quantum effect occurs at boundaries, the discretized system is completed with the following boundary conditions Finally, we impose the interface conditions (64)-(65). Using the Boltzmann statistics, the term containing the quantum quasi-Fermi energy variables in (62) enters as a multiplication factor. Consequently, the quantum density (58) and the quantum current (61) can be written in a simpler form. Thanks to the computed quantities (66) and (67), the interface conditions are thus expressed explicitly: The full nonlinear system (68)-(75) is solved by means of a Newton algorithm. All the entries of the Jacobian matrix can be evaluated explicitly. The resulting linear system is solved by a GMRES solver with incomplete LU factorization for preconditioning.
• Finally, thanks to A q (x I 1 ) and A q (x I 2 ) that we just computed, we can complete the quantum density N n s (58) using (62). We are now able to define, for each energy band n, the charge density of the hybrid model In preparation for self-consistent computations, we can also compute a 3D density N 3D similarly to (50).

Coupling with the Poisson equation
We now consider that the macroscopic potential V is solution of the 3D Poisson equation (51). Because of the highly nonlinear coupling between the density and the potential equations, we use an iterative method of Gummel type [17] that we quickly remind here. Let V old P be a given electrostatic potential. We define the potential energy V = −qV old P and we compute the 3D density N 3D coming from the hybrid Schrödinger Quantum Drift-Diffusion approach described in the previous Subsection. Then, we solve the 3D Poisson equation using piecewise linear finite elements on a prismatic mesh (see [19] for details), modified according to the Gummel iteration algorithm, as in [10], that is We repeat it until the quantity V old P − V new P L ∞ becomes sufficiently small. For the device presented in the next Section, a potential is applied between the two electron reservoirs (Source and Drain). We first consider the whole iterative system at thermal equilibrium (zero applied Drain-Source voltage V DS ). Once the convergence is reached, we increment V DS and start a new iteration. For the simulations presented in this paper, a reasonable increment step is 0.02 V.

Presentation of numerical results
In this paper, we compare the numerical results obtained with five different approaches: -approach S: the Schrödinger model proposed in [9] is used in the entire domain (x L , x R ), -approach DD: the Drift-Diffusion model proposed in [20] is used in the entire domain (x L , x R ), -approach QDD: the Quantum Drift-Diffusion proposed in Section 3 is used in the entire domain (x L , x R ), -approach S-DD: the hybrid approach detailed in [19] couples spatially the Schrödinger system with the Drift-Diffusion model, -approach S-QDD: the hybrid approach described in Section 4 couples spatially the Schrödinger system with the Quantum Drift-Diffusion model.

A Gate-all-around Carbon Nanotube Field-Effect Transistor
Numerical simulations are carried out for a Carbon Nanotube Field-Effect Transistor (CNT-FET). A section along the transport direction (x-axis) is presented in Fig.1. It contains a (10,0) zig-zag single-walled CNT (represented in green and yellow in Fig.1) surrounded by a layer of dielectric SiO 2 ( r,ox = 3.9) of 1.4 nm thickness acting as an insulator (in red in Fig.1). In Fig.2, a representation of carbon atom positions is presented for the (10,0) zig-zag CNT. They are placed on a circle of 0.78 nm diameter. The relative permittivity of a CNT is a device dependent quantity, that here has been taken as r,c = 13. In the last part of this section, we will also consider the case of an anisotropic permittivity.  . Finally, a Gate is imposed all-around the transversal structure to modulate the number of free electrons. Since the effect of changing the gate voltage is as expected in the experiments that we did, we present here only results corresponding to a gate voltage V GS equal to V GS = −0.1 V and we refer to [19] for a discussion on the gate voltage impact.
Along this paper, the transport-Poisson problems are solved for the 9 first conduction bands. As described precisely in [19], this choice is widely sufficient since for our (10,0) zig-zag CNT, the energy levels increase quickly and consequently only the first bands give a significant contribution to the total current. These energy bands are not all non degenerated as required by Assumption 2.1. However, it numerically turns out that, for the device under consideration with the gate-all-around, the off-diagonal terms in the potential matrix are negligible and a decoupled system can still be considered.
Finally, we would like to make some comments about the electron mobility constant µ that is related to the diffusive constant D (expressed in scaled form in (15)) by the Einstein relation. For strongly confined structures, it is a device dependent physical parameter and a well established value is not found in the literature. Therefore, we first discuss how this value affects the computed current in our simulations. Solving the transport problem with the QDD model provides a current which is proportional to the mobility constant. It is not the case with the S-QDD model as we can see in Fig.3 where the output current-voltage characteristics obtained with the S-QDD model are presented for different mobilities. We notice that the two typical regimes (an ohmic regime for small values of V DS and then a quasi-saturation regime) are always observed. But, above all, we observe that for mobilities larger than 5 × 10 −3 m 2 .V −1 .s −1 the current is virtually not modified.
Since our goal is to discuss the effects of the quantum correction term in the QDD and S-QDD models rather than to give a quantitative description of IV curves, we fix in the sequel the electron mobility value. As done in [19], we have chosen µ = 0.5 × 10 −4 m 2 .V −1 .s −1 . The output current-voltage characteristics obtained with this choice are presented in Fig.4 for the five different approaches. Due to this small mobility constant (that corresponds to a highly collisional regime), the DD and the QDD currents are much smaller than the Schrödinger current. The two hybrid approaches (S-DD and S-QDD) that take into account both electronphonon collisions in the reservoirs and ballistic quantum effects in the active zone, give an intermediate current that is however closer to the quantum one.

Effects of the quantum correction term
First, we start to present 2D slices along the transport direction (x-axis in the pictures), cutting the cross-section in the middle. Fig.5 represents the self-consistent potential at thermal equilibrium and Figs.6-7 represent the density in logarithm scale respectively at equilibrium and for an applied voltage V DS = 0.2 V. The left pictures are obtained with the DD model in-   stead the right ones with the QDD approach. Clearly, confined electron channels are apparent between Source and Drain in the carbon regions whereas the oxide layer acts as insulator. We observe that the quantum correction term added in the QDD approach has a non negligible effect close to the doping discontinuities (x=10 nm and x=20 nm). In particular, the density is smoother with QDD than with DD.
In order to perform a better comparison, we also present 1D profiles. It corresponds to averaged quantities resulting from an integration of the 3D quantities over the 2D wire section divided by the wire section area. The potential (on the left) and the density (on the right) are presented at thermal equilibrium in Fig.8 and for an applied voltage V DS = 0.2 V in Fig.9. For a better visualization, it is also interesting to plot the inverse of the density (Fig.10). In all figures, three approaches are compared : the Schrödinger model in dashed lines, the Drift-Diffusion one in dotted lines and the Quantum-Drift-Diffusion one in solid lines. The results confirm that the transport in such confined structures is strongly governed by quantum effects. The DD model that is purely classical gives not enough accurate results.

Hybrid strategy
We are now interested in the spatial hybrid strategies. In Figs.11-12, again three approaches are compared : the Schrödinger model in dashed lines, the hybrid Schrödinger-Drift-Diffusion one in dotted lines and the hybrid Schrödinger-Quantum-Drift-Diffusion one in solid lines. The potential (on the left) and the inverse of the density (on the right) is presented at thermal equilibrium in Fig.11 and for an applied voltage V DS = 0.2 V in Fig.12. The clear difference observed between the S-DD and the S-QDD curves emphasizes that the quantum correction plays an important role not only in the active zone but also in the collisional reservoirs. Finally, we also would like to highlight that the computational cost of the S-QDD approach is much cheaper than the one of the full Schrödinger model. In particular, the large number of Schrödinger equations at each iteration step (one for each energy band and each wave vector) are performed on a smaller domain (divided by 3 for our device).

Interface positions
For the S-QDD model, we now study the influence of the interface positions. In previous figures, the interfaces were located at x I 1 = 10 nm and x I 2 = 20 nm (at doping discontinuities). In Figs.13-15, in the left pictures, x I 2 is fixed at 20 nm and we move x I 1 , and inversely, in the right pictures, x I 1 is fixed at 10 nm and we move x I 2 . The potential is presented in Fig.13, the inverse of density in Fig.14 and the Current-Voltage characteristics in Fig.15. On the one hand, we observe perceptible differences when one interface is placed inside the active zone (see Fig.14 for instance). On the other hand, in Fig.15, the saturation current stays almost unchanged when the interface x I 2 goes to Drain, instead the saturation current increases up to the quantum value when the interface x I 1 reaches Source. These results are qualitatively similar to those presented in [19] for the S-DD model and they confirm that modeling electron transport (collisional vs ballistic) is important, specially in the Source.

Anisotropic permittivity
To finish, we would like to say few words about the relative permittivity of the CNT. Up to here, we have taken r,c = 13. However, the polarizability in a CNT is a complex phenomenon [23], that depends on the geometry and that has different values in the longitudinal and in the transverse directions. In Figs.16-17, we present the potential, the density and the IV curves obtained with S-QDD for the isotropic case (solid lines) and for an anisotropic case that corresponds to a longitudinal permittivity r,c,l = 142 and a transversal one r,c,t = 10.9 (dashed lines). As we can see, the qualitative behavior of the results is preserved but a significant increase of the current value is noticed due to the larger permittivity in the longitudinal direction. A similar behavior is observed for full Schrödinger. The 3D electrostatic effects are important and have to be investigated more deeply in the future.

Conclusion
In this work, we derive, using an entropy minimization technique, a QDD model in the context of strongly confined nanostructures integrating the atomistic information of the transversal section through the definition of effective quantities. Numerical simulations performed for a gate-all-around CNTFET allow to verify that the quantum correction appearing in this QDD model improves significatively the results compared to those obtained with the analogous DD model [20].
In a second part, we use this QDD model into a hybrid strategy, spatially coupling it with the effective mass Schrödinger model proposed in [9] and imposing the continuity of the electron density and the current at interfaces. Numerical simulations show a clear difference between the S-DD and the S-QDD results, showing that the novel quantum correction plays an important role also in the collisional reservoirs.
IV curves obtained with a hybrid approach are less sensitive to the device dependent mobility constant than those corresponding to full macroscopic models (DD or QDD). Additionally, in comparison with the full quantum ballistic model [9], the computational cost of the hybrid strategy is much cheaper, since the large number of Schrödinger equations is performed on a smaller domain and a single macroscopic equation is used elsewhere.