Exploitation of the Maximum Entropy Principle in Mathematical Modeling of Charge Transport in Semiconductors

Abstract

In the last two decades, the Maximum Entropy Principle (MEP) has been successfully employed to construct macroscopic models able to describe the charge and heat transport in semiconductor devices. These models are obtained, starting from the Boltzmann transport equations, for the charge and the phonon distribution functions, by taking—as macroscopic variables—suitable moments of the distributions and exploiting MEP in order to close the evolution equations for the chosen moments. Important results have also been obtained for the description of charge transport in devices made both of elemental and compound semiconductors, in cases where charge confinement is present and the carrier flow is two- or one-dimensional.

1. Introduction

In the design of new generation electronic devices, the modern micro- and nanoelectronics industry has an increasing need for mathematical models to simulate devices before they are realized in the laboratory. In particular, the accurate modeling of energy transport in semiconductors is necessary in order to describe high-field phenomena, such as hot electron propagation, impact ionization and heat generation in the bulk material. Semiclassically, the most effective way to describe these phenomena makes use of the semiclassical Boltzmann Transport Equation (BTE) [1]. However, solving the BTE is a daunting computational task, both using an indirect stochastic approach by Monte Carlo methods and direct numerical schemes based on discontinuous Galerkin methods or finite differences [2]. This is the reason why, in many cases, it is desirable to have simpler macroscopical models, known as hydrodynamical-like models, which are highly useful for computer aided design (CAD) purposes. These models are obtained from the infinite set of moment equations of the BTE by a suitable truncation procedure. It is well-known that a closure assumption is required in order to have a closed system of evolution equations. In the past, various closure assumptions have been made for the semiconductor transport moment systems, leading to various classes of hydrodynamical models, see e.g., [3,4,5]. However, these various closure assumptions are, at best, only phenomenological and often a consistent physical and mathematical justification is lacking. Lately, a closure assumption based on the Maximum Entropy Principle of extended thermodynamics [6,7] has been successfully applied, both in the parabolic and non-parabolic band approximation, to various types of semiconductors [8,9,10,11,12,13]. The resulting models, which differ for the choice of the moments to assume as field variables, are, in fact, able to describe charge transport due both to electrons and holes and also heat transport due to phonons. All the main scattering mechanisms between carriers and phonons and among phonons themselves are taken into account. The models also have nice mathematical properties, being symmetric hyperbolic. In particular, this assures the well-posedness of the Cauchy problems and the finite velocity of the propagation of disturbances [14].

Due to the increasing shrinking of modern device dimensions, quantum effects are beginning to play a relevant role in charge transport. In the framework of the moment method and of Maximum Entropy Principle (MEP), a strategy to take into account these effects has been proposed in [15] and consists of using the moment system arising from the Wigner equation. The starting point is to expand the Wigner function and the relative transport equation with respect to the squared reduced Planck constant ${\hslash }^2$. The zero-order part of the collision operator is supposed to be the same as the semiclassical one, while the first-order contribution is supposed to act only on the ${\hslash }^2$ correction of the Wigner function and is modeled in a relaxation form. Therefore, at zero order, the Wigner function is given by the solution of the semiclassical Boltzmann equation, which is approximated with the standard maximum entropy method, while the ${\hslash }^2$ order correction is obtained with a Chapman–Enskog expansion in the high field scaling.

In the description of charge transport in some devices, such as double gate metal oxide semiconductor field effect transistors (DG-MOSFETs), where quantum effects are relevant only along one direction, called the confinement or transversal direction, another strategy can be used. One can adopt a quasi-static description along the confining direction based on a coupled Schrödinger–Poisson system which leads to a subband decomposition of the electron energy levels, while the transport along the longitudinal directions can be described by a semiclassical Boltzmann equation for each subband [16]. Therefore, a complete description can be done in terms of a coupled Schrödinger–Poisson–Boltzmann system. However, the numerical integration of the transport part, which has been performed by employing Monte Carlo or deterministic methods [16,17], is also in this case very expensive, from a computational point of view, for CAD purposes. Consequently, it can be convenient to substitute the Boltzmann equations with macroscopic models again obtained by using the moment method and closing the moment equations with MEP [18,19,20].

In this paper, we will give a review of all the above-mentioned models according to the following plan. In Section 2, the 3D semiclassical macroscopic models are presented, showing how they are closed by using MEP and in Section 3 their quantum correction is discussed. In Section 4, the case of quantum confinement is described. Eventually in Section 5, some numerical simulations are sketched.

2. The 3D Semiclassical Macroscopic Models

In the semiclassical description, roughly speaking, the charge and heat transport in semiconductors can be described by modeling them as consisting of a mixture of gases of electrons and phonons. As regards electrons, those which mainly contribute to the current flow occupy states near to the lowest minima of the lowest conductions bands and to the highest maxima of the highest valence bands. The latter contribution can be conveniently treated in terms of pseudo-particles called holes (missing electrons), which have electric charge, crystal momentum and energy opposite to those of the missing electrons in the considered valence band [1,21,22]. An electron/hole population for each neighborhood of the lowest/highest minima/maxima of the lowest/highest conduction/valence bands is taken into account. The neighborhoods of the minima/maxima are usually called valleys. The electron/hole state is characterized by the index of the band it occupies, by its wave vector and spin state.

The dependence of the carrier energy on the wave vector in each valley can be analytically approximated by non-parabolic dispersion relations of the following form [1,9]

$${\mathcal{E}}_A({\mathbf{k}}_A)=\frac{{\hslash }^2{\left|{\mathbf{k}}_A\right|}^2}{2m_e}{\gamma }_A({\mathcal{E}}_A){\psi }_A^{-1}({\mathbf{n}}_A),$$

(1)

where ${\mathcal{E}}_A$ is the carrier energy in the A-th valley measured from the bottom of the valley, the index A running over the considered valleys, ${\mathbf{k}}_A$ is the carrier quasi-wave vector, referred, for each valley, to the minimum or maximum of the valley [1], ${\mathbf{n}}_A:=\frac{{\mathbf{k}}_A}{|{\mathbf{k}}_A|}$, ${\gamma }_A$ is the non-parabolicity factor, and $m_e$ is the free electron mass.

For electrons, if one uses the ellipsoidal approximation, the dependence of ${\psi }_A$ on ${\mathbf{n}}_A$ is given by

$${\psi }_A^{-1}=\frac{{({({\mathbf{n}}_A)}_1)}^2}{{(m_A^{*})}_1}+\frac{{({({\mathbf{n}}_A)}_2)}^2}{{(m_A^{*})}_2}+\frac{{({({\mathbf{n}}_A)}_3)}^2}{{(m_A^{*})}_3},$$

where ${(m_A^{*})}_i^{-1},$ $i=1,2,3,$ are the diagonal elements (eigenvalues) of the inverse effective mass tensor of the A-th valley, multiplied by $m_e$, referred to an orthonormal basis of the tensor.

Analogously, for holes, one has

$${\psi }_A^{-1}=\left|A_A\right|\mp \sqrt{B_A^2+\frac{C_A^2}{4}\left({\sin }^4\vartheta {\sin }^{2}2\phi +{\sin }^{2}2\vartheta \right)},$$

in the case of warped bands, with φ and θ the azimuthal and polar angle of the wave vector respectively, and $A_A$, $B_A$ and $C_A$ inverse valence band parameters.

In the analytic approximations, for each valley, the domain of the wave-vector is extended to all ${\mathbb{R}}^3$ and the volume element in the $\mathbf{k}$-space can be written as (henceforth we omit the valley index unless there is a possibility of confusion)

$$d^3\mathbf{k}=\frac{m_e\sqrt{2m_e}}{{\hslash }^3}\sqrt{\frac{\mathcal{E}}{{\gamma }^5(\mathcal{E})}}\left(\gamma (\mathcal{E})-\mathcal{E}{\gamma }^{\prime }(\mathcal{E})\right){\psi }^{\frac{3}{2}}(\phi ,\theta )d\mathcal{E}d\mathsf{\Omega },$$

where the prime denotes a derivative with respect to the argument of the function, and $d\mathsf{\Omega }$ is the solid angle element. The charge carrier velocity, given by $\mathbf{v}=\frac{1}{\hslash }{\nabla }_{\mathbf{k}}\mathcal{E}$, has the following expression in terms of the energy and the angular variables

$$\begin{array}{ccc}\hfill v_i& =& \frac{\hslash }{2m_e}\frac{{\gamma }^2(\mathcal{E})}{\gamma (\mathcal{E})-\mathcal{E}{\gamma }^{\prime }(\mathcal{E})}{l}_i,l_i:=\frac{\partial }{\partial k_i}{(|\mathbf{k}|}^2{\psi }^{-1})=g(\mathcal{E}){\eta }_i(\phi ,\theta )n_i,\hfill \\ \hfill g:& =& \frac{2\sqrt{2m_e\mathcal{E}}}{\hslash \sqrt{\gamma (\mathcal{E})}},{\eta }_i:=\frac{\sqrt{\psi }}{m_i^{*}}.\hfill \end{array}$$

Similarly, for phonons that are used to describe the crystal vibrations and strongly influence charge transport, one population for each branch has to be considered [23]. The number of branches is equal to $3\nu$, ν being the total number of atoms per unit cell of the crystal lattice, and a is the lattice constant. Among the most used analytic approximations for the phonon dispersion relation, we report the one which is isotropic and quadratic in the phonon wave vector $\mathbf{q}$

$${ϵ}_p={ϵ}_0^p+\hslash v_s^p\left|\mathbf{q}\right|+\hslash c^p{\left|\mathbf{q}\right|}^2,\left|\mathbf{q}\right|\in \left[0,\frac{2\pi }{a}\right],$$

(2)

where ${ϵ}_p$ is the phonon energy, $v_s^p$ and $c^p$, $p=1,\cdots ,3\nu$, are suitable constants depending on the material under consideration, and p varies over the $3\nu$ phonon branches. The Einstein and the Debye approximations are particular cases of (2), respectively obtained for $v_s^p=c^p=0$ and $c^p=0,v_s^p\ne 0$.

At the kinetic level, the state of the electrons, holes and phonons can be described by their one-particle distribution functions $f_{A_{\beta }}$, $\beta =e$ (electron), h (hole), $A_{\beta }$ indicating the valleys occupied by β-carriers ($\left\{A_e\right\}\cup \left\{A_h\right\}=\left\{A\right\}$), and $g_p$, with $p=1,\cdots ,3\nu$, whose time evolution is determined by the Boltzmann–Bloch–Peierls (BBP) system (see [24] and references therein)

$$\frac{\partial f_{A_{\beta }}}{\partial t}+{\mathbf{v}}_{A_{\beta }}·{\nabla }_{\mathbf{x}}{f}_{A_{\beta }}+\frac{q_{\beta }}{\hslash }\mathbf{E}·{\nabla }_{\mathbf{k}}{f}_{A_{\beta }}={\mathcal{C}}_{im}(f_{A_{\beta }})+{\sum }_{A_{{\beta }^{\prime }}^{\prime },p}{\mathcal{C}}_p^{A_{\beta }{A}_{{\beta }^{\prime }}^{\prime }}(f_{A_{\beta }},f_{A_{{\beta }^{\prime }}^{\prime }},g_p)+{\sum }_{A_{\overline{\beta }}^{\prime },p}{\mathcal{I}}_p^{A_{\beta }{A}_{\overline{\beta }}^{\prime }}[f_{A_{\beta }},f_{A_{\overline{\beta }}^{\prime },g_p}],$$

(3A)

$$\frac{\partial g_p}{\partial t}+{\mathbf{v}}_p·{\nabla }_{\mathbf{x}}{g}_p=\sum _{A_{\beta },A_{{\beta }^{\prime }}^{\prime }}{\mathcal{C}}_{A_{\beta },A_{{\beta }^{\prime }}^{\prime }}^p(g_p,f_{A_{\beta }},f_{A_{{\beta }^{\prime }}^{\prime }})+\sum _{\eta }{\mathcal{C}}_{\eta }^p(g_p)+\sum _{A_{\beta },A_{\overline{\beta }}^{\prime }}{\mathcal{I}}_{A_{\beta }{A}_{\overline{\beta }}^{\prime }}^p[g_p,f_{A_{\beta }},f_{A_{\overline{\beta }}^{\prime }}],$$

(3B)

$$\mathrm{\Delta }({\epsilon }_{s}V)=q\left(n(\mathbf{x},t)-p(\mathbf{x},t)+N_a(\mathbf{x})-N_d(\mathbf{x})\right),$$

(3C)

where $q_{\beta }$ is the charge of the particles, $q_{e,h}=\mp q$ with q absolute value of the elementary charge, $\overline{\beta }$ is the complement of β with respect to the set $\{e,h\}$ ${\mathbf{v}}_p=\frac{1}{\hslash }{\nabla }_{\mathbf{q}}{ϵ}_p$ is the phonon group velocity, ${\epsilon }_s$ is the material permittivity, V and $\mathbf{E}$ are the electric potential and field, η labels the type of scattering among phonons themselves (see below), $N_a$ and $N_d$ are the acceptor and donor concentrations, while n is the total electron density

$$n=\sum _{A_e}{f}_{A_e},$$

(4)

and p the total hole density

$$p=\sum _{A_h}{f}_{A_h}.$$

(5)

The Boltzmann Equation (3A) is coupled among them through the Poisson Equation (3C) and some of the collision operators which appear at the right-hand side of (3A) and (3B).

The operator ${\mathcal{C}}_{im}$ [1] takes into account scatterings between carriers and impurities and is elastic and intravalley, which means that the initial and the final state of the carrier lie in the same valley. Its form is

$${\mathcal{C}}_{im}(f_{A_{\beta }})={\int }_{{\mathbb{R}}^3}[w_{im}({\mathbf{k}}^{\prime },\mathbf{k})f_{A_{\beta }}({\mathbf{k}}^{\prime })-w_{im}(\mathbf{k},{\mathbf{k}}^{\prime })f_{A_{\beta }}(\mathbf{k})]d{\mathbf{k}}^{\prime },$$

the impurity scattering transition rate being given by

$$\begin{array}{ccc}\hfill w_{im}(\mathbf{k},{\mathbf{k}}^{\prime })& =& {\mathcal{K}}^{(im)}\frac{1}{{\left[|\mathbf{k}-{\mathbf{k}}^{\prime }{|}^2+{\lambda }_D^2\right]}^2}\delta (\mathcal{E}({\mathbf{k}}^{\prime })-\mathcal{E}(\mathbf{k})),\hfill \end{array}$$

with ${\lambda }_D=\sqrt{\frac{(N_{a/d})q^2}{{ϵ}_s{k}_B{T}_L}}$ the inverse Debye length and ${\mathcal{K}}^{(im)}=\frac{Z_{a/d}^2{N}_{a/d}{q}^4}{4\pi \hslash {ϵ}_s^2}$, where Z is the impurity charge number, $T_L$ the lattice temperature and $k_B$ the Boltzmann constant. According to whether interactions with donors or acceptors are considered, $N_d$ or $N_a$, and $Z_d$ or $Z_a$ have to be taken.

The collision operators ${\mathcal{C}}_p^{A{A}^{\prime }}(f_{A,f_{A^{\prime }},g_p)}$ describe interactions between carriers of the same type and phonons. These scatterings can be intravalley (intraband) ($A=A^{\prime }$), or intervalley (intraband or interband) ($A\ne A^{\prime }$) [1], which means that the initial and final state belong to different valleys. These operators read [24]

$$\begin{array}{c}{\mathcal{C}}_p^{A{A}^{\prime }}(f_A,f_{A^{\prime }},g_p)={\int }_{S_{\frac{2\pi }{a}}}{\int }_{{\mathbb{R}}^3}[w_p^{A{A}^{\prime }+}(\mathbf{k},{\mathbf{k}}^{\prime },\mathbf{q}){\kappa }_1^{A{A}^{\prime }}(g_p,f_A,f_{A^{\prime }})\hfill \\ +w_p^{A_{\beta }{A}^{\prime }-}(\mathbf{k},{\mathbf{k}}^{\prime },\mathbf{q}){\kappa }_2^{A{A}^{\prime }}(g_p,f_A,f_{A^{\prime }})]d{\mathbf{k}}^{\prime }d\mathbf{q},\hfill \end{array}$$

where $S_{\frac{2\pi }{a}}$ is the sphere of radius $\frac{2\pi }{a}$, which approximates the first Brillouin zone $\mathcal{B}$ in the case of isotropic approximations for the phonon dispersion relations, and, for example, the integrands are given by

$$\begin{array}{c}{\kappa }_1^{A_e{A}_e^{\prime }}(g_p,f_{A_e},f_{A_e^{\prime }}):=\left(\frac{g_p}{y_p}+1\right)f_{A_e^{\prime }}({\mathbf{k}}^{\prime })-\frac{g_p}{y_p}{f}_{A_e}(\mathbf{k}),\hfill \\ {\kappa }_2^{A_e{A}_e^{\prime }}(g_p,f_{A_e},f_{A_e^{\prime }}):=\frac{g_p}{y_p}{f}_{A_e^{\prime }}({\mathbf{k}}^{\prime })-\left(\frac{g_p}{y_p}+1\right)f_{A_e(\mathbf{k})},\hfill \\ w_p^{A_e{A}_e^{\prime }\pm }(\mathbf{k},{\mathbf{k}}^{\prime },\mathbf{q}):=s_p^{A_e{A}_e^{\prime }}(\mathbf{q})\delta \left({\mathcal{E}}_{A_e^{\prime }}({\mathbf{k}}^{\prime })-{\mathcal{E}}_{A_e}(\mathbf{k})\mp {ϵ}_p\right)\delta \left({\mathbf{k}}^{\prime }-\mathbf{k}\mp \mathbf{q}+\mathbf{G}\right),\hfill \end{array}$$

$y_p$ being the p-phonon density of states, δ the Dirac delta function, and $\mathbf{G}$ a vector of the Brillouin zone, whose presence is due to the fact that the total wave vector is conserved up to it. The scattering functions $s_p^{A{A}^{\prime }}$ are characteristic of the type of interaction of carriers, for example with acoustic and non-polar optical phonons for elemental semiconductors such as Si and Ge, and also with polar optical phonons for compound semiconductors, such as GaAs and SiC. The previous expressions of the scattering rates w are valid in the case when only conduction electrons are involved; for all the other cases as well as for the expressions of the scattering functions, we refer the interested reader to [1,8,9].

The generation-recombination collision operators ${{\mathcal{I}}_p}^{A_{\beta }{A}_{\overline{\beta }}^{\prime }}[f_{A_{\beta }},f_{A_{\overline{\beta }}^{\prime }},g_p]$ include several mechanisms, among which the most important are the Auger and the Schockley–Read–Hall (SRH) processes. In a relaxation time approximation and in the simpler case of a single conduction and a single valence band, the operators read [25]

$$\mathcal{I}[f_A,f_{\overline{A}}]=-\left({\Gamma }_{n}n+{\Gamma }_{h}p\right)\left(\tilde{n}{f}_A-{n_i}^2{\mathcal{M}}_A\right)-\frac{\tilde{n}{f}_A-{n_i}^2{\mathcal{M}}_A}{{\tau }_h(n+n_i)+{\tau }_n(p+n_i)},$$

(6)

where $\tilde{n}=n$ or p according to whether the A-th valley is populated by electrons or holes, the ${\mathcal{M}}_A$’s are the Maxwellians normalized to unit, ${\Gamma }_n$ and ${\Gamma }_h$ are the Auger electron and hole constant coefficients, ${\tau }_n$ and ${\tau }_h$ are the carrier lifetimes, and $n_i$ is the intrinsic concentration [1].

The collision operators which appear in the Boltzmann–Peierls equations for phonons, relative to their interactions with carriers, read [24]

$$\begin{array}{c}\hfill {\mathcal{C}}_{A_{\beta }{A}_{{\beta }^{\prime }}^{\prime }}^p(g_p,f_{A_{\beta }},f_{A_{{\beta }^{\prime }}^{\prime }})={\int }_{{\mathbb{R}}^3}{\int }_{{\mathbb{R}}^3}{w}_p^{A_{\beta }{A}_{{\beta }^{\prime }}^{\prime }+}(\mathbf{k},{\mathbf{k}}^{\prime },\mathbf{q}){\kappa }_1^{A_{\beta }{A}_{{\beta }^{\prime }}^{\prime }}(g_p,f_{A_{\beta }},f_{A_{{\beta }^{\prime }}^{\prime }})d{\mathbf{k}}^{\prime }d\mathbf{k}.\end{array}$$

The other interactions of phonons can be distinguished into intrinsic ones, arising from the anharmonicity of the interatomic forces, and extrinsic ones, due to phonon scatterings at the boundaries of the crystal and at various types of crystal defects and imperfections. In their turn, the anharmonic scatterings can be normal processes (N-processes), in which the phonon total momentum after the interaction is conserved, and umklapp processes (U-processes) for which the total momentum changes by a reciprocal lattice vector multiplied by ℏ after a collision. On the other hand, all extrinsic processes do not conserve the total momentum after a collision and, together with the umklapp ones, are usually called resistive processes. For the expressions of the corresponding collision operators, we refer the interested readers to [23] and references therein.

Macroscopic models can be constructed starting from the BBP system by taking a suitable number of moments of the distribution functions [14]. The most common choice is that of considering the following functions of the carrier and phonon wave vectors ${\psi }_{A_{\beta }}(\mathbf{k}):=\left(1,{\mathbf{v}}_{A_{\beta }},{\mathcal{E}}_{A_{\beta }},{\mathcal{E}}_{A_{\beta }}{\mathbf{v}}_{A_{\beta }}\right)$ and $\left({\psi }_p(\mathbf{q})\right):=\left({ϵ}_p,{ϵ}_p{\mathbf{v}}_p\right)$, to which the following macroscopic state variables correspond:

$$\left(\begin{array}{c}{n}_{A_{\beta }}\\ W_{A_{\beta }}\end{array}\right)={\int }_{{\mathbb{R}}^3}\left(\begin{array}{c}1\\ \frac{{\mathcal{E}}_{A_{\beta }}}{n_{A_{\beta }}}\end{array}\right)f_{A_{\beta }}d\mathbf{k},\left(\begin{array}{c}{\mathbf{V}}_{A_{\beta }}\\ {\mathbf{S}}_{A_{\beta }}\end{array}\right)=\frac{1}{n_{A_{\beta }}}{\int }_{{\mathbb{R}}^3}\left(\begin{array}{c}{\mathbf{v}}_{A_{\beta }}\\ {\mathcal{E}}_{A_{\beta }}{\mathbf{v}}_{A_{\beta }}\end{array}\right)f_{A_{\beta }}d\mathbf{k},$$

(7)

$$\left(\begin{array}{c}{W}_p\\ {\mathbf{Q}}_p\end{array}\right)={\int }_{S_{{}_{\frac{2\pi }{a}}}}\left(\begin{array}{c}{ϵ}_p\\ {ϵ}_p{\mathbf{v}}_p\end{array}\right)g_{p}d\mathbf{q},$$

(8)

which are the carrier number densities, the average energies, velocities and energy-fluxes per carrier, and the phonon average energies and energy-fluxes, respectively. The above choice involves the minimal number of moments necessary for describing the thermal energy transport, but this number, if required by the physical problem under study, can be easily extended to cover, for example, an arbitrary number of scalar and vector moments both for carriers and phonons, by taking into account, e.g., higher microscopic energy powers in the functions ψ [7,26].

The evolution equations for the state variables (7) and (8) can be obtained directly from the Boltzmann equations by integration:

$$\frac{\partial }{\partial t}\left(\begin{array}{c}{n}_{A_{\beta }}\\ n_{A_{\beta }}{W}_{A_{\beta }}\end{array}\right)+{\nabla }_{\mathbf{x}}·\left(n_{A_{\beta }}\left(\begin{array}{c}{\mathbf{V}}_{A_{\beta }}\\ {\mathbf{S}}_{A_{\beta }}\end{array}\right)\right)-n_{A_{\beta }}{q}_{\beta }\left(\begin{array}{c}0\\ \mathbf{E}·{\mathbf{V}}_{A_{\beta }}\end{array}\right)=n_{A_{\beta }}\left(\begin{array}{c}{C}_{n_{A_{\beta }}}\\ C_{W_{A_{\beta }}}\end{array}\right),$$

(9)

$$\begin{array}{ccc}\frac{\partial }{\partial t}\left(n_{A_{\beta }}\left(\begin{array}{c}{\mathbf{V}}_{A_{\beta }}\\ {\mathbf{S}}_{A_{\beta }}\end{array}\right)\right)& +& {\nabla }_{\mathbf{x}}\left(n_{A_{\beta }}\left(\begin{array}{c}{F}_{A_{\beta }}^{(0)}\\ F_{A_{\beta }}^{(1)}\end{array}\right)\right)\hfill \\ -n_{A_{\beta }}{q}_{\beta }\left(\begin{array}{c}{G}_{A_{\beta }}^{(0)}\\ G_{A_{\beta }}^{(1)}\end{array}\right)\mathbf{E}& =& n_{A_{\beta }}\left(\begin{array}{c}{C}_{{\mathbf{V}}_{A_{\beta }}}\\ C_{{\mathbf{S}}_{A_{\beta }}}\end{array}\right)\hfill \end{array}$$

(10)

$$\frac{\partial }{\partial t}\left(\begin{array}{c}{W}_p\\ {\mathbf{Q}}_p\end{array}\right)+{\nabla }_{\mathbf{x}}·\left(\begin{array}{c}{\mathbf{Q}}_p\\ T_p\end{array}\right)=\left(\begin{array}{c}{C}_{W_p}\\ C_{{\mathbf{Q}}_p}\end{array}\right).$$

(11)

In the above equations, the extra-fluxes and production terms relative to charges respectively read

$$\begin{array}{c}\left(\begin{array}{c}{F}_{A_{\beta }}^{(0)}\\ F_{A_{\beta }}^{(1)}\end{array}\right)=\frac{1}{n_{A_{\beta }}}\int \left(\begin{array}{c}1\\ {\mathcal{E}}_{A_{\beta }}\end{array}\right){\mathbf{v}}_{A_{\beta }}\otimes {\mathbf{v}}_{A_{\beta }}{f}_{A_{\beta }}d\mathbf{k},\left(\begin{array}{c}\text{velocity flux},\\ \text{flux of the energy-flux}\end{array}\right)\hfill \\ \left(\begin{array}{c}{G}_{A_{\beta }}^{(0)}\\ G_{A_{\beta }}^{(1)}\end{array}\right)=\frac{1}{\hslash n_{A_{\beta }}}\int \left(\begin{array}{c}{\nabla }_{\mathbf{k}}{\mathbf{v}}_{A_{\beta }}\\ {\nabla }_{\mathbf{k}}\left({\mathcal{E}}_{A_{\beta }}{\mathbf{v}}_{A_{\beta }}\right)\end{array}\right)f_{A_{\beta }}d\mathbf{k},\hfill \end{array}$$

$$\begin{array}{cc}& n_{A_{\beta }}{C}_{M_{{\psi }_{A_{\beta }}}}=\int {\psi }_{A_{\beta }}(\mathbf{k})[{\mathcal{C}}_{im}(f_{A_{\beta }})+\sum _{A_{{\beta }^{\prime }}^{\prime },p}{\mathcal{C}}_p^{A_{\beta }{A}_{{\beta }^{\prime }}^{\prime }}(f_{A_{\beta }},f_{A_{{\beta }^{\prime }}^{\prime },},g_p)\hfill \\ +& \sum _{A_{\overline{\beta }}^{\prime },p}{\mathcal{I}}_p^{A_{\beta }{A}_{\overline{\beta }}^{\prime }}[f_{A_{\beta }},f_{A_{\overline{\beta }}^{\prime },g_p}]]d\mathbf{k},M_{{\psi }_{A_{\beta }}}\text{-production},\hfill \end{array}$$

with $M_{{\psi }_{A_{\beta }}}:=\left(n_{A_{\beta }},{\mathbf{V}}_{A_{\beta }},W_{A_{\beta }},{\mathbf{S}}_{A_{\beta }}\right)$, while for the phonons, we have

$$\begin{array}{c}{T}_p=\int {ϵ}_p{\mathbf{v}}_p\otimes {\mathbf{v}}_p{g}_p(\mathbf{q})d\mathbf{q},\mathrm{flux}\mathrm{of}\mathrm{the} \text{energy-flux},\hfill \\ \left(\begin{array}{c}{C}_{W_p}\\ C_{{\mathbf{Q}}_p}\end{array}\right)=\int \left(\begin{array}{c}{ϵ}_p\\ {ϵ}_p{\mathbf{v}}_p\end{array}\right)\sum _{A_{\beta },A_{{\beta }^{\prime }}^{\prime }}{\mathcal{C}}_{A_{\beta },A_{{\beta }^{\prime }}^{\prime }}^p(g_p,f_{A_{\beta }},f_{A_{{\beta }^{\prime }}^{\prime }})+\sum _{\eta }{\mathcal{C}}_{\eta }^p\hfill \\ +\sum _{A_{\beta },A_{\overline{\beta }}^{\prime }}{\mathcal{I}}_{A_{\beta }{A}_{\overline{\beta }}^{\prime }}^p[g_p,f_{A_{\beta }},f_{A_{\overline{\beta }}^{\prime }}]d\mathbf{q},\left(\begin{array}{c}\mathrm{energy}\mathrm{production}\\ \text{energy-flux production}\end{array}\right).\hfill \end{array}$$

(12)

In the evolution equations, the number of the unknowns is greater than that of the equations, therefore constitutive equations are needed for the extra-variables $F_{ij}^{A(0)},G_{ij}^{A(0)},F_{ij}^{A(1)},$ $G_{ij}^{A(1)},T_{ij}^p$, $C_{n_A},C_{W_A},C_{{\mathbf{V}}_A},C_{{\mathbf{S}}_A}$, $C_{W_p},$ $C_{{\mathbf{Q}}_p}$. A systematic way to find these relations is founded on a universal physical principle: the maximum entropy principle [6,14,27,28]. It states that, if a certain number of moments is known, then the least biased distribution functions, which can be used for evaluating the unknown moments, are those maximizing the total entropy functional under the constraint that they reproduce the known moments. In the case under consideration, neglecting the mutual interactions among the subsystems and degeneration effects (the degeneracy case can be tackled in the same way; we present the non-degenerate case only for the sake of simplicity), the total entropy is

$$\begin{array}{c}\mathcal{S}=-k_B\{\sum _{A_{\beta }}{\int }_{{\mathbb{R}}^3}\left(f_{A_{\beta }}\ln \frac{f_{A_{\beta }}}{y_{A_{\beta }}}-f_{A_{\beta }}\right)d\mathbf{k}+\sum _p{\int }_{S_{\frac{2\pi }{a}}}[g_p\ln \frac{g_p}{y_p}\hfill \\ -y_p\left(1+\frac{g_p}{y_p}\right)\ln \left(1+\frac{g_p}{y_p}\right)\left]d\mathbf{q}\right\},\hfill \end{array}$$

with $y_{A_{\beta }}=\frac{2}{{(2\pi )}^3}$ the charge density of states, while the constraints are given by (7) and (8). Let us introduce the functional spaces

$${\mathcal{F}}_{{\psi }_{A_{\beta }}}=\left\{h:{\mathbb{R}}^3\mapsto \mathbb{R}s.t.{\int }_{{\mathbb{R}}^3}\left|h(\mathbf{k}){\psi }_j(\mathbf{k})\right|d\mathbf{k}<+\infty ,\forall {\psi }_j\in \{{\psi }_{A_{\beta }}\}\right\},$$

$${\mathcal{F}}_{{\psi }_p}=\left\{h:{\mathbb{R}}^3\mapsto \mathbb{R}s.t.{\int }_{S_{\frac{2\pi }{a}}}\left|h(\mathbf{q}){\psi }_j(\mathbf{q})\right|d\mathbf{q}<+\infty ,\forall {\psi }_j\in \{{\psi }_p\}\right\}.$$

Given the values of the moments $M_{{\psi }_{A_{\beta }}}$ and $W_p,{\mathbf{Q}}_p$, MEP amounts to solve the following optimization problem:

$$\underset{f_{A_{\beta }}\in {\mathcal{F}}_{{\psi }_{A_{\beta }}},g_p\in {\mathcal{F}}_{{\psi }_p}}{max}\mathcal{S}$$

under the constraints

$${\int }_{{\mathbb{R}}^3}{\psi }_{A_{\beta }}(\mathbf{k})f_{A_{\beta }}d\mathbf{k}=M_{{\psi }_{A_{\beta }}},{\int }_{S_{\frac{2\pi }{a}}}{g}_p\left(\begin{array}{c}{ϵ}_p\\ {ϵ}_p{\mathbf{v}}_p\end{array}\right)d\mathbf{q}=\left(\begin{array}{c}{W}_P\\ {\mathbf{Q}}_p\end{array}\right)$$

where time and position are considered fixed and omitted for the sake of simplifying the notation. The solution of this maximization problem is given by [29]

$$\begin{array}{ccc}\hfill f_{A_{\beta }}^{ME}& =& \exp (-{\Lambda }_{A_{\beta }}-{\Lambda }_{W_{A_{\beta }}}{\mathcal{E}}_{A_{\beta }}-{\mathbf{v}}_{A_{\beta }}·({\Lambda }_{{\mathbf{V}}_{A_{\beta }}}+{\mathcal{E}}_{A_{\beta }}{\Lambda }_{{\mathbf{S}}_{A_{\beta }}})),\hfill \\ \hfill g_p^{ME}& =& {\displaystyle \frac{y_p}{\exp ({ϵ}_p{\Lambda }_{W_p+{ϵ}_p{\mathbf{v}}_p·{\Lambda }_{{\mathbf{Q}}_p}})-1},}\hfill \end{array}$$

which, linearized with respect to the vector variables (small anisotropy assumption) [8,9,27], becomes

$$\begin{array}{ccc}\hfill f_{A_{\beta }}^{ME}& =& \exp (-{\Lambda }_{A_{\beta }}-{\Lambda }_{W_{A_{\beta }}}{\mathcal{E}}_{A_{\beta }})\left(1-{\mathbf{v}}_{A_{\beta }}·({\Lambda }_{{\mathbf{V}}_{A_{\beta }}}+{\mathcal{E}}_{A_{\beta }}{\Lambda }_{{\mathbf{S}}_{A_{\beta }}})\right),\hfill \\ \hfill g_p^{ME}& =& {\displaystyle \frac{y_p}{\exp ({ϵ}_p{\Lambda }_{W_p})-1}-\frac{y_p{ϵ}_p\exp ({ϵ}_p{\Lambda }_{{}_{W_p}})}{{(\exp ({ϵ}_p{\Lambda }_{W_p})-1)}^2}{\mathbf{v}}_p·{\Lambda }_{{\mathbf{Q}}_p},}\hfill \end{array}$$

where the Λ’s are Lagrange multipliers, related to the state variables by means of the constraint relations (7) and (8). Linearization is made in order to simplify the inversion of the constraints, otherwise it has to be done numerically [29]. After inversion, the dependence of the distribution functions on $(\mathbf{x},t)$ will only be through the state variables, and substituting the distributions into the integrals defining the extra-variables, the needed closure relations can be obtained. It can be proven that the balance equations with the closure relations given by MEP form a quasilinear hyperbolic system in the relevant physical region of the field variables.

Following this procedure, mathematical models for the description of charge transport in unipolar and bipolar silicon devices have been constructed, see for example [10,30], while results relative to Si thermal properties can be found in [23,31]. The procedure has also been applied to compound semiconductors such as GaAs, GaN and SiC [8,9].

3. Quantum Corrections to the Semiclassical Models

Based on the previous considerations, a natural way to get a quantum macroscopic model is to use MEP in a quantum framework to close the moment system arising from the Wigner equations. The general guidelines can be found in [32,33]. This approach has been followed in [34] (see also [35] and references therein). However, one has to deal with operatorial equations which are very complex and therefore hard to be solved numerically. Moreover, drastic simplifications of the collision terms are introduced in order to make the problem tractable.

Another strategy can be adopted to close the moment system arising from the Wigner equation. According to this strategy, the Wigner function and the relative transport equation are expanded with respect to ${\hslash }^2$. The zero-order part of the collision operator is supposed to be the same as the semiclassical one, while the first-order contribution is supposed to act only on the ${\hslash }^2$ correction of the Wigner function and is modeled in a relaxation form. Therefore, at the zero order, the Wigner function is given by the solution of the semiclassical Boltzmann equation, which is approximated with the standard maximum entropy method, while the ${\hslash }^2$ order correction is obtained with a Chapman–Enskog expansion in the high field scaling.

The typical physical situation which can be described in this way is that when the main contribution to the charge transport is semiclassical, the quantum effects enter as small perturbations. For example, this is reasonable for devices such as MOSFETs of characteristic length of about ten nanometers subjected to strong electric fields.

The main assumption is that there is a balance between the ${\hslash }^2$ drift and collision terms. This can be motivated by the fact that the collision frequency of the semiclassical scatterings tends to increase as the energy rises. Therefore, quantum effects are expected to be relevant at high fields; in such conditions, there should be high energies with consequently high collision frequencies.

The starting point for our derivation of the quantum corrections to the semiclassical model is the single particle Wigner–Poisson system, which represents the quantum analogue of the semiclassical Boltzmann–Poisson system (for the sake of simplicity, only the case of a single valley in the conduction band is considered),

$$\frac{\partial w}{\partial t}+S[\mathcal{E}]w+\frac{q}{m^{*}}\Theta [V]w=\mathcal{C}[w],$$

(13)

$$\mathrm{div}\left({ϵ}_s{\nabla }_{\mathbf{x}}{V}_S\right)=-q(N_d-N_a-n).$$

(14)

where $m^{*}$ is the electron effective mass, the potential V usually is the sum of a self-consistent part $V_S$, solution of the Poisson equation, and an additional part which models the potential barriers in heterostructures. The unknown function $w(\mathbf{x},\mathbf{p},t)$, depending on the position $\mathbf{x}$, crystal momentum $\mathbf{p}=\hslash \mathbf{k}$ and time t, is the Wigner quasi-distribution, defined as

$$\begin{array}{cc}\hfill w(\mathbf{x},\mathbf{p},t)=& {\mathcal{F}}^{-1}[\rho (\mathbf{x}+\frac{\hslash }{2m^{*}}\mathit{\eta },\mathbf{x}-\frac{\hslash }{2m^{*}}\mathit{\eta },t)](\mathbf{p})\hfill \\ & =\frac{1}{{(2\pi )}^3}{\int }_{{\mathbb{R}}^3}\rho (\mathbf{x}+\frac{\hslash }{2m^{*}}\mathit{\eta },\mathbf{x}-\frac{\hslash }{2m^{*}}\mathit{\eta },t)e^{i\mathbf{p}·\mathit{\eta }}d\mathit{\eta }.\hfill \end{array}$$

(15)

Here, $\rho (\mathbf{x},\mathbf{y},t)$ is the density matrix, which is related to the wave function $\psi (\mathbf{x},t)$ by

$$\begin{array}{c}\hfill \rho (\mathbf{x},\mathbf{y},t)=\overline{\psi (\mathbf{x},t)}\psi (\mathbf{y},t).\end{array}$$

(16)

The operator $\mathcal{F}$ denotes the Fourier transform, given for a function $g(\mathbf{p})\in L^1({\mathbb{R}}^3)$ by

$$\mathcal{F}[g(\mathbf{p})](\mathit{\eta })={\int }_{{\mathbb{R}}_{\mathbf{p}}^3}g(\mathbf{p})e^{-i\mathbf{p}·\mathit{\eta }}d\mathbf{p},$$

with i the imaginary unit, and ${\mathcal{F}}^{-1}$ denotes the inverse Fourier transform

$${\mathcal{F}}^{-1}[h(\mathit{\eta })]=\frac{1}{{(2\pi )}^3}{\int }_{{\mathbb{R}}_{\mathit{\eta }}^3}h(\mathit{\eta })e^{i\mathbf{p}·\mathit{\eta }}d\mathit{\eta }.$$

The terms $S[\mathcal{E}]$ and $\Theta [V]$ represent the pseudo-differential operators

$$\begin{array}{ccc}\hfill S[\mathcal{E}]w(\mathbf{x},\mathbf{p},t)& =& \frac{i}{\hslash {(2\pi )}^3}{\int }_{{\mathbb{R}}_{{\mathbf{x}}^{\prime }}^3\times {\mathbb{R}}_{\nu }^3}\left[\mathcal{E}\left(\mathbf{p}+\frac{\hslash }{2}\mathit{\nu },t\right)-\mathcal{E}\left(\mathbf{p}+\frac{\hslash }{2}\mathit{\nu },t\right)\right]\hfill \\ & & w({\mathbf{x}}^{\prime },\mathbf{p},t)e^{-i({\mathbf{x}}^{\prime }-\mathbf{x})·\mathit{\nu }}d{\mathbf{x}}^{\prime }d\mathit{\nu }\hfill \\ \hfill \Theta [V]w(\mathbf{x},\mathbf{p},t)& =& \frac{i{m}^{*}}{\hslash {(2\pi )}^3}{\int }_{{\mathbb{R}}_{{\mathbf{p}}^{\prime }}^3\times {\mathbb{R}}_{\mathit{\eta }}^3}\left[V\left(\mathbf{x}+\frac{\hslash }{2m^{*}}\mathit{\eta },t\right)-V\left(\mathbf{x}-\frac{\hslash }{2m^{*}}\mathit{\eta },t\right)\right]\hfill \\ & & w(\mathbf{x},{\mathbf{p}}^{\prime },t)e^{-i({\mathbf{p}}^{\prime }-\mathbf{p})·\mathit{\eta }}d{\mathbf{p}}^{\prime }d\mathit{\eta }.\hfill \end{array}$$

As is well known, $w(\mathbf{x},\mathbf{p},t)$ is not in general positive definite. However, it is possible to calculate the macroscopic quantities of interest as expectation values (moments) of $w(\mathbf{x},\mathbf{p},t)$ in the same way as in the semiclassical case, e.g.,

$$n(\mathbf{x},t)={\int }_{{\mathbb{R}}^3}w(\mathbf{x},\mathbf{p},t)d\mathbf{p}(\mathrm{density}),$$

(17)

$$\mathbf{V}(\mathbf{x},t)=\frac{1}{n(\mathbf{x},t)}{\int }_{{\mathbb{R}}^3}\mathbf{v}w(\mathbf{x},\mathbf{p},t)d\mathbf{p}(\mathrm{velocity}),$$

(18)

$$W(\mathbf{x},t)=\frac{1}{n(\mathbf{x},t)}{\int }_{{\mathbb{R}}^3}\mathcal{E}(\mathbf{p})w(\mathbf{x},\mathbf{p},t)d\mathbf{p}(\mathrm{energy}),$$

(19)

$$\mathbf{S}(\mathbf{x},t)=\frac{1}{n(\mathbf{x},t)}{\int }_{{\mathbb{R}}^3}\mathcal{E}(\mathbf{p})\mathbf{v}w(\mathbf{x},\mathbf{p},t)d\mathbf{p}(\text{energy-flux}).$$

(20)

The operator $\mathcal{C}[w]$ represents the quantum collision term. Its formulation is itself an open problem. Some attempts can be found in [36,37], but a derivation suitable for applications in electron devices is still lacking. In [15], an expression has been proposed which is a perturbation of the semiclassical collision term, useful for the formulation of macroscopic models.

As a general guideline, $\mathcal{C}[w]$ should drive the system towards the equilibrium. Let us consider the electrons in a thermal bath at the lattice temperature $T_L=1/k_B\beta$. The equilibrium Wigner function $w_{eq}$ for an arbitrary energy band, up to ${\hslash }^2$ terms, reads [38,39]

$$\begin{array}{c}\hfill w_{eq}(\mathbf{x},\mathbf{p},\beta )=\frac{n(\mathbf{x},t)e^{-\beta \mathcal{E}}}{A_0(\beta ,m^{*})}\left\{1+{\hslash }^2\left[\left(\frac{q{\beta }^2}{8}\frac{{\partial }^{2}V}{\partial x_i\partial x_j}+\frac{q^2{\beta }^3}{24}\frac{\partial V}{\partial x_i}\frac{\partial V}{\partial x_j}\right)\right.\right.\\ \hfill \left.\left.\left(\frac{{\partial }^2\mathcal{E}}{\partial p_i\partial p_j}-\frac{A_{ij}(\beta ,m^{*})}{A_0(\beta ,m^{*})}\right)-\frac{q{\beta }^3}{24}\frac{{\partial }^{2}V}{\partial x_i\partial x_j}\left(v_i{v}_j-\frac{B_{ij}(\beta ,m^{*})}{A_0(\beta ,m^{*})}\right)\right]\right\}+\mathcal{O}({\hslash }^4),\end{array}$$

where

$$\begin{array}{ccc}\hfill A_0(\beta ,m^{*})& =& {\int }_{{\mathbb{R}}^3}{e}^{-\beta \mathcal{E}}d\mathbf{p},A_{ij}(\beta ,m^{*})={\int }_{{\mathbb{R}}^3}{e}^{-\beta \mathcal{E}}\frac{{\partial }^2\mathcal{E}}{\partial p_i\partial p_j}d\mathbf{p},\hfill \\ \hfill B_{ij}(\beta ,m^{*})& =& {\int }_{{\mathbb{R}}^3}{e}^{-\beta \mathcal{E}}{v}_i{v}_{j}d\mathbf{p}.\hfill \end{array}$$

We suppose that the expansion

$$w=w^{(0)}+{\hslash }^2{w}^{(1)}+\mathcal{O}({\hslash }^4)$$

(21)

holds. By proceeding in a formal way, as $\hslash \to 0$ the Wigner equation gives the semiclassical Boltzmann equation

$$\frac{\partial w^{(0)}}{\partial t}+\mathbf{v}·{\nabla }_{\mathbf{x}}{w}^{(0)}+\frac{q}{m^{*}}{\nabla }_{\mathbf{x}}V·{\nabla }_{\mathbf{v}}{w}^{(0)}={\mathcal{C}}^{(0)}[w^{(0)}]$$

(22)

While, at the first order in ${\hslash }^2$, we have (Einstein’s convention is used: summation with respect to repeated dummy indices is understood.)

$$\begin{array}{ccc}{\displaystyle \frac{\partial w^{(1)}}{\partial t}}& +& \mathbf{v}·{\nabla }_{\mathbf{x}}{w}^{(1)}-{\displaystyle \frac{1}{24}}{\displaystyle \frac{{\partial }^3\mathcal{E}}{\partial p_i\partial p_j\partial p_k}}{\displaystyle \frac{{\partial }^3{w}^{(1)}}{\partial x_i\partial x_j\partial x_k}}\hfill \\ & +& {\displaystyle \frac{q}{m^{*}}}{\nabla }_{\mathbf{x}}V·{\nabla }_{\mathbf{v}}{w}^{(1)}-{\displaystyle \frac{q}{24m^{*3}}}{\displaystyle \frac{{\partial }^{3}V}{\partial x_i\partial x_j\partial x_k}}{\displaystyle \frac{{\partial }^3{w}^{(0)}}{\partial v_i\partial v_j\partial v_k}}={\mathcal{C}}^{(1)}\hfill \end{array}$$

(23)

with ${\mathcal{C}}^{(1)}$ to be modeled. We make the following

Assumption 1.

$$\begin{array}{ccc}\hfill \mathcal{C}[w]& =& {\mathcal{C}}^{(0)}[w^{(0)}]+{\hslash }^2{\mathcal{C}}^{(1)}[w^{(1)}]={\mathcal{C}}_C[w^{(0)}]-{\hslash }^2\nu \left(w^{(1)}-w_{eq}^{(1)}\right)+\mathcal{O}({\hslash }^4),\hfill \\ & & \mathrm{with}{\mathcal{C}}_C[w^{(0)}]\mathrm{the}\mathrm{semiclassical}\mathrm{collision}\mathrm{operator}\hfill \\ & & \mathrm{and}\nu >0\mathrm{quantum}\mathrm{collision}\mathrm{frequency}\hfill \end{array}$$

(24)

Remark 1.

At variance with other approaches, only the ${\hslash }^2$ correction to the collision term has a relaxation form. This assures that as $\hslash \to 0$ one gets the semiclassical scattering of electrons with phonons and impurities. We note that $w^{(0)}>0$ and therefore the semiclassical expression of the collision term makes sense.

The value of the quantum collision frequency ν is a fitting parameter that can be determined by comparing the results with the experimental data.

We require that $\mathcal{C}[w]$ conserves the electron density, that is we make the following

Basic Assumption 2.

$${\int }_{{\mathbb{R}}^3}\mathcal{C}[w]d\mathbf{p}=0.$$

(25)

Proposition 1.

The collision operator $\mathcal{C}[w]$ of the form (24) satisfies up to terms $\mathcal{O}({\hslash }^4)$, the following properties:

1. 

$Ker\left(\mathcal{C}[w]\right)$ is given by the quantum Maxwellian

$$w_{(eq)}=w_{eq}^{(0)}+{\hslash }^2{w}_{eq}^{(1)},$$

with $w_{eq}^{(0)}$ the classical Maxwellian.

2. 

$$-k_B{\int }_{{\mathbb{R}}^3}{\mathcal{C}}^{(0)}[w^{(0)}]\ln \frac{w^{(0)}}{\exp (-\frac{\beta m^{*}{\left|\mathbf{v}\right|}^2}{2})}d\mathbf{p}=-k_B{\int }_{{\mathbb{R}}^3}\left[\ln w^{(0)}+\frac{\beta m^{*}{\left|\mathbf{v}\right|}^2}{2}\right]{\mathcal{C}}^{(0)}d\mathbf{p}\ge 0,$$

3. 

$$-\frac{1}{2}{\mathcal{C}}^{(1)}[w^{(1)}]\left(w^{(1)}-w_{eq}^{(1)}\right)\ge 0$$

Moreover, the equality holds if and only if w is the quantum Maxwellian, defined above.

Properties 1 and 3 are straightforward. Property 2 is based on the proof in [40,41,42] valid in the classical case.

3.1. Quantum Corrections in the High Field Approximation

In the case of high electric fields and effective mass approximation for the dispersion relation

$$\mathcal{E}=\frac{{\left|\mathbf{p}\right|}^2}{2m^{*}},\mathbf{v}=\frac{\mathbf{p}}{m^{*}},$$

it is possible to find an approximation for $w^{(1)}$ by a suitable Chapman–Enskog expansion. Let us introduce the dimensionless variables

$$\tilde{x}=\frac{\left|\mathbf{x}\right|}{l_0},\tilde{t}=\frac{t}{t_0},\tilde{v}=\frac{\left|\mathbf{v}\right|}{v_0},$$

with $l_0$, $t_0$ and $v_0=l_0/t_0$ the typical length, time and velocity. Let $l_V$ be the characteristic length of the electrical potential and $1/t_C$ the characteristic collision frequency. After scaling the collision frequency according to $\tilde{\nu }=t_C\nu$, Equation (23) can be rewritten as

$$\begin{array}{ccc}\hfill \frac{1}{t_0}\frac{\partial w^{(1)}}{\partial t}+\frac{v_0}{l_0}\mathbf{v}·{\nabla }_{\mathbf{x}}{w}^{(1)}& +& \frac{v_0}{l_V}\left[\frac{q}{m^{*}}{\nabla }_{\mathbf{x}}V·{\nabla }_{\mathbf{v}}{w}^{(1)}-\frac{q}{24m^{*3}}\frac{{\partial }^{3}V}{\partial x_i\partial x_j\partial x_k}\frac{{\partial }^3{w}^{(0)}}{\partial v_i\partial v_j\partial v_k}\right]\hfill \\ & =& -\frac{1}{t_C}\nu \left(w^{(1)}-w_{eq}^{(1)}\right).\hfill \end{array}$$

We will continue to denote the scaled variables as the unscaled ones in order to simplify the notation.

Let us introduce the characteristic length associated with the quantum correction of the collision term (a kind of mean free path in a semiclassical context)

$$l_C=v_0{t}_C.$$

Let us assume that the quantum effects occur in the high field and collision dominated regime, where drift and collision mechanisms have the same characteristic length. Therefore, we set formally

$$\frac{l_C}{l_V}=1$$

and observe that, in the high frequency regime, the Knudsen number

$$\alpha =\frac{l_C}{l_0}$$

is a small parameter. Substituting it in the previous equation, we get

$$\begin{array}{ccc}\hfill \alpha \frac{\partial w^{(1)}}{\partial t}+\alpha \mathbf{v}·{\nabla }_{\mathbf{x}}{w}^{(1)}+\frac{q}{m^{*}}{\nabla }_{\mathbf{x}}V·{\nabla }_{\mathbf{v}}{w}^{(1)}& -& \frac{q}{24m^{*3}}\frac{{\partial }^{3}V}{\partial x_i\partial x_j\partial x_k}\frac{{\partial }^3{w}^{(0)}}{\partial v_i\partial v_j\partial v_k}\hfill \\ & =& -\nu \left(w^{(1)}-w_{eq}^{(1)}\right).\hfill \end{array}$$

(26)

The zero order in α gives

$$\frac{q}{m^{*}}{\nabla }_{\mathbf{x}}V·{\nabla }_{\mathbf{v}}{w}^{(1)}-\frac{q}{24m^{*3}}\frac{{\partial }^{3}V}{\partial x_i\partial x_j\partial x_k}\frac{{\partial }^3{w}^{(0)}}{\partial v_i\partial v_j\partial v_k}=-\nu \left(w^{(1)}-w_{eq}^{(1)}\right),$$

(27)

and by Fourier transforming, one has

$$\begin{array}{c}\hfill w^{(1)}(\mathbf{x},\mathbf{v},t)={\mathcal{F}}^{-1}\left\{\frac{1}{\nu +\frac{iq}{m^{*}}\mathit{\eta }·{\nabla }_{\mathbf{x}}V}\left[-\frac{iq}{24m^{*3}}\frac{{\partial }^{3}V}{\partial x_i\partial x_j\partial x_k}{\eta }_i{\eta }_j{\eta }_k\mathcal{F}{w}^{(0)}(\mathit{\eta })\right.\right.\\ \hfill +\nu \mathcal{F}{w}_{eq}^{(1)}(\mathit{\eta })]\}(\mathbf{x},\mathbf{v},t).\end{array}$$

(28)

Approximating $w^{(0)}$ with MEP distribution function, we obtain

$$\begin{array}{c}\hfill w(\mathbf{x},\mathbf{v},t)\approx f_{ME}(\mathbf{x},\mathbf{v},t)+{\hslash }^2{w}^{(1)}(\mathbf{x},\mathbf{v},t).\end{array}$$

(29)

For example, in the 8-moment case, using the $f_{ME}$ found in [10], the following explicit approximation for the Wigner function is obtained

$$\begin{array}{c}w(\mathbf{x},\mathbf{v},t)\approx {\displaystyle \frac{n\exp (-\frac{3m^{*}{\left|\mathbf{v}\right|}^2}{4W})}{{\left({\displaystyle \frac{4}{3}}\pi m^{*}W\right)}^{3/2}}}\left\{1-\left[\left(-{\displaystyle \frac{21m^{*}}{4W}}+\mathcal{E}{\displaystyle \frac{9m^{*}}{4W^2}}\right){\mathbf{V}}^{(0)}+\left({\displaystyle \frac{9m^{*}}{4W^2}}-\mathcal{E}{\displaystyle \frac{27m^{*}}{20W^3}}\right){\mathbf{S}}^{(0)}\right]·\mathbf{v}\right\}\hfill \\ +{\hslash }^2{\mathcal{F}}^{-1}\left\{{\displaystyle \frac{1}{\nu +{\displaystyle \frac{iq}{m^{*}}}\mathit{\eta }·{\nabla }_{\mathbf{x}}V}}\left[-{\displaystyle \frac{iq}{24m^{*3}}}{\displaystyle \frac{{\partial }^{3}V}{\partial x_i\partial x_j\partial x_k}}{\eta }_i{\eta }_j{\eta }_k\mathcal{F}{w}^{(0)}(\mathit{\eta })+\nu \mathcal{F}{w}_{eq}^{(1)}(\mathit{\eta })\right]\right\}(\mathbf{x},\mathbf{v},t)\hfill \end{array}$$

(30)

which can be used for evaluating the unknown quantities in the moment system, associated with the Wigner equation. The vectors ${\mathbf{V}}^{(0)}$ and ${\mathbf{S}}^{(0)}$ represent the velocity and the energy-flux at zero order in ${\hslash }^2$.

3.2. The Quantum Moment Equations

In analogy with the semiclassical case, multiplying (13) by suitable weight functions ψ, depending, in the physically relevant cases, on the velocity $\mathbf{v}$, and integrating over the velocity space, one has the balance equations for the macroscopic quantities of interest

$$\begin{array}{c}\hfill \frac{\partial }{\partial t}{\int }_{{\mathbb{R}}^3}w(\mathbf{x},\mathbf{v},t)\psi (\mathbf{v})d\mathbf{v}+{\nabla }_{\mathbf{x}}·{\int }_{{\mathbb{R}}^3}\psi (\mathbf{v})\mathbf{v}wd\mathbf{v}+\frac{q}{m^{*}}{\int }_{{\mathbb{R}}^3}\psi (\mathbf{v})\Theta [V]wd\mathbf{v}\\ \hfill ={\int }_{{\mathbb{R}}^3}\psi (\mathbf{v})\mathcal{C}[w]d\mathbf{v}.\end{array}$$

(31)

In the 8-moment model, the basic variables are the macroscopic density, velocity, energy and energy-flux, that are the moments relative to the weight functions 1, $\mathbf{v}$, $\frac{1}{2}{m}^{*}{\left|\mathbf{v}\right|}^2$, $\frac{1}{2}{m}^{*}{\left|\mathbf{v}\right|}^2\mathbf{v}$.

By evaluating (31) for $\psi =1$, under the assumption that the necessary moments of $w^{(1)}(\mathbf{x},\mathbf{v},t)$ and $\frac{{\partial }^3{w}^{(0)}}{\partial v_i{v}_j{v}_k}$ with respect to $\mathbf{v}$ exist, one has

$$\begin{array}{c}\hfill \frac{q}{m^{*}}{\int }_{{\mathbb{R}}^3}\Theta [V]wd\mathbf{v}=\frac{q}{m^{*}}{\nabla }_{\mathbf{x}}·{\int }_{{\mathbb{R}}^3}{\nabla }_{\mathbf{v}}{w}^{(0)}d\mathbf{v}+{\hslash }^2\left[\frac{q}{m^{*}}{\nabla }_{\mathbf{x}}V·{\int }_{{\mathbb{R}}^3}{\nabla }_{\mathbf{v}}{w}^{(1)}d\mathbf{v}\right.\\ \hfill \left.-\frac{q}{24m^3}\frac{{\partial }^{3}V}{\partial x_i\partial x_j\partial x_k}{\int }_{{\mathbb{R}}^3}\frac{{\partial }^3{w}^{(0)}}{\partial v_i\partial v_j\partial v_k}d\mathbf{v}\right]=0,\end{array}$$

obtaining the continuity equation

$$\frac{\partial }{\partial t}n+\frac{\partial (n{V}_i)}{\partial x^i}=0.$$

(32)

In order to get the other moment equations, we observe that from (27) it follows that

$$\begin{array}{c}\hfill {\displaystyle \frac{q}{m^{*}}}{\nabla }_{\mathbf{x}}V·{\int }_{{\mathbb{R}}^3}\psi (\mathbf{v}){\nabla }_{\mathbf{v}}{w}^{(1)}d\mathbf{v}-\frac{q}{24m^3}{\displaystyle \frac{{\partial }^{3}V}{\partial x_i\partial x_j\partial x_k}}{\int }_{{\mathbb{R}}^3}\psi (\mathbf{v}){\displaystyle \frac{{\partial }^3{w}^{(0)}}{\partial v_i\partial v_j\partial v_k}}d\mathbf{v}\\ \hfill +\nu {\int }_{{\mathbb{R}}^3}\psi (\mathbf{v})\left(w^{(1)}-w_{eq}^{(1)}\right)d\mathbf{v}=0,\end{array}$$

(33)

for any weight function $\psi (\mathbf{v})$ such that the integrals exist.

By taking into account (33), multiplying Equation (13) by the weight functions $\mathbf{v}$, $\frac{1}{2}{m}^{*}{\left|\mathbf{v}\right|}^2$, $\frac{1}{2}{m}^{*}{\left|\mathbf{v}\right|}^2\mathbf{v}$, and after integration, one finds the balance equations for velocity, energy and energy-flux

$$\frac{\partial }{\partial t}(n{V}_i)+\frac{\partial (n{U}_{ij})}{\partial x_j}+nq{E}_i=n{C}_{V_i}\left(W^{(0)},V_i^{(0)},S_i^{(0)}\right),$$

(34)

$$\frac{\partial }{\partial t}(nW)+\frac{\partial (n{S}_j)}{\partial x_j}+nq{V}_k^{(0)}{E}^k=n{C}_W(W^{(0)}),$$

(35)

$$\frac{\partial }{\partial t}(n{S}_i)+\frac{\partial (n{F}_{ij})}{\partial x_j}+\frac{5}{3}n\frac{q}{m^{*}}{E}_i{W}^{(0)}=n{C}_{S_i}\left(W^{(0)},V_i^{(0)},S_i^{(0)}\right).$$

(36)

Here, $V_i^{(0)}$, $W^{(0)}$ and $P_i^{(0)}$ are the zero-order components of the average velocity, energy and energy-flux. The components of the flux of the velocity and the flux of the energy-flux are defined as

$$\begin{array}{ccc}\hfill U_{ij}& =& \frac{1}{n(\mathbf{x},t)}{\int }_{{\mathbb{R}}^3}{v}_i{v}_{j}w(\mathbf{x},\mathbf{v},t)d\mathbf{v},\hfill \end{array}$$

(37)

$$\begin{array}{ccc}\hfill F_{ij}& =& \frac{1}{n(\mathbf{x},t)}{\int }_{{\mathbb{R}}^3}\frac{1}{2}{m}^{*}{v}_i{v}_j{\left|\mathbf{v}\right|}^{2}w(\mathbf{x},\mathbf{v},t)d\mathbf{v}.\hfill \end{array}$$

(38)

The production terms are defined as

$$\begin{array}{ccc}\hfill n{C}_{V_i}& =& {\int }_{{\mathbb{R}}^3}{v}_i\mathcal{C}[w]d\mathbf{v},\hfill \end{array}$$

(39)

$$\begin{array}{ccc}\hfill n{C}_W& =& {\int }_{{\mathbb{R}}^3}\frac{1}{2}{m}^{*}{v}^2\mathcal{C}[w]d\mathbf{v},\hfill \end{array}$$

(40)

$$\begin{array}{ccc}\hfill n{C}_{S_i}& =& {\int }_{{\mathbb{R}}^3}\frac{1}{2}{m}^{*}{\left|\mathbf{v}\right|}^2{v}_i\mathcal{C}[w]d\mathbf{v}.\hfill \end{array}$$

(41)

Remark 2.

The quantum corrections affect only the free streaming part, while the drift and production terms appear only at the zero order.

Therefore, $C_W(W^{(0)})$, $C_{V_i}\left(W^{(0)},V_i^{(0)},S_i^{(0)}\right)$ and $C_{S_i}\left(W^{(0)},V_i^{(0)},S_i^{(0)}\right)$ are the same as in the semiclassical case.

In the system (32), (34)–(36) are not closed for the presence of the unknown quantities $U_{ij}$, $F_{ij}$, $C_{V_i}$, $C_W$ and $C_{S_i}$. We solve the closure problem with the approximation (29), assuming a collision dominated high field regime for the quantum effects.

In order to evaluate the unknown quantities present in the moment equations, the following formal lemmas are useful

Lemma 1.

$$\begin{array}{c}\hfill {\int }_{{\mathbb{R}}^3}{v}_{i_1}\cdots v_{i_k}{w}^{(1)}d\mathbf{v}={\left[i^k\frac{\partial }{\partial {\eta }_{i_1}}\cdots \frac{\partial }{\partial {\eta }_{i_k}}\mathcal{F}{w}^{(1)}(\eta )\right]}_{\eta =0}.\end{array}$$

(42)

Lemma 2.

$$\mathcal{F}{w}_{eq}^{(1)}(0)=\frac{\partial }{\partial {\eta }_{i_k}}\mathcal{F}{w}_{eq}^{(1)}(0)=\frac{\partial }{\partial {\eta }_i}\frac{\partial }{\partial {\eta }_j}\frac{\partial }{\partial {\eta }_k}\mathcal{F}{w}_{eq}^{(1)}(0)=0$$

(43)

$$\frac{\partial }{\partial {\eta }_i}\frac{\partial }{\partial {\eta }_j}\mathcal{F}{w}_{eq}^{(1)}(0)=n(\mathbf{x},t)\frac{{\hslash }^2\beta q}{12m^{*2}}\frac{{\partial }^{2}V}{\partial x_i\partial x_j}$$

(44)

$$\frac{\partial }{\partial {\eta }_i}\frac{\partial }{\partial {\eta }_j}{\mathrm{\Delta }}_{\eta }\mathcal{F}{w}_{eq}^{(1)}(0)=-n(\mathbf{x},t)\frac{{\hslash }^{2}q}{12m^{*3}}\left[\mathrm{\Delta }V{\delta }_{ij}+5\frac{{\partial }^{2}V}{\partial x_i\partial x_j}\right].$$

(45)

The proof follows by a simple computation. With the aid of these lemmas, we get the following closure relations:

Proposition 2.

$$\begin{array}{ccc}\hfill J_i& =& n{V}_i=n{V}_i^{(0)}+\mathcal{O}({\hslash }^4),\hfill \end{array}$$

(46)

$$\begin{array}{ccc}\hfill W& =& W^{(0)}-\frac{{\hslash }^2\beta q}{24m^{*}}\mathrm{\Delta }V+\mathcal{O}({\hslash }^4),\hfill \end{array}$$

(47)

$$\begin{array}{ccc}\hfill U_{ij}& =& \frac{2}{3}{\displaystyle \frac{W^{(0)}}{m^{*}}}{\delta }_{ij}-\frac{{\hslash }^2\beta q}{12{(m^{*})}^2}\frac{{\partial }^{2}V}{\partial x_i\partial x_j}+\mathcal{O}({\hslash }^4),\hfill \end{array}$$

(48)

$$\begin{array}{ccc}\hfill S_i& =& S_i^{(0)}-\frac{{\hslash }^2{\beta }^2{q}^2}{24m^{*2}\nu }\left(2\frac{{\partial }^{2}V}{\partial x_i\partial x_r}\frac{\partial V}{\partial x_r}+\frac{\partial V}{\partial x_i}\mathrm{\Delta }V\right)-\frac{{\hslash }^{2}q}{8m^{*2}\nu }\frac{\partial }{\partial x_i}\mathrm{\Delta }V+\mathcal{O}({\hslash }^4),\hfill \end{array}$$

(49)

$$\begin{array}{ccc}{F}_{ij}& =& {\displaystyle \frac{10}{9m^{*}}}{(W^{(0)})}^2-{\displaystyle \frac{{\hslash }^2\beta q^3}{3m^{*3}{\nu }^2}}{\displaystyle \frac{\partial V}{\partial x_{(i}}}{\displaystyle \frac{{\partial }^{2}V}{\partial x_{j)}\partial x_r}}\frac{\partial V}{\partial x_r}\hfill \\ & & -{\displaystyle \frac{{\hslash }^2{q}^2}{4m^{*3}{\nu }^2}}{\displaystyle \frac{{\partial }^{3}V}{\partial x_i\partial x_j\partial x_r}}{\displaystyle \frac{\partial V}{\partial x_r}}-{\displaystyle \frac{{\hslash }^2\beta q^3}{12m^{*3}{\nu }^2}}\left({\displaystyle \frac{\partial V}{\partial x_i}}{\displaystyle \frac{\partial V}{\partial x_j}}\mathrm{\Delta }V+{|\nabla V|}^2{\displaystyle \frac{{\partial }^{2}V}{\partial x_i\partial x_j}}\right)\hfill \\ & & -{\displaystyle \frac{{\hslash }^2{q}^2}{4m^{*3}{\nu }^2}}{\displaystyle \frac{\partial \mathrm{\Delta }V}{\partial x_{(i}}}{\displaystyle \frac{\partial V}{\partial x_{j)}}}-{\displaystyle \frac{{\hslash }^{2}q}{24m^{*2}}}\left(\mathrm{\Delta }V{\delta }_{ij}+5{\displaystyle \frac{{\partial }^{2}V}{\partial x_i\partial x_j}}\right)\hfill \\ & & -{\displaystyle \frac{{\hslash }^{2}q}{4m^{*2}\nu }}\left({\displaystyle \frac{\partial \mathrm{\Delta }V}{\partial x_{(i}}}{V}_{j)}+{\displaystyle \frac{{\partial }^{3}V}{\partial x_{(i}\partial x_j\partial x_{k)}}}{V}_k\right)+\mathcal{O}({\hslash }^4).\hfill \end{array}$$

(50)

In the previous relationships, parentheses indicate symmetrization, e.g.,

$$\begin{array}{ccc}\hfill A_{i(jk)}& =& \frac{1}{2}\left(A_{ijk}+A_{ikj}\right),\hfill \\ \hfill A_{(ijk)}& =& \frac{1}{3!}\left(A_{ijk}+A_{ikj}+A_{jik}+A_{jki}+A_{kij}+A_{kji}\right).\hfill \end{array}$$

For the proof, see [15].

Remark 3.

The zero order in ${\hslash }^2$ is the same as that obtained in [10,30].

Remark 4.

Since

$${\int }_{{\mathbb{R}}^3}{w}^{(1)}d\mathbf{v}={\int }_{{\mathbb{R}}^3}{w}_{eq}^{(1)}d\mathbf{v}=0,$$

we have

$${\int }_{{\mathbb{R}}^3}\mathcal{C}[w]d\mathbf{v}=0,$$

that is the assumption (25) is fulfilled.

Remark 5.

In the limit of high frequency $\nu \to \infty$ one has the simplified model

$$\begin{array}{ccc}\hfill J_i& =& n{V}_i=n{V}_i^{(0)}+\mathcal{O}({\hslash }^4),\hfill \end{array}$$

(51)

$$\begin{array}{ccc}\hfill W& =& W^{(0)}-\frac{{\hslash }^2\beta q}{24m^{*}}\mathrm{\Delta }V+\mathcal{O}({\hslash }^4)\hfill \end{array}$$

(52)

$$\begin{array}{ccc}\hfill U_{ij}& =& \frac{2}{3}\frac{W^{(0)}}{m^{*}}{\delta }_{ij}-\frac{{\hslash }^2\beta q}{12{(m^{*})}^2}\frac{{\partial }^{2}V}{\partial x_i\partial x_j}+\mathcal{O}({\hslash }^4)\hfill \end{array}$$

(53)

$$\begin{array}{ccc}\hfill S_i& =& S_i^{(0)}+\mathcal{O}({\hslash }^4),\hfill \end{array}$$

(54)

$$\begin{array}{ccc}\hfill F_{ij}& =& \frac{10}{9m^{*}}{(W^{(0)})}^2-\frac{{\hslash }^{2}q}{24m^{*2}}\left(\mathrm{\Delta }V{\delta }_{ij}+5\frac{{\partial }^{2}V}{\partial x_i\partial x_j}\right)\hfill \end{array}$$

(55)

From Equation (28), one sees that in the limit $\nu \to \infty$, $w^{(1)}$ reduces to the quantum correction of the equilibrium Wigner function $w_{eq}^{(1)}$. The resulting quantum corrections to the tensor $U_{ij}$ are the same as those obtained in [43] by using a shifted Wigner function, but with the semiclassical contribution which also contains a heat flux, not added ad hoc.

Remark 6.

Under the assumption (see [43]) that the equilibrium relation

$$\nabla V={(q\beta )}^{-1}\nabla \ln n+\mathcal{O}({\hslash }^2)$$

(56)

is valid, one can recover a formulation of the quantum corrections of density gradient type.

The density gradient version is more familiar because the quantum corrections take a form similar to the Bohm potential arising in the Madelung model of quantum fluids in the zero temperature limit. Moreover, some numerical experiments [43] lead to it being considered more robust than the original formulation in terms of the derivatives of the electric potential. In the limit $\nu \to \infty$, the closure relations in the density gradient form explicitly read as

$$\begin{array}{ccc}\hfill J_i& =& n{V}_i=n{V}_i^{(0)}+\mathcal{O}({\hslash }^4),\hfill \end{array}$$

(57)

$$\begin{array}{ccc}\hfill W& =& W^{(0)}-\frac{{\hslash }^2}{24m^{*}}\mathrm{\Delta }\ln n+\mathcal{O}({\hslash }^4),\hfill \end{array}$$

(58)

$$\begin{array}{ccc}\hfill U_{ij}& =& \frac{2}{3}\frac{W^{(0)}}{m^{*}}{\delta }_{ij}-\frac{{\hslash }^2}{12{(m^{*})}^2}\frac{{\partial }^2\ln n}{\partial x_i\partial x_j}+\mathcal{O}({\hslash }^4),\hfill \end{array}$$

(59)

$$\begin{array}{ccc}\hfill S_i& =& S_i^{(0)}+\mathcal{O}({\hslash }^4),\hfill \end{array}$$

(60)

$$\begin{array}{ccc}\hfill F_{ij}& =& \frac{10}{9m^{*}}n{(W^{(0)})}^2-\frac{{\hslash }^2}{24m^{*2}\beta }\left(\mathrm{\Delta }\ln n{\delta }_{ij}+5\frac{{\partial }^2\ln n}{\partial x_i\partial x_j}\right).\hfill \end{array}$$

(61)

4. Two-Dimensional Electron Gases: The Case of Quantum Confinement

In this section, we consider the case of an ensemble of electrons (the treatment of holes would be analogous) confined along one dimension, called quasi 2Delectron gas (2DEG) [16,18,19,20,44]. This situation arises when the length scale in one (the confined) space direction of the semiconductor device under study is of the order of de Broglie wavelength of electrons, while the nonconfined directions have a much bigger length scale. In other words, electrons are in a quantum regime in the confined direction and exhibit a classical behaviour in the nonconfined ones [16] as shown in Figure 1.

Let us suppose that electrons are quantized along one direction, which we choose as the z direction, and free to move in the orthogonal x-y plane. Let $\mathsf{\Omega }=[0,L_x]\times [0,L_y]\times [0,L_z]\subset {\mathbb{R}}^3$ be the spatial domain. Here, $L_x$, $L_y$, $L_z>0$ are fixed. More general cases with a variable confining direction can be easily incorporated.

Let us introduce the parameter σ defined as the ratio between the transversal (perpendicular to the x-y plane) typical length scale $L_t$ and the longitudinal typical length scale $L_l$

$$\sigma =\frac{L_t}{L_l}$$

and assume that $\sigma \ll 1$. Moreover, as is customary, let us assume the following ansatz about the electron wave function

$$\psi (\mathbf{k},\mathbf{r})=\psi (k_x,k_y,k_z,x,y,z)=\frac{1}{\sqrt{\mathcal{A}}}\phi (\sigma {\mathbf{r}}_{\left|\right|},z)e^{i{\mathbf{k}}_{\left|\right|}·\sigma {\mathbf{r}}_{\left|\right|}}\mathrm{in}\mathsf{\Omega },$$

(62)

with ${\mathbf{k}}_{\left|\right|}=(k_x,k_y)$ and ${\mathbf{r}}_{\left|\right|}=(x,y)$ denoting the longitudinal components of the wave-vector $\mathbf{k}$ and the position vector $\mathbf{r}$, respectively, and $\mathcal{A}$ being the area of the x-y cross-section. The previous wave function is inserted into the stationary Schrödinger equation in the effective mass approximation

$$\left[-\frac{{\hslash }^2}{2m^{*}}\mathrm{\Delta }+E_C(\mathbf{r})\right]\psi =E\psi ,$$

(63)

where $m^{*}$ is the electron effective mass, $E_C$ is the conduction band minimum, $E_C=-q(V_C+V)$, with $V_C(z)$ the confining potential and $V({\mathbf{r}}_{\left|\right|},z)$ the self-consistent electrostatic potential.

Under the assumption that the confining potential is so high that it gives rise to a barrier which can be considered infinite, it enters into the equations only through the boundary conditions

$$\psi ({\mathbf{r}}_{\left|\right|},z)=0atz=0\mathrm{and}z=L_z,$$

and $E_C=-qV$ is taken.

Therefore, introducing the slowly varying variable ${\tilde{\mathbf{r}}}_{\left|\right|}=\sigma {\mathbf{r}}_{\left|\right|}$ and inserting (62) into (63), one has

$$\begin{array}{c}\hfill -\frac{{\hslash }^2}{2m^{*}}{\sigma }^2\left[{\mathrm{\Delta }}_{{\mathbf{r}}_{\left|\right|}}\phi ({\tilde{\mathbf{r}}}_{\left|\right|},z)+2i{\mathbf{k}}_{\left|\right|}·{\nabla }_{{\mathbf{r}}_{\left|\right|}}\phi ({\tilde{\mathbf{r}}}_{\left|\right|},z)-{\mathbf{k}}_{\left|\right|}^2\phi ({\tilde{\mathbf{r}}}_{\left|\right|},z)\right]\\ \hfill -\frac{{\hslash }^2}{2m^{*}}\frac{\partial }{\partial z^2}\phi ({\tilde{\mathbf{r}}}_{\left|\right|},z)-qV({\tilde{\mathbf{r}}}_{\left|\right|},z)\phi ({\tilde{\mathbf{r}}}_{\left|\right|},z)=E\phi ({\tilde{\mathbf{r}}}_{\left|\right|},z).\end{array}$$

In the limit $\sigma \to 0^{+}$, one formally gets that the envelope function φ must solve the reduced Schrödinger equation

$$\left[-\frac{{\hslash }^2}{2m^{*}}\frac{d^2}{d{z}^2}-qV({\mathbf{r}}_{\left|\right|},z)\right]\phi (\sigma {\mathbf{r}}_{\left|\right|},z)=\epsilon \phi (\sigma {\mathbf{r}}_{\left|\right|},z),0\le z\le L_z,$$

(64)

with boundary conditions

$$\phi =0\mathrm{at}z=0\mathrm{and}z=L_z.$$

(65)

The solution of (64) and (65) depends only parametrically on ${\mathbf{r}}_{\left|\right|}$ (and in a slow way), and more in general also on time t if a non-steady-state solution is considered. In fact, electrons as waves achieve equilibrium along the confined direction in a time which is much shorter than the typical transport time, so that one can adopt a quasi-static description along the z-direction. The couple (64) and (65) constitutes a one dimensional Sturm–Liouville problem, which admits a countable set of eigen-pairs (subbands) $({\phi }_{\nu }({\mathbf{r}}_{\left|\right|},z),{\epsilon }_{\nu }({\mathbf{r}}_{\left|\right|})),\nu =1,2...$, whose eigenvalues are simple and do not cross. The potential V is obtained from the Poisson equation

$$\nabla ·({ϵ}_s\nabla V)=-q(N_d(\mathbf{r})-n),$$

(66)

where the electron density n is given by

$$n(\mathbf{r})=\sum _{\nu =1}^{+\infty }{\rho }_{\nu }({\mathbf{r}}_{\left|\right|}){\left|{\phi }_{\nu }({\mathbf{r}}_{\left|\right|},z,)\right|}^2,$$

with ${\rho }_{\nu }$ the (areal) density of electrons of the νth subband and for simplicity an n–doped material has been considered. Of course, the Schrödinger and Poisson equations are coupled and must be solved simultaneously.

For devices with longitudinal characteristic length of a few tens of nanometers, the transport of electrons in the longitudinal direction is semi-classical within a good approximation. Therefore, electrons in each subband can be considered as different populations, the state of each of them being described by a semiclassical distribution function $f_{\nu }(x,y,k_x,k_y,t)$, and electron transport along the longitudinal direction being determined by adding to the Schrödinger–Poisson system the following system of coupled Boltzmann equations

$$\frac{\partial f_{\nu }}{\partial t}+\frac{1}{\hslash }{\nabla }_{{\mathbf{k}}_{\left|\right|}}{E}_{\nu }·{\nabla }_{{\mathbf{r}}_{\left|\right|}}{f}_{\nu }-\frac{q}{\hslash }{\mathbf{E}}_{\nu }^{eff}·{\nabla }_{{\mathbf{k}}_{\left|\right|}}{f}_{\nu }=\sum _{\mu =1}^{\infty }{C}_{\nu ,\mu }[f_{\nu },f_{\mu }],\nu =1,2,\cdots$$

(67)

where ${\displaystyle {\mathbf{E}}_{\nu }^{eff}=\frac{1}{q}{\nabla }_{{\mathbf{r}}_{\left|\right|}}{\epsilon }_{\nu }({\mathbf{r}}_{\left|\right|})}$, and the electron energy in each subband $E_{\nu }$ is the sum of a transversal contribution ${\epsilon }_{\nu }$ and a longitudinal one

$$E_{\nu }({\mathbf{r}}_{\left|\right|},{\mathbf{k}}_{\left|\right|}):={\epsilon }_{\nu }({\mathbf{r}}_{\left|\right|})+{\epsilon }_{\left|\right|}({\mathbf{k}}_{\left|\right|})={\epsilon }_{\nu }({\mathbf{r}}_{\left|\right|})+\frac{1}{2\alpha }\left(\sqrt{1+4\alpha \frac{{\hslash }^2}{2m^{*}}\left(k_x^2+k_y^2\right)}-1\right),$$

where α is the non-parabolicity parameter. Consequently, the longitudinal electron velocity is

$${\mathbf{v}}_{\left|\right|}=\frac{1}{\hslash }{\nabla }_{{\mathbf{k}}_{\left|\right|}}{\epsilon }_{\left|\right|}=\frac{\hslash {\mathbf{k}}_{\left|\right|}}{m^{*}\left(1+2\alpha {\epsilon }_{\left|\right|}\right)}.$$

(68)

A formal justification of Equation (67) can be found in [45]. Regarding the collision operator, in the non-degenerate approximation, each contribution has the general form [18,19,20]

$$C_{\nu ,\mu }[f_{\nu },f_{\mu }]=\frac{1}{{(2\pi )}^2}{\int }_{B_2}\left[w_{\mu \nu }({\mathbf{k}}_{\left|\right|}^{\prime },{\mathbf{k}}_{\left|\right|})f_{\mu }^{\prime }-w_{\nu \mu }({\mathbf{k}}_{\left|\right|},{\mathbf{k}}_{\left|\right|}^{\prime })f_{\nu }\right]d^2{\mathbf{k}}_{\left|\right|}^{\prime },$$

where $B_2$ is the two-dimensional Brillouin zone. When $\mu =\nu$, the scattering is intra-subband; when $\mu \ne \nu$, it is the inter-subband. We recall that inter-subband transitions may also happen in the absence of phonons, by quantum mechanical tunneling induced by the presence of a static field as has been shown for example in [46]. The term $w_{\mu \nu }({\mathbf{k}}_{\left|\right|}^{\prime },{\mathbf{k}}_{\left|\right|})$ is the transition rate from the longitudinal state with wave-vector ${\mathbf{k}}_{\left|\right|}^{\prime }$, belonging to the μth subband, to the longitudinal state with wave-vector ${\mathbf{k}}_{\left|\right|}$, belonging to the νth subband, and $f_{\mu }^{\prime }\equiv f_{\mu }({\mathbf{r}}_{\left|\right|},{\mathbf{k}}_{\left|\right|}^{\prime },t)$. In Si, the relevant 2D scattering mechanisms are the acoustic phonon scattering, and the non-polar optical phonon scattering. Their expressions are the 2D version of those shown in Section 2 and can be found in [18,19,20,47]. Attempts to directly solve the Schrödinger–Poisson–Boltzmann system can be found, for example, in [16,17,48] where numerical schemes based on finite differences have been employed. However, the direct numerical approach is a daunting computational task and requires very long computing times. This has again prompted the development of simpler macroscopic models for CAD purposes. These models can be obtained as moment equations of the Boltzmann transport equations under suitable closure relations [18,19,20]. The moment of the ν-th subband distribution with respect to a weight function $a({\mathbf{k}}_{\left|\right|})$ reads

$$M_{a,\nu }={\int }_{B_2}a({\mathbf{k}}_{\left|\right|})f_{\nu }({\mathbf{r}}_{\left|\right|},{\mathbf{k}}_{\left|\right|},t)d^2{\mathbf{k}}_{\left|\right|}.$$

In particular, analogously to the previous cases, we take as basic moments

$$\mathrm{the}\mathrm{areal}\mathrm{density}{\rho }_{\nu }={\int }_{B_2}{f}_{\nu }({\mathbf{r}}_{\left|\right|},{\mathbf{k}}_{\left|\right|},t)d^2{\mathbf{k}}_{\left|\right|},$$

(69)

$$\mathrm{the}\mathrm{longitudinal}\mathrm{mean}\mathrm{velocity}{\mathbf{V}}_{\nu }=\frac{1}{{\rho }_{\nu }}{\int }_{B_2}{\mathbf{v}}_{\left|\right|}{f}_{\nu }({\mathbf{r}}_{\left|\right|},{\mathbf{k}}_{\left|\right|},t)d^2{\mathbf{k}}_{\left|\right|},$$

(70)

$$\mathrm{the}\mathrm{longitudinal}\mathrm{mean}\mathrm{energy}{W}_{\nu }=\frac{1}{{\rho }_{\nu }}{\int }_{B_2}{\epsilon }_{\left|\right|}{f}_{\nu }({\mathbf{r}}_{\left|\right|},{\mathbf{k}}_{\left|\right|},t)d^2{\mathbf{k}}_{\left|\right|},$$

(71)

$$\mathrm{the}\mathrm{longitudinal}\mathrm{mean} \text{energy-flux}{\mathbf{S}}_{\nu }=\frac{1}{{\rho }_{\nu }}{\int }_{B_2}{\epsilon }_{\left|\right|}{\mathbf{v}}_{\left|\right|}{f}_{\nu }({\mathbf{r}}_{\left|\right|},{\mathbf{k}}_{\left|\right|},t)d^2{\mathbf{k}}_{\left|\right|}.$$

(72)

The corresponding moment system reads

$$\begin{array}{ccc}\hfill \frac{\partial {\rho }_{\nu }}{\partial t}+{\nabla }_{{\mathbf{r}}_{\left|\right|}}·({\rho }_{\nu }{\mathbf{V}}_{\nu })& =& {\rho }_{\nu }\sum _{\mu }{C}_{\rho ,\nu \mu },\hfill \end{array}$$

(73)

$$\begin{array}{ccc}\hfill \frac{\partial ({\rho }_{\nu }{\mathbf{V}}_{\nu })}{\partial t}+{\nabla }_{{\mathbf{r}}_{\left|\right|}}·({\rho }_{\nu }{\mathbf{F}}_{\nu }^{(0)})+{\rho }_{\nu }{\mathbf{G}}_{\nu }^{(0)}{\nabla }_{{\mathbf{r}}_{\left|\right|}}{\epsilon }_{\nu }& =& {\rho }_{\nu }\sum _{\mu }{C}_{\mathbf{V},\nu \mu },\hfill \end{array}$$

(74)

$$\begin{array}{ccc}\hfill \frac{\partial ({\rho }_{\nu }{W}_{\nu })}{\partial t}+{\nabla }_{{\mathbf{r}}_{\left|\right|}}·({\rho }_{\nu }{\mathbf{S}}_{\nu })+{\rho }_{\nu }{\nabla }_{{\mathbf{r}}_{\left|\right|}}{\epsilon }_{\nu }·{\mathbf{V}}_{\nu }& =& {\rho }_{\nu }\sum _{\mu }{C}_{W,\nu \mu },\hfill \end{array}$$

(75)

$$\begin{array}{ccc}\hfill \frac{\partial ({\rho }_{\nu }{\mathbf{S}}_{\nu })}{\partial t}+{\nabla }_{{\mathbf{r}}_{\left|\right|}}·({\rho }_{\nu }{\mathbf{F}}_{\nu }^{(1)})+{\rho }_{\nu }{\mathbf{G}}_{\nu }^{(1)}{\nabla }_{{\mathbf{r}}_{\left|\right|}}{\epsilon }_{\nu }& =& {\rho }_{\nu }\sum _{\mu }{C}_{\mathbf{S},\nu \mu },\hfill \end{array}$$

(76)

where

$$\begin{array}{ccc}\hfill \left(\begin{array}{c}{\mathbf{F}}_{\nu }^{(0)}\\ {\mathbf{F}}_{\nu }^{(1)}\end{array}\right)& =& \frac{1}{{\rho }_{\nu }}{\int }_{B_2}\left(\begin{array}{c}1\\ {\epsilon }_{\left|\right|}\end{array}\right){\mathbf{v}}_{\left|\right|}\otimes {\mathbf{v}}_{\left|\right|}{f}_{\nu }({\mathbf{r}}_{\left|\right|},{\mathbf{k}}_{\left|\right|},t)d^2{\mathbf{k}}_{\left|\right|},\hfill \end{array}$$

(77)

$$\begin{array}{ccc}\hfill \left(\begin{array}{c}{\mathbf{G}}_{\nu }^{(0)}\\ {\mathbf{G}}_{\nu }^{(1)}\end{array}\right)& =& \frac{1}{{\rho }_{\nu }}{\int }_{B_2}\left(\begin{array}{c}{\displaystyle \frac{1}{\hslash }{\nabla }_{{\mathbf{k}}_{\left|\right|}}{\mathbf{v}}_{\left|\right|}}\\ {\displaystyle \frac{1}{\hslash }{\nabla }_{{\mathbf{k}}_{\left|\right|}}\left({\epsilon }_{\left|\right|}{\mathbf{v}}_{\left|\right|}\right)}\end{array}\right)f_{\nu }({\mathbf{r}}_{\left|\right|},{\mathbf{k}}_{\left|\right|},t)d^2{\mathbf{k}}_{\left|\right|},\hfill \end{array}$$

(78)

$$\begin{array}{ccc}\hfill \left(\begin{array}{c}{C}_{\rho ,\nu \mu }\\ C_{W,\nu \mu }\end{array}\right)& =& \frac{1}{{\rho }_{\nu }}{\int }_{B_2}\left(\begin{array}{c}1\\ {\epsilon }_{\left|\right|}\end{array}\right)\left[w_{\mu \nu }({\mathbf{k}}_{\left|\right|}^{\prime },{\mathbf{k}}_{\left|\right|})f_{\mu }^{\prime }-w_{\nu \mu }({\mathbf{k}}_{\left|\right|},{\mathbf{k}}_{\left|\right|}^{\prime })f_{\nu }\right]d^2{\mathbf{k}}_{\left|\right|}^{\prime }{d}^2{\mathbf{k}}_{\left|\right|},\hfill \end{array}$$

(79)

$$\begin{array}{ccc}\hfill \left(\begin{array}{c}{C}_{\mathbf{V},\nu \mu }\\ C_{\mathbf{S},\nu ,\mu }\end{array}\right)& =& \frac{1}{{\rho }_{\nu }}{\int }_{B_2}\left(\begin{array}{c}{\mathbf{v}}_{\left|\right|}\\ {\epsilon }_{\left|\right|}{\mathbf{v}}_{\left|\right|}\end{array}\right)\left[w_{\mu \nu }({\mathbf{k}}_{\left|\right|}^{\prime },{\mathbf{k}}_{\left|\right|})f_{\mu }^{\prime }-w_{\nu \mu }({\mathbf{k}}_{\left|\right|},{\mathbf{k}}_{\left|\right|}^{\prime })f_{\nu }\right]d^2{\mathbf{k}}_{\left|\right|}^{\prime }{d}^2{\mathbf{k}}_{\left|\right|}.\hfill \end{array}$$

(80)

Also in this case, the above written fluxes and production terms are extra-variables for which closure relations are needed. In order to employ MEP, a suitable expression of the entropy for the system under consideration has to be found. Neglecting the mutual interaction of electrons in different subbands, considering the phonon gas as a thermal bath and assuming the electron gas is sufficiently dilute, in each subband the expression of the entropy is the semiclassical limit of that arising from the Fermi statistics. Eventually, for obtaining the total entropy, one has to consider that electrons have a quantum behaviour along the z-direction, therefore we define the entropy density of the system as

$$\mathcal{S}[f_1,f_2,\cdots ]=-k_B\sum _{\nu =1}^{+\infty }\underset{\mathrm{quantum}\mathrm{effect}}{\underbrace{|{\phi }_{\nu }({\mathbf{r}}_{\left|\right|},z,t){|}^2}}\underset{\mathrm{semiclassical}\mathrm{contribution}}{\underbrace{{\int }_{B_2}\left(f_{\nu }\log (2{\pi }^2{f}_{\nu })-f_{\nu }\right)d^2{\mathbf{k}}_{\left|\right|}}}.$$

(81)

Remark 7.

The proposed expression of the entropy has been introduced for the first time in [20,47]. It combines quantum effects and semiclassical transport along the longitudinal direction, weighting the contribution of each subband with the square modulus of the ${\phi }_{\nu }(z,t)$’s arising from the Schrödinger–Poisson block.

According to MEP, the $f_{\nu }$’s are estimated by means of the distributions $f_{\nu }^{MEP}$’s that $\forall {\mathbf{r}}_{\left|\right|}\in [0,L_x]\times [0,L_y]$, $t>0,$ solve the problem:

$$\begin{array}{c}{\displaystyle \underset{\begin{array}{c}{f}_1({\mathbf{r}}_{\left|\right|},·,t)\in {\mathcal{F}}_{2D}\\ f_2({\mathbf{r}}_{\left|\right|},·,t)\in {\mathcal{F}}_{2D}\\ \cdots \end{array}}{\max }\mathcal{S}[f_1,f_2,\cdots ]}\hfill \\ \mathrm{under}\mathrm{the}\mathrm{constraints}{M}_{a_A,\nu }={\int }_{B_2}{a}_A({\mathbf{k}}_{\left|\right|})f_{\nu }({\mathbf{r}}_{\left|\right|},{\mathbf{k}}_{\left|\right|},t)d^2{\mathbf{k}}_{\left|\right|},\hfill \end{array}$$

where the $M_{a_A,\nu }$’s are the basic moments (69)–(72) and ${\mathcal{F}}_{2D}$ is the space of the summable function with respect to ${\mathbf{k}}_{\left|\right|}$ such that the moments appearing in the system (73)–(76) exist.

With the above choice of the functions ${\displaystyle a_A({\mathbf{k}}_{\left|\right|})=(1,{\mathbf{v}}_{\left|\right|},{\epsilon }_{\left|\right|},{\epsilon }_{\left|\right|}{\mathbf{v}}_{\left|\right|})}$, the resulting maximum entropy distribution functions read (the factor $k_B$ has been included into the multipliers)

$$f_{\nu }^{MEP}=\exp \left[-\left({\lambda }_{\nu }+{\lambda }_{\mathbf{V},\nu }·{\mathbf{v}}_{\left|\right|}+\left({\lambda }_{W,\nu }+{\lambda }_{\mathbf{S},\nu }·{\mathbf{v}}_{\left|\right|}\right){\epsilon }_{\left|\right|}\right)\right].$$

In order to complete the procedure, similarly to the previous case, one has to insert the $f_{\nu }^{MEP}$’s into the constraints (69)–(72) and express the Lagrange multipliers as functions of the basic moments ${\rho }_{\nu }$, ${\mathbf{V}}_{\nu }$, $W_{\nu }$, ${\mathbf{S}}_{\nu }$. In general, in the case considered in this section as well as in that considered in Section 2, it is not possible to assure that such a procedure can be accomplished. The solvability of the MEP problem depends on the band structure and on the choice of the moments, as well known already in gas-dynamics [14,49]. For example, if the parabolic case (where the energy is a quadratic function of the wave vector) is considered, the same drawbacks of classical gas-dynamics arise regarding the lack of integrability of the MEP distribution function. In the case of the Kane dispersion relation, the solvability of the MEP problem is guaranteed by the following property, proven in [49].

Proposition 3.

When the Kane model for the energy band is employed, the corresponding maximum entropy models are symmetric hyperbolic systems with convex domains of definition and the equilibria are interior points, guaranteeing the validity of expansions around equilibrium states.

The proof has been given for a 3D electron gas, but it can be extended to a 2DEG in a straightforward way.

For the 2DEG case, the final expressions of the extra-fluxes and production terms, obtained by using linearized MEP distribution functions, can be found in [18,47] for the parabolic approximation and in [19,20] for the non-parabolic one. The resulting moment system has been used to simulate double-gate MOSFETs in [18,19].

5. Some Numerical Simulations

In this last section, we present some of the numerical results obtained by implementing the models sketched in the previous sections. The information concerning the numerical schemes will be skipped. The interested readers are referred to the indicated papers. For optimization of electron devices using the models presented above, see [50] and references therein.

5.1. p-n Junction

The full model for holes and electrons, presented in Section 2, has been used for the numerical simulation of a Si p-n junction diode with a doping profile given by [22]

$$N_d(x)-N_a(x)=\left\{\begin{array}{c}-c_0\mathrm{if}0<x<L_0\hfill \\ c_0\mathrm{if}{L}_0<x<L,\hfill \end{array}\right.$$

with an abrupt junction at $L_0$. The physical parameters of the device are reported in

Table 1. Physical parameters of the p-n junction.

Parameter	Physical Meaning	Numerical Value
L	device length	${10}^{-3}$cm
${\mathrm{L}}_0$	length of p-region	0.7 × ${10}^{-3}$cm
${\mathrm{c}}_0$	doping concentration	${10}^{15}$ ${\mathrm{cm}}^{-3}$
${\mathrm{n}}_i$	intrinsic concentration	1.075 × ${10}^{10}$ ${\mathrm{cm}}^{-3}$

. The numerical solution has been obtained by resorting to an extension [51] of the traditional central differencing scheme for one-dimensional balance laws with (possibly stiff) source terms, which has been developed on the basis of the Nessyhau and Tadmor scheme [52] for homogeneous hyperbolic systems.

The bias voltage at the contacts is the superposition of the thermal equilibrium boundary potential (the built in potential $V_{bi}$) and the applied potential $V_a$. In Figure 2 and Figure 3, the stationary carrier densities and the electrostatic potential are shown for the following applied voltages: $V_a=-0.2,0,0.75,1$ V, the sign + referring to direct polarization.

5.2. Simulation of a MESFET

In this subsection, we show the simulation of a bi-dimensional Metal Semiconductor Field Effect Transistor (MESFET) made with the MEP model. The numerical method is based on the discretization proposed in [53].

The shape of the device is pictured in Figure 4. The device has a 0.4 $\mathsf{\mu }$m channel. The source and drain lengths are 0.1 $\mathsf{\mu }$m and the contact at the gate is 0.2 $\mathsf{\mu }$m long. The distance between the gate and the other two contacts is 0.1 $\mathsf{\mu }$m. The lateral subdiffusion of the source and the drain region is about 0.05 $\mathsf{\mu }$m. The same doping concentration $C(x)$ as in [53] is considered

$$C(x)=\left\{\begin{array}{cc}{n}_{+}=3\times {10}^{17}{\mathrm{cm}}^{-3}\hfill & \mathrm{in}\mathrm{the}{n}^{+}\mathrm{region}\hfill \\ n_{-}={10}^{17}{\mathrm{cm}}^{-3}\hfill & \mathrm{in}\mathrm{the}n\mathrm{region}\hfill \end{array}\right.$$

with abrupt junctions. We take a reference frame with axes parallel to the edges of the device. The numerical domain representing the MESFET is

$$\mathsf{\Omega }=[0,0.6]\times [0,0.2]$$

where the unit length is the micron.

The regions of high doping $n^{+}$ are the subset

$$[0,0.1]\times [0.15,0.2]\cup [0.5,0.6]\times [0.15,0.2].$$

${\Gamma }_D$ represents the source, gate and drain

$${\Gamma }_D=\left\{(x,y):y=0.2,0\le x\le 0.1,0.2\le x\le 0.4,0.5\le x\le 0.6\right\}.$$

The remaining part of $\partial \mathsf{\Omega }$ is ${\Gamma }_N$. The boundary conditions are assigned as follows. We have ohmic contacts at the source and drain:

$$\begin{array}{c}\hfill n=n_{+},W=\frac{3}{2}{k}_B{T}_L,V=\left\{\begin{array}{cc}{\Phi }_b& \mathrm{at}\mathrm{drain}\hfill \\ 0& \mathrm{at}\mathrm{source}.\hfill \end{array}\right.\end{array}$$

(82)

At the gate we have a Schottky contact

$$n=n_g,W=\frac{3}{2}{k}_B{T}_L,V={\Phi }_g.$$

(83)

Indeed, the potential at the contacts should include the built-in potential and the density at the gate should be related to the potential. Here, we do not enter into the details of the modeling of the Schottky contacts and, by using the invariance of the electric field with respect to changes of the potential for additive constants, we choose as in [54]

$$n_g=3.9\times {10}^5{\mathrm{cm}}^{-3},{\Phi }_g=-0.8V$$

In the remaining part ${\Gamma }_N$ of the boundary, we have

$$\begin{array}{c}\hfill \mathit{\nu }·\nabla n=0,\mathit{\nu }·\nabla W=0,\mathit{\nu }·\nabla V=0.\end{array}$$

(84)

Here, ∇ is the bi-dimensional gradient operator while ν is the unit outward normal vector to $\partial \mathsf{\Omega }$ in the considered points.

A uniform mesh has been used with 33 × 97 grid points. The stationary solution is reached after about 5 picoseconds. The numerical results are shown in Figure 5, Figure 6, Figure 7 and Figure 8. The typical depletion region close to the gate with the presence of steep gradients is numerically well described.

5.3. Simulation of a MOSFET

The Metal Oxide Semiconductor Field Effect Transistor (MOSFET) is a largely used device in modern electronics. The same mathematical model and numerical scheme as for the MESFET is adopted; the shape of the device is pictured in Figure 9. The device has a 0.2 $\mathsf{\mu }$m channel. The source and drain lengths are 0.1 $\mathsf{\mu }$m and the contact at the gate is 0.15 $\mathsf{\mu }$m long. The distance between the gate and the other two contacts is 0.025 $\mathsf{\mu }$m. The lateral subdiffusion of the source and the drain region is about 0.05 $\mathsf{\mu }$m. The gate oxide is 0.15 $\mathsf{\mu }$m long and 20 nm thick. Of course, smaller devices can be considered without any substantial modification to the used approach.

The doping concentration is

$$C(x)=\left\{\begin{array}{c}{n}_{+}={10}^{17}{\mathrm{cm}}^{-3}\mathrm{in}\mathrm{the}{n}^{+}\mathrm{regions}\hfill \end{array}{p}_{-}={10}^{14}{\mathrm{cm}}^{-3}\mathrm{in}\mathrm{the}p\mathrm{region}\right.$$

with abrupt junctions.

At variance with MESFET, there are different built-in potentials which we explicitly take into account by using the simple model

$${\Phi }_b^D={\Phi }_b^S=\frac{k_B{T}_L}{q}\ln \frac{n_{+}}{n_i}$$

at drain and source,

$${\Phi }_b^B=-\frac{k_B{T}_L}{q}\ln \frac{p_{-}}{n_i}$$

at bulk contact. Here, $n_i$ is the intrinsic electron concentration (${10}^{10}$ ${\mathrm{cm}}^{-3}$). The reference axes are chosen parallel to the edges of the device. The silicon part of the MOSFET is represented by the numerical domain

$$[0,0.4]\times [0,0.4]$$

and at the top of the silicon part the silicon oxide domain is

$$[0.125,0.275]\times [0.4,0.42]$$

where the unit length is the micron.

The regions of high-doping $n^{+}$ are the subset

$$[0,0.1]\times [0.35,0.4]\cup [0.3,0.4]\times [0.35,0.4].$$

A uniform mesh of 64 × 64 grid points has been used in the silicon part. The Poisson equation is solved on the entire (silicon and oxide) domain. Of course, in the oxide, the Poisson equation becomes the Laplace one.

We have assumed ohmic contacts at the source, drain, gate and bulk contacts, to be homogeneous Neumann conditions on the remaining part of the boundary. The surface charge at the oxide interface is neglected and the continuity of the electric potential and electric field is imposed. The values of density and energy at the interface are obtained by the interior grid points with a linear interpolation in the direction orthogonal to the boundary.

In order to reach the desired bias, we have needed to resort to a continuation method on applied potential. First, we iterate with respect to the difference of the built-in potential between drain and bulk contacts, keeping at zero $V_{DS}$. Then, we iterate with respect to the drain–gate potential and finally we increase the drain–source potential.

All the main features of the electron dynamics are well described, in particular the charge accumulation beside the oxide, and the pronounced depletion at the drain contact due to the strong electric field.

Again, the density current presents a singularity at the first edge of the drain and therefore we evaluated the current by considering, as for MESFET, the regularization from the interior.

As for the MESFET device, current gain and miniaturization are crucial tasks to take into account when designing this kind of device. However, miniaturization plays a crucial role for MOSFETs for various reasons; first, a shorter MOSFET helps the current flow, since a small device decreases the resistance. Second, a small MOSFET helps in reducing gate capacity and, in general, to have a faster switch on/off operation. The numerical results are presented in Figure 10, Figure 11 and Figure 12.

5.4. Simulation of a DG-MOSFET

The 2DEG model, illustrated in Section 4, has been employed to simulate the nanoscale silicon DG-MOSFET represented in Figure 13, including the non-parabolicity effects [18,19]. The length of the diode is $L_x$ = 40 nm, the width of the silicon layer is $L_z=8$ nm and the thickness of each oxide layer is $t_{ox}$ = 1 nm. The highly doped $n^{+}$ regions are 10 nm long. The gate contacts have the same length as the n region and are above it. The device is supposed to be infinite in the y direction.

The doping in the $n^{+}$ regions is $N_d(x)$ = $N_d^{+}$ = ${10}^{20}$ ${\mathrm{cm}}^{-3}$ and in the n region is $N_d(x)$ = $N_d^{-}$ = ${10}^{15}$ ${\mathrm{cm}}^{-3}$, with a regularization at the two junctions by a hyperbolic tangent profile

$$N_d(x)=N_d^{+}-d\left(\tanh \frac{x-x_1}{s}-\tanh \frac{x-x_2}{s}\right),$$

where $s=0.2$ nm, ${\displaystyle d=\frac{N_d^{+}}{2}\left(1-\frac{N_d^{-}}{N_d^{+}}\right)}$, $x_1=10$ nm and $x_2=30$ nm.

Due to the symmetries and dimensions of the device, the transport is, within a good approximation, one dimensional and along the longitudinal direction with respect to the two oxide layers, while the electrons are quantized in the transversal direction. The oxide gives rise, with a good approximation, to an infinitely deep potential barrier; in fact, realistic values of the potential barrier are more than 3 eV high and it is very unlikely to find electrons with such an energy in the device under consideration. Six equivalent valleys are considered with a single effective mass $m^{*}=0.32m_e$. For details about the appropriate boundary and initial conditions, as well as the numerical method, the interested reader is referred to [18,19]. The numerical experiments indicate that it is sufficient to take into account only the first three subbands.

As the first case, a symmetric situation with $V_D=0.5$ V and $V_{gl}=V_{gu}=-3$ V is considered, where $V_D$ is the voltage applied at the drain with respect to that at the source, and $V_{gl}$ and $V_{gu}$ are the voltages applied at the lower and the upper gate, respectively. In Figure 14 and Figure 15, the steady state density and the potential are plotted. The solution does not present any spurious oscillation or boundary layer and reflects the symmetry of the problem.

As the second case, $V_D=0.5$ V, $V_{gl}=-3$ V and $V_{gu}=3$ V are taken. In Figure 16 and Figure 17, the density and the potential are respectively plotted, while in Figure 18b the first three subband bottoms are shown. One can note the depletion region beneath the upper gate.

In Figure 18, the first three subband bottoms are shown. There is a good qualitative agreement with the other numerical simulations known in the literature [17].

From Figure 19, a very accurate current conservation is evident, proving the robustness of the numerical method. In the second case, the current is reduced by one half, due to the gate voltage in agreement with the behaviour of the density. The areal density, the average velocity and the energy measured from the subband bottom, and the current in the first three subbands are shown in Figure 20, Figure 21 and Figure 22. The areal density is not symmetrical between the source and drain within each subband but the total areal density is so.

The drift (mean) velocity has been evaluated according to the formula

$$\mathbf{V}=\frac{{\sum }_{\nu }{\rho }^{\nu }{\mathbf{V}}^{\nu }}{{\sum }_{\mu }{\rho }^{\mu }}.$$

Similarly, the global longitudinal mean energy has been evaluated taking as reference value the bottom of the first subband according to the formula

$$W=\frac{{\sum }_{\nu }{\rho }^{\nu }\left(W^{\nu }+{\epsilon }_{\nu }-{\epsilon }_1\right)}{{\sum }_{\mu }{\rho }^{\mu }}.$$

The maximum drift velocity in the channel is one and half times the saturation velocity when $V_{gl}=V_{gu}=-3$ V, while it is about two times the saturation velocity when $V_{gl}=-3$ V and $V_{gu}=3$ V. Moreover, in the first case, the velocity in the first subband is lower than that in the first subband in the second case. Instead, the velocity in the second and the third subbands is higher in the first case, but with a resulting total longitudinal current which is lower. The energy has an evidently different value between the source and drain as happens in the semiclassical case.

Eventually, the characteristic curves are shown in Figure 23 by fixing $V_{gl}=-3$ V and varying $V_{gu}$ from −3 V to 3 V. With increasing $V_{gu}$, the average longitudinal current increases as a consequence of the controlling effect of the gate voltage on the electric characteristics of the device.

6. Conclusions

A review on the exploitation of the Maximum Entropy Principle in the formulation of macroscopic models able to describe the charge and heat transport in semiconductor devices has been presented. The models are obtained starting from the Boltzmann transport equations for the charge and the phonon distribution functions, by taking—as macroscopic variables—suitable moments of the distributions and resorting to MEP in order to close the evolution equations for the chosen moments. Firstly, semiclassical models have been reviewed and eventually an approach to also take into account quantum effects has been shown. The latter consists of using the moment system arising from the Wigner equation.

In the cases where quantum effects are relevant only along one direction, another strategy has been introduced. A quasi-static description has been adopted, based on a coupled Schrödinger–Poisson system, along the confining direction, leading to a subband decomposition of the electron energy levels, while the transport along the longitudinal directions is described by a semiclassical Boltzmann equation for each subband.

Several physical situations, which have been investigated in the last two decades, have been reported in the paper. For the sake of conciseness, some further applications have been omitted. A MEP based model for 2D–3D electron gases, of particular interest in the simulation of nanoscale MOSFET, can be found in [55]. A hydrodynamical model for charge transport in graphene has been formulated in [56,57].

Apparently, MEP is revealed to be a sound method for devising, in a systematic way, macroscopic models, e.g., energy–transport and hydrodynamical ones, which are less complex to tackle numerically but still retain a good accuracy with respect to the results based on the direct solutions of the transport equations.