Next Article in Journal
Chemical Reactions Using a Non-Equilibrium Wigner Function Approach
Next Article in Special Issue
Exploitation of the Maximum Entropy Principle in Mathematical Modeling of Charge Transport in Semiconductors
Previous Article in Journal
From Tools in Symplectic and Poisson Geometry to J.-M. Souriau’s Theories of Statistical Mechanics and Thermodynamics
Previous Article in Special Issue
Maximum Entropy Closure of Balance Equations for Miniband Semiconductor Superlattices
 
 
Font Type:
Arial Georgia Verdana
Font Size:
Aa Aa Aa
Line Spacing:
Column Width:
Background:
Article

A Hydrodynamic Model for Silicon Nanowires Based on the Maximum Entropy Principle

by
Orazio Muscato
*,† and
Tina Castiglione
Dipartimento di Matematica e Informatica, Università di Catania, Catania 95125, Italy
*
Author to whom correspondence should be addressed.
These authors contributed equally to this work.
Entropy 2016, 18(10), 368; https://doi.org/10.3390/e18100368
Submission received: 14 July 2016 / Revised: 6 September 2016 / Accepted: 30 September 2016 / Published: 19 October 2016
(This article belongs to the Special Issue Maximum Entropy Principle and Semiconductors)

Abstract

:
Silicon nanowires (SiNW) are quasi-one-dimensional structures in which the electrons are spatially confined in two directions, and they are free to move along the axis of the wire. The spatial confinement is governed by the Schrödinger–Poisson system, which must be coupled to the transport in the free motion direction. For devices with the characteristic length of a few tens of nanometers, the transport of the electrons along the axis of the wire can be considered semiclassical, and it can be dealt with by the multi-sub-band Boltzmann transport equations (MBTE). By taking the moments of the MBTE, a hydrodynamic model has been formulated, where explicit closure relations for the fluxes and production terms (i.e., the moments on the collisional operator) are obtained by means of the maximum entropy principle of extended thermodynamics, including the scattering of electrons with phonons, impurities and surface roughness scattering. Numerical results are shown for a SiNW transistor.

1. Introduction

As integrating smaller-sized high performance electronic devices has become a key point over the past few decades, innovative transistor designs have emerged. Chip-to-wafer density has to remain high to be competitive, so size reduction is still a burning issue, and non-planar devices are part of these brand new products. As a consequence, multi-gate field-effect transistors (FET) have been developed. In such a new gate architecture, silicon nanowires (SiNWs) FETs seem to be one of the most attractive choices. They enable avoiding problems due to the short channel effect, such as threshold voltage or drain-induced barrier lowering. More precisely, surrounding-gate transistors where the gate encircles the nanowire channel allow a better electrostatic gate control. SiNW transistors of diameters even down to 3 nm have already been prepared, during these years, by various experimental groups [1,2]. Silicon nanowires are quasi-one-dimensional structures in which the electrons are spatially confined in two directions, and they are free to move along the axis of the wire. The nonequilibrium Green’s function (NEGF) formalism is a more advanced transport model for the simulation of SiNW, as it takes into account the wave nature of electrons. The NEGF formalism necessitates rather intensive computational efforts, since it requires detailed information on the propagation of the electron wave packet injected in the device, and microscopic scattering mechanisms other than electron-phonon interactions are difficult to incorporate because the corresponding self-energy terms are usually nonlocal functions of the spatial coordinates.
However, the main quantum transport phenomena in SiNW at room temperature, the source-to-drain tunneling and the conductance fluctuation induced by the quantum interference become significant only when the channel lengths of the SiNW are smaller than 10 nm [3]. Therefore, for longer channels, semiclassical formulations based on the 1D multi-sub-band Boltzmann transport equation (MBTE) can give reliable terminal characteristics when it is solved self-consistently with the 3D Poisson and 2D Schrödinger equations in order to get the self-consistent potential, sub-band energies and wavefunctions. In the literature, the numerical solution of the MBTE has been obtained using Monte Carlo simulations [4,5,6,7,8] or deterministic solvers [6,9] at the expense of large computational effort and statistical noise [10,11,12].
Another alternative is to obtain from the MBTE hydrodynamic models that are a good engineering-oriented approach. This can be achieved by taking the moments of the MBTE and by closing the obtained hierarchy of balance equations, as well as modeling the production terms (i.e., the moments on the collisional operator). The closure problem can be tackled by means of the maximum entropy principle (MEP) of extended thermodynamics [13,14,15], in which the distribution function used to calculate the higher-order moments and the production terms is assumed to be that which maximizes the entropy under the constraints of the given set of moments. We want to underline that the distribution function obtained with the MEP is an approximation of the real one, but from the other side, this distribution is useful to determine analytically without any fitting procedure the higher-order moments and the production terms. In this way, a hydrodynamic model has been obtained with a simplified band structure [16,17,18,19].
The goal of this paper is to set up a consistent hydrodynamic model compatible with the real band structure dictated by atomistic simulation models. The plan of the paper is the following: in Section 2, the electronic band structure is introduced; in Section 3, the quantum confinement is discussed, and the main scattering mechanisms are treated in Section 4; in Section 5, we introduce the kinetic and hydrodynamic models; in Section 6 closure relations are obtained using the MEP; finally, in Section 7, the low-field mobility for a gate-all-around SiNW transistor has been evaluated, and conclusions are drawn in Section 8.

2. Electronic Band Structure

In bulk silicon (Si), the lowest conduction band is formed by six equivalent valleys near the X-point of the Brillouin zone. In this case, we have an indirect band-gap of 1.143 eV at ± 0.85 · 2 π / a 0 in the Δ direction with m t * = 0.19, m l * = 0.98 (units electron mass) and lattice parameter a 0 = 5.43 Å [20]. In SiNW, the band structure is altered with respect to the bulk case depending on the cross-section wire dimension, the atomic configuration and the crystal orientation. The atomistic modeling is able to capture the nanowire band structure, including information about band coupling and mass variations as functions of quantization [3,21,22,23,24]. Atomistic simulations are more appropriate for nanowires of a few nanometer cross-sectional sides, due to the high variability of the involved parameters. However, for SiNW with cross-sections greater than 3 nm, atomistic simulations using the tight-binding (TB) approach show that the band structure is more stable [21]. For a rectangular SiNW with longitudinal direction along the [ 100 ] crystal orientation, confined in the plane ( y , z ) , the 1D Brillouin zone is 1/2 as long as the length of the bulk Si Brillouin zone along the Δ line (i.e., π / a 0 ). The six equivalent Δ conduction valleys of the bulk Si are split into two groups because of the quantum confinement (see Figure 1). The sub-bands related to the four unprimed valleys Δ 4 ( [ 0 ± 10 ] and [ 00 ± 1 ] orthogonal to the wire axis) are projected into a unique valley in the Γ point of the one-dimensional Brillouin zone. Therefore, a SiNW is a direct band-gap semiconductor. The sub-bands related to the primed valleys Δ 2 ( [ ± 100 ] along the wire axis) are found at higher energies and exhibit a minimum, located at k x = ± 0 . 37 π / a 0 (see Figure 2), and the energy gap between the Δ 4 and Δ 2 bottom valley is 117 meV. From the energy dispersion relation E ( k ) obtained from the TB, one can evaluate the effective mass m * in the parabolic spherical band approximation, i.e.,
1 m * = 1 2 2 E k 2
obtaining m Δ 2 * = 0.94, m Δ 4 * = 0.27. Non-parabolic correction to Equation (1) can be introduced [3,25,26], but the fitting parameter depend heavily on the particular atomic configuration.

3. Quantum Confinement

In an ideal quantum wire, the system is assumed sufficiently long along the axis of the wire that translational invariance holds. For a quantum wire with linear expansion in the x-direction and confined in the plane y-z, the normed wave function ϕ ( x , y , z ) can be written in the form:
ϕ ( x , y , z ) = χ l μ ( y , z ) e i k x x L x
where μ is the valley index (one Δ 4 valley and two Δ 2 valleys), l = 1 , N s u b the sub-band index, χ l μ ( y , z ) is the wave function of the l-th sub-band and μ -th valley and the term e i k x x / L x describes an independent plane wave in the x-direction confined to the normalization length, L x 2 x L x 2 , with wave number k x . Let us consider a set of conduction electrons moving in the wire. Any one electron will react to the confining potential U ( y , z ) and to the presence of all other free electrons in the system. The simplest approximation that takes into account the presence of many electrons, called Hartree approximation, is to assume that the electrons as whole produce an average electrostatic potential φ , and that a given electron feels the resulting total potential energy V t o t , i.e.,
V t o t ( x , y , z ) = U ( y , z ) e φ ( x , y , z )
where e is the absolute value of the electric charge. The spatial confinement in the ( y , z ) plane is governed by the Schrödinger–Poisson system (SP):
{ (4a) H [ φ ] χ l x μ [ φ ] = ε l x μ [ φ ] χ l x μ [ φ ] (4b) H [ φ ] = ħ 2 2 m μ * ( 2 y 2 + 2 z 2 ) + V t o t ( x , y , z ) (4c) · [ ϵ 0 ϵ r φ ( x , y , z ) ] = e ( N D N A n [ φ ] ) (4d) n [ φ ] ( x , y , z , t ) = μ l ρ l μ ( x , t ) | χ l x μ [ φ ] ( x , z , t ) | 2
Equation (4a) is the Schrödinger equation in the effective mass approximation (EMA), which has been proved adequate down to the 3-nm wire width [27]. Equation (4c) is the Poisson equation, where ϵ 0 is the absolute dielectric constant (8.854 × 10 12 C/(V·m)); ϵ r is the relative dielectric constant (11.7 for Si, 3.9 for S i O 2 ; and N D , N A are the assigned doping profiles (due to donors and acceptors). The electron density n [ ϕ ] is given by Equation (4d), where ρ l μ ( x , t ) is the linear density in the μ -valley and l-sub-band, which must be evaluated by the transport model (hydrodynamic/kinetic) in the free movement direction. In the following, we shall assume that the cross-section Ω of the wire is surrounded by an oxide, which gives rise to a deep potential barrier,
U ( y , z ) = 0 i f ( y , z ) Ω 3 . 1 e V otherwise .
where 3.1 is the height of the energy barrier at the Si- SiO 2 interface.
The SP system forms a coupled nonlinear Partial Differential Equations, which it is usually solved by an iteration between Poisson and Schrödinger equations. Since a simple iteration by itself does not converge, it is necessary to introduce an adaptive iteration scheme [28]. The solution gives the electrostatic potential φ , as well as the eigenvalues (or sub-band energies) ε l x μ and the eigenvectors (or electron envelope wavefunctions) χ l x μ as a function of the unconfined x direction.
Finally, the total electron energy in the μ -valley and l-sub-band is:
{ (6a) E l μ = ε l x μ + ε ( k x ) + E c (6b) ε ( k x ) [ 1 + α ε ( k x ) ] = ħ 2 k x 2 2 m μ *
where E c is the bottom of the conduction band and α is the non-parabolicity factor (zero in the parabolic case).

4. Carrier Scattering

Essential for the description of any particle transport with a semiclassical approach is the definition of the transition rate w(k, k’), which represents the probability that an electron with wave number vector b , due to a scattering, passes into a state with wave number vector k in the unit time. In the SiNW case:
w ( k , k ) = w ( k x , μ , l , k x , μ , l ) = w ( k x , k x )
representing the probability per unit time that the state ϕ becomes the state ϕ :
ϕ ( x , y , z ) = χ l μ ( y , z ) e i k x x L x ϕ ( x , y , z ) = χ l μ ( y , z ) e i k x x L x .
The transitions rate fulfill Fermi’s golden rule, i.e.,
w ( k x , k x ) = 2 π | ϕ | H p | ϕ | 2 δ E ( ϕ ) E ( ϕ ) ω q
where H p is the perturbation potential, E ( ϕ ) = E l μ is the total energy (6a) and it includes thus the possibility of inter-sub-band scattering processes. The term ± ω q enables that a particle of wave vector q is involved in the scattering process (the plus stands for absorption and the minus for emission).

4.1. Electron-Phonon Scattering

The perturbation potential H p is caused by a deformation potential and writes [29]:
H p = q K q e ± i q · x 2 ρ V c ω ( q ) g ( q ) + 1 2 1 2
where ρ is the mass density, V c the crystal volume, g ( q ) is the phonon distribution function involved in the scattering, ω ( q ) the phonon energy obtained from the phonon dispersion relation for Si (bulk in our case) and:
First order deformation potential:   K q = | q | D 1 Zero order deformation potential:   K q = D o
The scattering mechanism can drive one electron in the same valley (intra-valley) or in a different valley (inter-valley). In bulk silicon, the intra-valley scattering involves only acoustic phonons (two types, LA and TA), and they are evaluated using the first order approximation (with the Debye approximation ω a c ( q ) = v s | q | [20]). The inter-valley scattering is due to six types of phonons: three of the g-type, when electrons scatter between valleys on the same axis, or of the f-type, when the scattering occurs between valleys on perpendicular axes, and they are evaluated in the zero order approximation (using the Einstein approximation with ω ( q ) = constant).
For SiNW, we shall consider bulk phonons and follow the bulk Si scattering selection rules. According to bulk Si scattering selection rules, the elastic processes (due to elastic phonons, surface roughness scattering, impurity scattering) are only intra-valley, whereas the inelastic ones (due to inelastic phonons) are only inter-valley [24].
Since in SiNW, the six conduction valleys reduce to a three-valley model, the valley index μ assumes the value Δ 4 (for the unprimed valleys) and Δ 2 , Δ 2 (for the two primed valleys). However, the two valleys Δ 2 , Δ 2 are symmetric with the same mass (see Figure 2), and then, they can be considered equivalent. Hence, in the following, we shall deal with a two-valley model, the valley A = Δ 4 and the valley B = Δ 2 = Δ 2 . This equivalence in the valleys introduces the following scattering rules
valley A intra - valley , acoustic , elastic A A intra - valley , inelastic A A ( f s c a t t e r i n g , Z i v = 2 ) inter - valley , inelastic A B ( f s c a t t e r i n g , Z i v = 2 )
valley B intra - valley , acoustic , elastic B B inter - valley , inelastic B A ( f s c a t t e r i n g , Z i v = 2 ) intra - valley , inelastic B B ( g s c a t t e r i n g , Z 0 = 1 )
The intra-valley elastic scattering rate (acoustic in the equipartioned case) is given by [30]:
w a c ( k x , μ , l , k x , μ , l ) = s a c G l l μ μ δ E l μ E l μ , s a c = 2 π D a c 2 k B T L ρ v s 2 L x , μ = ( A , B )
For inter-valley scattering, supposing the phonons to be in thermal equilibrium, one obtains:
w i v ( k x , μ , l , k x , μ , l ) = s i v g i v + 1 2 1 2 G l l μ μ δ E l μ E l μ Δ μ μ ω i v
s i v = π D i v 2 ρ ω i v L x Z i v , Δ μ μ = ε μ 0 ε μ 0 , μ , μ = ( A , B )
where ε μ 0 is the bottom of the energy valley, Z i v is the the number of possible final equivalent valleys for the type of inter-valley scattering under consideration, g i v is the Bose–Einstein distribution, i.e.,
g i v = 1 exp ω i v k B T L 1
where T L is the lattice temperature, and G l l μ μ is the form factor:
G l l μ μ = | χ l μ ( y , z ) | 2 | χ l μ ( y , z ) | 2 d y d z .
For intra-valley inelastic scattering, Equation (13) is still valid supposing to have μ = μ , i.e.,
w 0 ( k x , μ , l , k x , μ , l ) = s 0 g 0 + 1 2 1 2 G l l μ μ δ E l μ E l μ ω 0 , s 0 = π D 0 2 ρ ω 0 L x Z 0
All parameters are listed in Table 1.

4.2. Scattering with Impurities

This process is dominant at low temperatures. For an impurity located in the position x i = ( x i , y i , z i ) , with charge Z i e , in the unscreened case, we have:
H p = Z i e 2 4 π ϵ s r , r 2 = ( x x i ) 2 + ( y y i ) 2 + ( z z i ) 2
Then, one obtains [30]:
ϕ | H p | ϕ = Z i e 2 2 π ϵ s L x e i q x x i d z d y ( χ l μ ) χ l μ K 0 | q x | ( y y i ) 2 + ( z z i ) 2
where q x = k x k x , and K 0 ( x ) is the modified Bessel function.
Since this process is elastic and intra-valley, from Equation (8) for the i-th impurity, we obtain:
w i m p ( k x , μ , l , k x , μ , l ) = 2 π Z i e 2 2 π ϵ s L x 2 | H i m p ( k x k x , y i , z i ) | 2 δ E l μ E l μ
where:
H i m p ( k x k x , y i , z i ) = d z d y ( χ l μ ) χ l μ K 0 | k x k x | ( y y i ) 2 + ( z z i ) 2 .
For the sake of simplicity, we shall assume that the impurities are distributed uniformly along the wire, i.e., they are located in ( x , 0 , 0 ) x [ L x / 2 , L x / 2 ] , and in the parabolic band approximation, we get:
w i m p ( k x , μ , l , k x , μ , l ) = Z 2 e 4 n i m μ * 2 π ϵ s 2 L x 3 | H i m p ( k x k x , 0 , 0 ) | 2 H ( a ) a δ ( k x a ) + δ ( k x + a )
a = k x 2 + 2 m μ * 2 ( ε l μ ε l μ )
where n i is the impurities number per unit length and H ( a ) the Heaviside function.

4.3. Surface Roughness Scattering

Surface roughness scattering (SRS) in silicon nanowires is the key scattering mechanism, as it yields a very strong dependence of the low-field electron mobility on the silicon body diameter, as well as on the effective field.
In the case of perfectly smooth Si- SiO 2 interface, the electron wavefunctions and energy level of each sub-band are obtained by solving the Equation (4a) in each x-cross section, with:
V t o t = V e f f + U ( y , z ) , V e f f = e φ ( x , y , z ) + V i m ( y , z ) + V s c ( y , z )
where U is the confining potential, V e f f the effective potential composed by φ the electrostatic potential energy (satisfying the Poisson equation), V i m the image potential due to the mismatch of the dielectric constants between Si and SiO 2 and V s c the exchange-correlation energy due to the electron-electron interaction. However, practically, one must take into account the roughness surface.
Let us consider the wire interface along the x-z plane (whose normal is the y-direction). Then, Δ ( x , z ) is a random function, which describes the deviation of the actual interface from the ideal flat interface. This fluctuation modifies directly the barrier potential U, and it induces a change in the other potentials. Therefore, the sub-band wave functions and energy bands depend explicitly on x.
For the 1D confinement (e.g., the quantum well), the first order complete theory can be found in [31]. In the case of infinite confining potential, the main results is (see Figure 4 in [31]) that for silicon thickness greater than 8 nm, the SRS mobility converges to the SRS bulk mobility, and moreover, the contribution of V i m , V s c can be neglected. Following [4,30], we shall consider a simpler SRS model, where:
  • V i m , V s c are neglected (which of course is true for silicon thickness greater than 8 nm);
  • SRS accounts for deformations only for the potential V e f f (not the wave functions);
  • fluctuation depends only on x, i.e., Δ ( x ) ;
  • the perturbation in the potential is only in the y-variable, i.e., V e f f [ x , y + Δ ( x ) , z ] .
If we expand the potential in the y variable (with x , z constants), up to the first order in Δ ( x ) , we get:
V e f f [ x , y + Δ ( x ) , z ] = V e f f [ x , y + Δ ( x ) , z ] Δ = 0 + V e f f ( x , y , z ) y Δ ( x ) + O ( 2 ) .
Then, finally, the SRS matrix element is:
H s r = V e f f [ x , y + Δ ( x ) , z ] V e f f ( x , y , z ) = e E y ( x , y , z ) Δ ( x )
E y ( x , y , z ) = 1 e V e f f y = φ y
The SRS is an intra-valley scattering mechanism, and it depends on the electric field normal to the surface. Since the cross-section of our wire is a rectangle, we have two contributions, one along the y-direction and the other along the z-direction. Hence,
w s r ( k x , μ , l , k x , μ , l ) = w s r ( k x , μ , l , k x , μ , l , E y ) + w s r ( k x , μ , l , k x , μ , l , E z )
Assuming exponentially correlated surface roughness, in the parabolic band approximation, one obtains [4]:
w s r ( k x , μ , l , k x , μ , l , E y ) = 4 2 e 2 m μ * 3 L x H ( a ) a [ F l l μ μ ( E y ) ] 2 λ s r Δ s r 2 ( k x k x ) 2 λ s r 2 + 2 δ k x a + δ k x + a
where Δ s r and λ s r are the rms (root mean square) height and the correlation length of the fluctuations at the Si- SiO 2 interface, respectively (see Table 1),
a = k x 2 + 2 m μ * ( ε l μ ε l μ ) , F l l μ μ ( E y ) = ( χ l μ ) ( y , z ) E y ( x , y , z ) χ l μ ( y , z ) d y d z
and χ l μ , ε l μ are given by Equation (27).
The contribution for w s r ( k , k , E z ) is similar to Equation (27), supposing to change E y E z .

5. Kinetic and Hydrodynamic Model

In low-dimensional systems, we consider the lattice with perfectly smooth boundaries, free of impurities or other random inhomogeneities. In such an ideal system the energy levels due to disorder are small, so that crystal momentum conservation is approximately preserved. In this framework, we can then construct a kinetic equation in which the distribution function evolves in time under the streaming motion of external forces and spatial gradient, and the randomizing influence of nearly point-like (in space-time) scattering events. For devices with a characteristic length of a few tens of nanometers, the transport of electrons along the axis of the wire can be considered semiclassical within a good approximation; otherwise, a quantum-kinetic approach must be used [32]. The distribution function for the electrons in a quantum wire, with linear expansion in x-direction, depends on the x-direction in real space, the wave vector in x-direction k x and the time t, i.e.,
f l μ = f l μ ( x , k x , t )
for the valley μ and sub-band l. The MBTE reads [30]:
f l μ t + v μ ( k x ) f l μ x e E e f f ( x , t , μ , l ) f l μ k x = η l C η [ f l μ , f l μ ] + η μ μ l C η [ f l μ , f l μ ]
where is the Planck constant divided by 2 π , v μ the electron group velocity and E e f f the effective field:
v μ = 1 E l μ k x = k x m μ * , E e f f ( x , t , μ , l ) = 1 e ε l μ x .
The RHS of Equation (29) is the collisional operator, which is split into two terms modeling respectively intra-valley (with μ = μ ) and inter-valley transitions (with μ μ ). In the low density approximation (not degenerate case), the collisional term for the η-th scattering rate writes:
C η [ f l μ , f l μ ] = L x 2 π d k x w η ( k x , μ , l , k x , μ , l ) f l μ ( x , k x , t ) w η ( k x , μ , l , k x , μ , l ) f l μ ( x , k x , t )
where η = { a c , s r , o } for the acoustic, SR and inelastic (intra-valley) scattering, η = i v for the inter-valley scattering. Starting from the kinetic Equation (29), one can obtain balance equations for macroscopic quantities associated to the flow. By multiplying (29) by a weight function ψ ( k x ) and integrating over 2 / ( 2 π ) d k x , one finds:
M l μ t + 2 ( 2 π ) R ψ v μ f l μ x d k x 2 ( 2 π ) e E e f f R ψ f l μ k x d k x = 2 ( 2 π ) η , l R ψ ( k x ) C η [ f l μ , f l μ ] d k x + 2 ( 2 π ) η , μ μ , l R ψ ( k x ) C η [ f l μ , f l μ ] d k x
where:
M l μ = 2 ( 2 π ) R ψ ( k x ) f l μ ( x , k x , t ) d k x
are the moments relative to the weight functions ψ . By integrating by parts and supposing that f tends to zero sufficiently fast as k , we obtain:
M l μ t + 2 ( 2 π ) x R ψ v μ f l μ d k x + 2 ( 2 π ) e E e f f R f l μ ψ ( k x ) k x d k x = 2 ( 2 π ) η , l R ψ ( k x ) C η [ f l μ , f l μ ] d k x + 2 ( 2 π ) η , μ μ , l R ψ ( k x ) C η [ f l μ , f l μ ] d k x
We have chosen a four-moments model with:
ψ = ( 1 , v μ , E , E v μ ) , v μ = k x m μ * , E = 2 k x 2 2 m μ *
and one obtains from Equation (34) the following balance equations in the unknown ( ρ l μ , V l μ , W l μ , S l μ ):
ρ l μ t + ( ρ l μ V l μ ) x = ρ l μ l C ρ ( μ , l , μ , l ) + ρ l μ l , μ μ C ρ ( μ , l , μ , l )
( ρ l μ V l μ ) t + 2 m μ * ( ρ l μ W l μ ) x + e m μ * ρ l μ E e f f = ρ l μ l C V ( μ , l , μ , l ) + ρ l μ l , μ μ C V ( μ , l , μ , l )
( ρ l μ W l μ ) t + ( ρ l μ S l μ ) x + ρ l μ e E e f f V l μ = ρ l μ l C W ( μ , l , μ , l ) + ρ l μ l , μ μ C W ( μ , l , μ , l )
( ρ l μ S l μ ) t + ( ρ l μ F l μ ) x + 3 e m μ * ρ l μ E e f f W l μ = ρ l μ l C S ( μ , l , μ , l ) + ρ l μ l , μ μ C S ( μ , l , μ , l )
where:
ρ l μ = 2 ( 2 π ) R f l μ ( x , k x , t ) d k x linear electron density , V l μ = 2 ( 2 π ) 1 ρ l μ R f l μ ( x , k x , t ) v μ d k x linear electron velocity , W l μ = 2 ( 2 π ) 1 ρ l μ R f l μ ( x , k x , t ) E d k x linear electron energy , S l μ = 2 ( 2 π ) 1 ρ l μ R f l μ ( x , k x , t ) E v μ d k x linear electron energy flux , F l μ = 2 ( 2 π ) 1 ρ l μ R f l μ ( x , k x , t ) v μ 2 E d k x flux of electron energy flux , C ρ ( μ , l , μ , l ) = 2 ( 2 π ) 1 ρ l μ R C a c [ f l μ , f l μ ] d k x + 2 ( 2 π ) 1 ρ l μ R C s r [ f l μ , f l μ ] d k x , + 2 ( 2 π ) 1 ρ l μ R C o [ f l μ , f l μ ] d k x . C ρ ( μ , l , μ , l ) = 2 ( 2 π ) 1 ρ l μ R C i v [ f l μ , f l μ ] d k x
The production terms for the velocity, energy and energy-flux C V ( μ , l , μ , l ) , C W ( μ , l , μ , l ) , C S ( μ , l , μ , l ) are obtained from C ρ ( μ , l , μ , l ) by multiplying the integrand function for { v μ , E , E v μ } , respectively. From the above definitions, we can introduce the following average quantities:
ρ = μ , l ρ l μ total linear density ,
V = μ , l ρ l μ V l μ ρ mean velocity ,
W = μ , l ρ l μ W l μ ρ mean energy ,
S = μ , l ρ l μ S l μ ρ mean energy flux .

6. Maximum Entropy Principle and Closure Relations

Since the number of unknowns exceeds the number of equations and the production terms are unknown, closure relations must be introduced. The MEP gives a systematic way for obtaining constitutive relations on the basis of information theory [13,14,15]. Such an approach has been used in the simulation of 2D nanoscale structures [33,34] and for simulating the 3D electron transport in sub-micrometric devices, in the case in which the lattice is considered as a thermal bath with constant temperature [35,36,37,38] or when the phonons are off-equilibrium [39,40,41,42,43,44,45]. We shall assume that the electron gas is sufficiently dilute, then the entropy density can be taken as the classical limit of the expression arising in the Fermi statistics, i.e.,
S e = μ , l | χ l μ ( y , z , t ) | 2 S e ( μ , l ) , S e ( μ , l ) = 2 ( 2 π ) k B R ( f l μ log f l μ f l μ ) d k x .
According to the MEP, if a given number of moments M A ( μ , l ) is known, the distribution functions f ^ l μ , which can be used to evaluate the unknown moments, correspond to the extremum of the total entropy density under the constraint that they yield the known moments, i.e.,
2 ( 2 π ) R ψ A ( k x ) f ^ l μ d k x = M A ( μ , l ) .
If we introduce a set of Lagrange multipliers λ A , the problem to maximize S e under the constraints (44) is equivalent to maximizing:
S = S e μ , l A λ A | χ l μ ( y , z , t ) | 2 2 ( 2 π ) R ψ A f ^ l μ d k x M A ( μ , l ) .
So doing, we shall obtain the following distribution function:
f ^ l μ = exp ( Σ ) , Σ = 1 k B A ψ A λ A
ψ A = ( 1 , v μ , E , E v μ ) , λ A = λ l μ , k B λ l μ V , k B λ l μ W , k B λ l μ S
By inserting the previous Equations in (44), we obtain:
M A = M A ( λ A )
which defines implicitly the Lagrange multipliers. In order to invert the above relations, we shall perform an expansion around the thermal equilibrium. In fact, at equilibrium, f ^ must reduce to the Maxwellian. This means:
λ l μ V | E = λ l μ S | E = 0 , λ l μ W | E = 1 k B T L .
Then, we consider the vanishing Lagrange multipliers of higher order with respect to equilibrium, by introducing the smallness parameter τ:
λ l μ V = τ λ ^ l V μ , λ S α = τ λ ^ l S μ .
The inversion problem has been tackled in [16] obtaining, up to the first order in τ (for simplicity, we shall omit the indexes μ , l ):
f ^ = exp λ k B λ W E 1 τ λ ^ V v + λ ^ S v E + O ( τ 2 ) .
The Lagrange multipliers are determined by imposing the constraint (4a):
λ l μ k B = log ρ l μ π 1 2 4 m μ * W l μ , λ l μ W = 1 2 W l μ
λ ^ l V μ = 5 m μ * 4 τ W l μ V l μ + m μ * 4 τ ( W l μ ) 2 S l μ , λ ^ l S μ = m μ * 4 τ ( W l μ ) 2 V l μ m μ * 12 τ ( W l μ ) 3 S l μ
whereas the higher order flux is:
F l μ = 6 ( W l μ ) 2 m μ * .
In order to close the system, we need functional relations for the production terms, which can be evaluated by using the MEP distribution Function (49).

6.1. Closure for the Electron Number Production

C ρ ( μ , l , μ , l ) = C ρ ( μ , l , μ , l ) ( a c ) + C ρ ( μ , l , μ , l ) ( s r ) + C ρ ( μ , l , μ , l ) ( o )
where:
C ρ ( μ , l , μ , l ) ( a c ) = 2 ( 2 π ) R C a c [ f l μ , f l μ ] d k x
C ρ ( μ , l , μ , l ) ( s r ) = 2 ( 2 π ) R C s r [ f l μ , f l μ ] d k x
C ρ ( μ , l , μ , l ) ( o ) = 2 ( 2 π ) R C o [ f l μ , f l μ ] d k x
and:
C ρ ( μ , l , μ , l ) ( i v ) = 2 ( 2 π ) R C i v [ f l μ , f l μ ] d k x

6.1.1. Evaluation of C ρ ( μ , l , μ , l ) ( o )

This is an intra-valley inelastic scattering mechanism, where the scattering rate is given by Equation (17). After long calculations, one obtains:
ρ l μ C ρ ( μ , l , μ , l ) ( o ) = s o p π L x 2 π g 0 G l l μ μ A 1 l l + μ μ ( Δ l l + μ μ ) + ( g 0 + 1 ) G l l μ μ A 1 l l μ μ ( Δ l l μ μ ) g 0 G l l μ μ A 2 l l μ μ ( Δ l l μ μ ) ( g 0 + 1 ) G l l μ μ A 2 l l + μ μ ( Δ l l + μ μ )
A 1 l l ± μ μ ( Δ ) = ρ l μ 2 π W l μ exp Δ 2 W l μ 0 + d k x H 2 ( k x ) 2 2 m μ * Δ 2 k x 2 2 m μ * Δ exp 2 k x 2 4 m μ * W l μ
A 2 l l ± μ μ ( Δ ) = ρ l μ 2 π W l μ 0 + H 2 ( k x ) 2 2 m μ * Δ l l ± μ μ 2 ( k x ) 2 2 m μ * Δ l l ± μ μ exp 2 k x 2 4 m μ * W l μ d k x
Δ l l ± μ μ = ε l μ ε l μ ± ω 0

6.1.2. Evaluation of C ρ ( μ , l , μ , l ) ( a c )

This is an intra-valley elastic scattering mechanism, where the scattering rate is given by Equation (12). We obtain:
ρ l μ C ρ ( μ , l , μ , l ) ( a c ) = s a c π L x 2 π G l l μ μ J 1 l l μ μ G l l μ μ J 2 l l μ μ
where:
J 1 l l μ μ = ρ l μ 2 π W l μ exp Δ l l μ μ 2 W l μ 0 + H 2 ( k x ) 2 2 m μ * Δ l l μ μ 2 k x 2 2 m μ * Δ l l μ μ exp 2 k x 2 4 m μ * W l μ d k x
J 2 l l μ μ = ρ l μ 2 π W l μ 0 + H 2 ( k x ) 2 2 m μ * Δ l l μ μ 2 ( k x ) 2 2 m μ * Δ l l μ μ exp 2 k x 2 4 m μ * W l μ d k x
and:
Δ l l μ μ = ε l μ ε l μ .

6.1.3. Evaluation of C ρ ( μ , l , μ , l ) ( s r )

This is an intra-valley elastic scattering mechanism, where the scattering rate is given by: Equation (27). We get:
ρ l μ C ρ ( μ , l , μ , l ) ( s r ) = 2 e 2 λ s r Δ s r 2 m μ * 2 2 π W l μ × | F l l μ μ | 2 ρ l μ I s ( Δ l l μ μ , W l μ ) + I s + ( Δ l l μ μ , W l μ ) | F l l μ μ | 2 ρ l μ I s ( Δ l l μ μ , W l μ ) + I s + ( Δ l l μ μ , W l μ )
where:
I s ± ( Δ l l μ μ , W l μ ) = 0 + H E + Δ l l μ μ 2 m μ * λ s r 2 2 E ± E + Δ l l μ μ 2 + 2 exp E 2 m μ * W l μ E [ E + Δ l l μ μ ] d E

6.1.4. Evaluation of C ρ ( μ , l , μ , l ) ( i m p )

This is an intra-valley scattering, where the scattering mechanism is given by Equation (22). We get:
ρ l μ C ρ ( μ , l , μ , l ) ( i m p ) = Ω 0 ( W l μ , W l μ ) + Ω 1 ( W l μ ) V l μ Ω 2 ( W l μ ) V l μ + Ω 3 ( W l μ ) S l μ Ω 4 ( W l μ ) S l μ
where:
Ω 0 ( W l μ , W l μ ) = Z 2 e 4 n i m μ * 2 π 3 ϵ s 2 3 A 1 l μ + A 1 l + μ B 1 l μ B 1 l + μ
Ω 1 ( W l μ ) = Z 2 e 4 n i m μ * 2 π 3 ϵ s 2 3 5 m μ * 4 W l μ A 2 l μ + A 2 l + μ + m μ * 4 ( W l μ ) 2 A 3 l μ + A 3 l + μ
Ω 2 ( W l μ ) = Z 2 e 4 n i m μ * 2 π 3 ϵ s 2 3 5 m μ * 4 W l μ B 2 l μ + B 2 l + μ + m μ * 4 ( W l μ ) 2 B 3 l μ + B 3 l + μ
Ω 3 ( W l μ ) = Z 2 e 4 n i m μ * 2 π 3 ϵ s 2 3 m μ * 4 ( W l μ ) 2 A 2 l μ + A 2 l + μ m μ * 12 ( W l μ ) 3 A 3 l μ + A 3 l + μ
Ω 4 ( W l μ ) = Z 2 e 4 n i m μ * 2 π 3 ϵ s 2 3 m μ * 4 ( W l μ ) 2 B 2 l μ + B 2 l + μ m μ * 12 ( W l μ ) 3 B 3 l μ + B 3 l + μ
A 1 l ± μ = ρ l μ π 4 m μ * W l μ + exp 2 k x 2 4 m μ * W l μ H ( a ) a H i m p ( k x ± a , 0 , 0 ) 2 d k x
A 2 l ± μ = ρ l μ π 4 m μ * W l μ + exp 2 k x 2 4 m μ * W l μ k z m μ * H ( a ) a H i m p ( k x ± a , 0 , 0 ) 2 d k x
A 3 l ± μ = ρ l μ π 4 m μ * W l μ + exp 2 k x 2 4 m μ * W l μ 3 k x 3 2 ( m μ * ) 2 H ( a ) a H i m p ( k x ± a , 0 , 0 ) 2 d k x
B 1 l ± μ = ρ l μ π 4 m μ * W l μ + exp 2 k x 2 4 m μ * W l μ H ( a ) a H i m p ( k x ± a , 0 , 0 ) 2 d k x
B 2 l ± μ = ρ l μ π 4 m μ * W l μ + exp 2 k x 2 4 m μ * W l μ k x m μ * H ( a ) a H i m p ( k x ± a , 0 , 0 ) 2 d k x
B 3 l ± μ = ρ l μ π 4 m μ * W l μ + exp 2 k x 2 4 m μ * W l μ 3 k x 3 2 ( m μ * ) 2 H ( a ) a H i m p ( k x ± a , 0 , 0 ) 2 d k x
H i m p ( k x ± a , 0 , 0 ) = d y d z ( χ l μ ) χ l μ K 0 | k x a | y 2 + z 2
a = k x 2 + 2 m μ * ( ε l μ ε l μ ) , a = k x 2 + 2 m μ * ( ε l μ ε l μ )

6.1.5. Evaluation of C ρ ( μ , l , μ , l ) ( i v )

This is an inter-valley inelastic scattering mechanism, where the scattering rate is given by Equation (13). The result is similar to Equation (58) (obtained for an intra-valley, inelastic scattering) with μ μ , but supposing to change Equation (61) into:
Δ ˜ l l ± μ μ = ε l μ ε l μ ± ω i v + Δ μ μ
where Δ μ μ is given in Equation (14). Then, we obtain:
ρ l μ C ρ ( μ , l , μ , l ) ( i v ) = s i v π L x 2 π g 0 G l l μ μ A 1 l l + μ μ ( Δ ˜ l l + μ μ ) + ( g 0 + 1 ) G l l μ μ A 1 l l μ μ ( Δ ˜ l l μ μ ) g 0 G l l μ μ A 2 l l μ μ ( Δ ˜ l l + μ μ ) ( g 0 + 1 ) G l l μ μ A 2 l l + μ μ ( Δ ˜ l l + μ μ )
A 1 l l ± μ μ ( Δ ) = ρ l μ 2 π W l μ exp Δ 2 W l μ 0 + d k x H 2 ( k x ) 2 2 m μ * Δ 2 k x 2 2 m μ * Δ exp 2 k x 2 4 m μ * W l μ
A 2 l l ± μ μ ( Δ ) = ρ l μ 2 π W l μ 0 + H 2 ( k x ) 2 2 m μ * Δ 2 ( k x ) 2 2 m μ * Δ exp 2 k x 2 4 m μ * W l μ d k x

6.2. Closure for the Production of Electron Energy

C W ( μ , l , μ , l ) = C W ( μ , l , μ , l ) ( a c ) + C W ( μ , l , μ , l ) ( s r ) + C W ( μ , l , μ , l ) ( o )
where:
C W ( μ , l , μ , l ) ( a c ) = 2 ( 2 π ) R E C a c [ f l μ , f l μ ] d k x
C W ( μ , l , μ , l ) ( s r ) = 2 ( 2 π ) R E C s r [ f l μ , f l μ ] d k x
C W ( μ , l , μ , l ) ( o ) = 2 ( 2 π ) R E C o [ f l μ , f l μ ] d k x
and:
C W ( μ , l , μ , l ) ( i v ) = 2 ( 2 π ) R E C i v [ f l μ , f l μ ] d k x
We observe that, with respect to Equations (54)–(57), there is an extra E , and by multiplying all of the previous integrals by E , similar results hold.

6.3. Closure for the Production of Electron Crystal Momentum

C V ( μ , l , μ , l ) = C V ( μ , l , μ , l ) ( a c ) + C V ( μ , l , μ , l ) ( s r ) + C V ( μ , l , μ , l ) ( o )
where:
C V ( μ , l , μ , l ) ( a c ) = 2 ( 2 π ) R k x m μ * C a c [ f l μ , f l μ ] d k x
C V ( μ , l , μ , l ) ( s r ) = 2 ( 2 π ) R k x m μ * C s r [ f l μ , f l μ ] d k x
C V ( μ , l , μ , l ) ( o ) = 2 ( 2 π ) R k x m μ * C o [ f l μ , f l μ ] d k x
and:
C V ( μ , l , μ , l ) ( i v ) = 2 ( 2 π ) R k x m μ * C i v [ f l μ , f l μ ] d k x

6.3.1. Evaluation of C V ( μ , l , μ , l ) ( o )

For this production term, the scattering rate is given by Equation (17). We get:
ρ l μ C V ( μ , l , μ , l ) ( o ) = Ω 1 l l μ μ ( W l μ ) ρ l μ V l μ + Ω 2 l l μ μ ( W l μ ) ρ l μ S l μ
where:
Ω 1 l l μ μ ( W l μ ) = s o p π L x 2 π G l l μ μ 5 4 W l μ ( I V l l μ μ + I V l l + μ μ ) + 1 4 ( W l μ ) 2 ( I S l l μ μ + I S l l + μ μ )
Ω 2 l l μ μ ( W l μ ) = s o p π L x 2 π G l l μ μ 1 4 ( W l μ ) 2 ( I V l l μ μ + I V l l + μ μ ) 1 12 ( W l μ ) 3 ( I S l l μ μ + I S l l + μ μ )
I V l l ± μ μ = m μ * 2 π W l μ g 0 + 1 2 ± 1 2 0 + k x 2 H 2 ( k x ) 2 2 m μ * Δ l l ± μ μ 2 ( k x ) 2 2 m μ * Δ l l ± μ μ exp 2 k x 2 4 m μ * W l μ d k x
I S l l ± μ μ = 3 2 m μ * 2 2 π W l μ g 0 + 1 2 ± 1 2 0 + k x 4 H 2 ( k x ) 2 2 m μ * Δ l l ± μ μ 2 ( k x ) 2 2 m μ * Δ l l ± μ μ exp 2 k x 2 4 m μ * W l μ d k x

6.3.2. Evaluation of C V ( μ , l , μ , l ) ( a c )

This is an intra-valley inelastic scattering rate, where the scattering rate is given by Equation (12). We get:
ρ l μ C V ( μ , l , μ , l ) ( a c ) = Δ 1 l l μ μ ( W l μ ) ρ l μ V l μ + Δ 2 l l μ μ ( W l μ ) ρ l μ S l μ
where:
Δ 1 l l μ μ ( W l μ ) = s a c π L x 2 π W l μ G l l μ μ 5 4 F l l μ μ + 1 4 W l μ L l l μ μ
Δ 2 l l μ μ ( W l μ ) = s a c π L x 2 π W l μ G l l μ μ 1 4 W l μ F l l μ μ 1 12 ( W l μ ) 2 L l l μ μ
F l l μ μ ( W l μ ) = 2 π W l μ 0 + H 2 ( k x ) 2 2 m μ * Δ l l μ μ 2 ( k x ) 2 2 m μ * Δ l l μ μ exp 2 k x 2 4 m μ * W l μ k z 2 m μ * d k x
L l l μ μ ( W l μ ) = 2 π W l μ 0 + H 2 ( k x ) 2 2 m μ * Δ l l μ μ 2 ( k x ) 2 2 m μ * Δ l l μ μ exp 2 k x 2 4 m μ * W l μ 3 k z 4 2 m μ * 2 d k x

6.3.3. Evaluation of C V ( μ , l , μ , l ) ( s r )

This is an intra-valley elastic scattering mechanism, and from Equation (27), we get:
C V ( μ , l , μ , l ) ( s r ) = m μ * e 2 λ s r Δ s r 2 6 2 π ( W l μ ) 5 | F l l μ μ | 2 Ω V V l μ + Ω S S l μ m μ * e 2 λ s r Δ s r 2 6 2 π ( W l μ ) 5 | F l l μ μ | 2 Ω ˜ V V l μ + Ω ˜ S S l μ
where:
Ω V = 15 8 π ( 2 W l μ ) 7 2 I S 2 ( Δ l l μ μ , W l μ , j = 0 ) + I S 2 + ( Δ l l μ μ , W l μ , j = 0 ) 3 4 π ( 2 W l μ ) 5 2 I S 2 ( Δ l l μ μ , W l μ , j = 1 ) + I S 2 + ( Δ l l μ μ , W l μ , j = 1 )
Ω ˜ V = 15 8 π ( 2 W l μ ) 7 2 I ˜ S 2 ( Δ l l μ μ , W l μ , j = 0 ) + I ˜ S 2 + ( Δ l l μ μ , W l μ , j = 0 ) 3 4 π ( 2 W l μ ) 5 2 I ˜ S 2 ( Δ l l μ μ , W l μ , j = 1 ) + I ˜ S 2 + ( Δ l l μ μ , W l μ , j = 1 )
Ω S = 3 4 π ( 2 W l μ ) 5 2 I S 2 ( Δ l l μ μ , W l μ , j = 0 ) + I S 2 + ( Δ l l μ μ , W l μ , j = 0 ) + π 2 ( 2 W l μ ) 3 2 I S 2 ( Δ l l μ μ , W l μ , j = 1 ) + I S 2 + ( Δ l l μ μ , W l μ , j = 1 )
Ω ˜ S = 3 4 π ( 2 W l μ ) 5 2 I ˜ S 2 ( Δ l l μ μ , W l μ , j = 0 ) + I ˜ S 2 + ( Δ l l μ μ , W l μ , j = 0 ) + π 2 ( 2 W l μ ) 3 2 I ˜ S 2 ( Δ l l μ μ , W l μ , j = 1 ) + I ˜ S 2 + ( Δ l l μ μ , W l μ , j = 1 )
I S 2 ± ( Δ l l μ μ , W l μ , j ) = 0 + E j + 1 2 H E Δ l l μ μ 2 m μ * λ s r 2 2 E ± E Δ l l μ μ 2 + 2 exp E 2 m μ * W l μ E Δ l l μ μ d E
I ˜ S 2 ± ( Δ l l μ μ , W l μ , j ) = 0 + E j H E Δ l l μ μ 2 m μ * λ s r 2 2 E ± E Δ l l μ μ 2 + 2 exp E 2 m μ * W l μ d E

6.3.4. Evaluation of C V ( μ , l , μ , l ) ( i m p )

This is an intra-valley scattering, and from Equation (22), we get:
ρ l μ C V ( μ , l , μ , l ) ( i m p ) = Ω ˜ 0 ( W l μ , W l μ ) + Ω ˜ 1 ( W l μ ) V l μ Ω ˜ 2 ( W l μ ) V l μ + Ω ˜ 3 ( W l μ ) S l μ Ω ˜ 4 ( W l μ ) S l μ
where:
Ω ˜ 0 ( W l μ , W l μ ) = Z 2 e 4 n i 4 π 3 ϵ s 2 2 A ˜ 1 l μ + A ˜ 1 l + μ B ˜ 1 l μ B ˜ 1 l + μ
Ω ˜ 1 ( W l μ ) = Z 2 e 4 n i 4 π 3 ϵ s 2 2 5 m μ * 4 W l μ A ˜ 2 l μ + A ˜ 2 l + μ + m μ * 4 ( W l μ ) 2 A ˜ 3 l μ + A ˜ 3 l + μ
Ω ˜ 2 ( W l μ ) = Z 2 e 4 n i 4 π 3 ϵ s 2 2 5 m μ * 4 W l μ B ˜ 2 l μ + B ˜ 2 l + μ + m μ * 4 ( W l μ ) 2 B ˜ 3 l μ + B ˜ 3 l + μ
Ω ˜ 3 ( W l μ ) = Z 2 e 4 n i 4 π 3 ϵ s 2 2 m μ * 4 ( W l μ ) 2 A ˜ 2 l μ + A ˜ 2 l + μ m μ * 12 ( W l μ ) 3 A ˜ 3 l μ + A ˜ 3 l + μ
Ω ˜ 4 ( W l μ ) = Z 2 e 4 n i 4 π 3 ϵ s 2 2 m μ * 4 ( W l μ ) 2 B ˜ 2 l μ + B ˜ 2 l + μ m μ * 12 ( W l μ ) 3 B ˜ 3 l μ + B ˜ 3 l + μ
A ˜ 1 l ± μ = ρ l μ π 4 m μ * W l μ + exp 2 k x 2 4 m μ * W l μ H ( a ) H i m p ( k x ± a , 0 , 0 ) 2 d k x
A ˜ 2 l ± μ = ρ l μ π 4 m μ * W l μ + k z m μ * exp 2 k x 2 4 m μ * W l μ H ( a ) H i m p ( k x ± a , 0 , 0 ) 2 d k x
A ˜ 3 l ± μ = ρ l μ π 4 m μ * W l μ + 3 k x 3 2 ( m μ * ) 2 exp 2 k x 2 4 m μ * W l μ H ( a ) H i m p ( k x ± a , 0 , 0 ) 2 d k x
B ˜ 1 l ± μ = ρ l μ π 4 m μ * W l μ + k x exp 2 k x 2 4 m μ * W l μ H ( a ) a H i m p ( k x ± a , 0 , 0 ) 2 d k x
B ˜ 2 l ± μ = ρ l μ π 4 m μ * W l μ + k x 2 m μ * exp 2 k x 2 4 m μ * W l μ H ( a ) a H i m p ( k x ± a , 0 , 0 ) 2 d k x
B ˜ 3 l ± μ = ρ l μ π 4 m μ * W l μ + 3 k x 4 2 ( m μ * ) 2 exp 2 k x 2 4 m μ * W l μ H ( a ) a H i m p ( k x ± a , 0 , 0 ) 2 d k x
H i m p ( k x ± a , 0 , 0 ) = d y d z ( χ l μ ) χ l μ K 0 | k x a | y 2 + z 2
a = k x 2 + 2 m μ * ( ε l μ ε l μ ) , a = k x 2 + 2 m μ * ( ε l μ ε l μ )

6.3.5. Evaluation of C V ( μ , l , μ , l ) ( i v )

This term is similar to Equation (96) (obtained for an intra-valley inelastic scattering), but with μ μ and supposing to change Equation (61) into:
Δ ˜ l l ± μ μ = ε l μ ε l μ ± ω i v + Δ μ μ
where Δ μ μ is given in Equation (14), we obtain:
ρ l μ C V ( μ , l , μ , l ) ( i v ) = Ω 1 l l μ μ ( W l μ ) ρ l μ V l μ + Ω 2 l l μ μ ( W l μ ) ρ l μ S l μ
where:
Ω 1 l l μ μ ( W l μ ) = s o p π L x 2 π G l l μ μ 5 4 W l μ ( I V l l μ μ + I V l l + μ μ ) + 1 4 ( W l μ ) 2 ( I S l l μ μ + I S l l + μ μ )
Ω 2 l l μ μ ( W l μ ) = s o p π L x 2 π G l l μ μ 1 4 ( W l μ ) 2 ( I V l l μ μ + I V l l + μ μ ) 1 12 ( W l μ ) 3 ( I S l l μ μ + I S l l + μ μ )
I V l l ± μ μ = m μ * 2 π W l μ g 0 + 1 2 ± 1 2 0 + k x 2 H 2 ( k x ) 2 2 m μ * Δ l l ± μ μ 2 ( k x ) 2 2 m μ * Δ l l ± μ μ exp 2 k x 2 4 m μ * W l μ d k x
I S l l ± μ μ = 3 2 m μ * 2 2 π W l μ g 0 + 1 2 ± 1 2 0 + k x 4 H 2 ( k x ) 2 2 m μ * Δ l l ± μ μ 2 ( k x ) 2 2 m μ * Δ l l ± μ μ exp 2 k x 2 4 m μ * W l μ d k x

6.4. Closure for the Production of Electron Energy-Flux

C S ( μ , l , μ , l ) = C S ( μ , l , μ , l ) ( a c ) + C S ( μ , l , μ , l ) ( s r ) + C S ( μ , l , μ , l ) ( o )
where:
C S ( μ , l , μ , l ) ( a c ) = 2 ( 2 π ) R v μ E C a c [ f l μ , f l μ ] d k x
C S ( μ , l , μ , l ) ( s r ) = 2 ( 2 π ) R v μ E C s r [ f l μ , f l μ ] d k x
C S ( μ , l , μ , l ) ( o ) = 2 ( 2 π ) R v μ E C o [ f l μ , f l μ ] d k x
and:
C S ( μ , l , μ , l ) ( i v ) = 2 ( 2 π ) R v μ E C i v [ f l μ , f l μ ] d k x
We observe that, with respect to Equations (92)–(95), there is an extra E , and by multiplying all of the previous integrals by E , similar results hold.

7. Case of Study

As a case of study, we introduce the so-called gate-all-around (GAA) SiNW transistor. This is a silicon nanowire with an added gate wrapped around it, in such a way that we have a three contact device with source, drain and gate (see Figure 3). Such devices have been designed during these years in order to maintain a good electrostatic control in the channel [1,2]. The gate electrode is assumed to be metallic, so that there is no potential drop inside the gate, and depletion effects are not considered. In the following, we shall consider a very simple SiNW transistor having the channel homogeneously doped to N D = 3 × 10 15 cm 3 and very long ( L x = 120 nm) with respect to the transversal dimensions ( L y = L z = 10 nm), where the oxide thickness ( SiO 2 ) is tox = 1 nm. In such a case, the moment system reduces to a set of ordinary differential equations with the time as the only independent variable to be coupled to the SP system Equation (4).
First of all, let us consider the thermal equilibrium regime where no voltage is applied to the contacts, i.e., V S = V D = V G = 0 and no current flows. Hence, the electron distribution function is the Maxwellian:
f l μ ( e q ) ( k x ) = N 0 exp 2 k x 2 2 m μ * + ε l μ + ε μ 0 ν k B T
where ν is the Fermi level, ε μ 0 the valley energy minimum, T the electron temperature, which we shall assume to be the same in each sub-band and equal to the lattice temperature. The condition of zero net current requires that the Fermi level must be constant throughout the sample, and it can be determined by imposing that the total electron number equals the total donor number in the wire. Then, the linear electron density at equilibrium writes:
ρ l μ ( e q ) ( x ) = N D L y L z m μ * Z ( e q ) exp ε l x μ ( e q ) ε μ 0 k B T , Z ( e q ) = μ , l m μ * exp ε l x μ ( e q ) ε μ 0 k B T
where the sub-band energies ε l x μ ( e q ) are obtained by solving the SP system Equation (4) using Equation (139) with V D = V S = V G = 0.
Now, we consider the quasi-equilibrium regime, where a very small axial electric field frozen along the channel (1000 V/cm) is applied, and we turn on the gate. The system is still in local thermal equilibrium, the distribution function is the Maxwellian, but some charge flows in the wire. The linear density writes:
ρ l μ = N D L y L z m μ * Z ( e q ) exp ε l x μ ε μ 0 k B T
where the only difference between Equations (139) and (140) is in the energy sub-bands ε l x μ , which now are obtained solving the SP system Equation (4) with V S = 0 . 012 V, V D = 0 V, V G = 1 V. Once the solution has been obtained, the energies ε l μ and wave functions χ l μ for each sub-band are fixed and exported into the hydrodynamic model. Moreover, the obtained linear density ρ l μ is used as the initial condition for this model. The other initial conditions for the balance Equations (35)–(38) are:
V l μ = 0 , W l μ = 1 2 k B T , S l μ = 0 .
Figure 4 shows the distribution of charge density (4d) and the total potential (3) along the cross-section at x = 48 nm and z = 0 nm. A surface inversion layer is formed, similar to a usual MOSFET channel with the electron density peak 1 nm from the oxide interface. The band-bending at the interface forms the quantum well for the electrons.
In Figure 5, the total linear density for the A and B valleys (i.e., ρ A = l ρ l A , ρ B = l ρ l B ) versus the simulation time is plotted. One can note a depletion of the higher B valley. In fact, the electrons change valley according to the inter-valley scattering mechanism given by Equation (13). The scattering from A to B happens if the electron energy is greater than a certain threshold (related to the energy gap between the two valleys, (i.e., 117 meV), whereas the converse is more probable. If we apply a small electric field (1000 V/cm) in the wire, the electrons, at the beginning of the simulation, cannot gain enough energy to jump from A to B, whereas it is more probable that the opposite will happen. As the simulation time increases, the electrons in A gain enough energy (i.e., they have a slightly increase of the mean energy) to activate the jump from A to B, and in the stationary regime, an equilibrium state is reached.
In Figure 6, Figure 7 and Figure 8, we show the mean velocity (40), energy (41) and energy-flux (42), respectively, versus the simulation time, obtained with and without the SRS mechanism. In these figures, one observes that the stationary regime is reached in a few picoseconds for the velocity and the energy flux, in about ten picoseconds for the energy, and the dependence on the SRS mechanism is clearly understood.
Finally, we have computed the low-field mobility μ l o w , which is a fundamental parameter for engineering applications. It is defined as the ratio between the average electron velocity, evaluated in the stationary regime, and the driving field ( E = 1 kV/cm), i.e.,
μ l o w = μ A ρ A + μ B ρ B ρ A + ρ B , μ A = l V l A E , μ B = l V l B E
where μ A , μ B are the mobilities in the respective valleys.
In Table 2, the results for the A and B valley mobility, obtained with and without the SRS mechanism, are presented. The difference between these values proves that the SRS is a key mechanism in the SiNW device performance. In the above table, we can notice that mobility in the A-valley is bigger with respect to that obtained in the B-valley. Since the mobility depends (inversely) on the effective mass, the valley splitting reduces the mobility along the axis of the wire (in the B-valley where the effective mass is 0.94), but quantum confinement increases mobility in the transverse direction (in the A-valley where the effective mass is 0.27).

8. Conclusions

Charge transport phenomena in SiNW can be treated by coupling the Schrödinger–Poisson system (governing the spatial confinement) with a hydrodynamic model (governing the transport along the free motion direction). The hydrodynamic model has been derived from the multi-sub-band Boltzmann transport equations using the MEP in order to obtain a closed system of PDEs. An appropriate electronic band structure obtained by tight-binding simulations has been included, as well as the main scattering mechanisms. This model has been used to evaluate the low-field mobility in a very simple gate-all-around SiNW transistor. The inclusion of high doping effects in the model and simulation of SiNW transistors will be the topics of future research.

Acknowledgments

We acknowledge the support of the Università degli Studi di Catania, FIR 2014 “Charge Transport in Graphene and Low dimensional Structures: modeling and simulation” and the National Group of Mathematical Physics (GNFM-INdAM), “Progetto giovani 2015”.

Author Contributions

The two authors contributed equally to this work. Both authors have read and approved the final manuscript.

Conflicts of Interest

The authors declare no conflict of interest.

References

  1. Singh, N.; Agarwal, A.; Bera, L.K.; Liow, T.Y.; Yang, R.; Rustagi, S.C.; Tung, C.H.; Kumar, R.; Lo, G.Q.; Balasubramanian, N.; et al. High-performance fully depleted silicon nanowire (diameter ≤ 5 nm) gate-all-around CMOS devices. IEEE Electron Device Lett. 2006, 27, 383–386. [Google Scholar] [CrossRef]
  2. Guerfi, Y.; Larrieu, G. Vertical Silicon Nanowire Field Effect Transistors with Nanoscale Gate-All-Around. Nanoscale Res. Lett. 2016, 11. [Google Scholar] [CrossRef] [PubMed]
  3. Wang, J.; Lundstrom, M. Does source-to-drain tunneling limit the ultimate scaling of MOSFETs? In Proceedings of the International Electron Devices Meeting, San Francisco, CA, USA, 8–11 December 2002; pp. 707–710.
  4. Ramayya, E.B.; Vasileska, D.; Goodnick, S.M.; Knezevic, I. Electron transport in silicon nanowires: The role of acoustic phonon confinement and surface roughness scattering. J. Appl. Phys. 2008, 104, 063711. [Google Scholar] [CrossRef]
  5. Ramayya, E.B.; Knezevic, I. Self-consistent Poisson-Schrödinger-Monte Carlo Solver: Electron Mobility in Silicon Nanowires. J. Comput. Electron. 2010, 9, 206–210. [Google Scholar] [CrossRef]
  6. Lenzi, M.; Palestri, P.; Gnani, E.; Reggiani, S.; Gnudi, A.; Esseni, D.; Selmi, L.; Baccarani, G. Investigation of the Transport Properties of Silicon Nanowires Using Deterministic and Monte Carlo Approaches to the Solution of the Boltzmann Transport Equation. IEEE Trans. Electron Devices 2008, 55, 2086–2096. [Google Scholar] [CrossRef]
  7. Davoody, A.H.; Ramayya, E.B.; Maurer, L.N.; Knezevic, I. Ultrathin GaN nanowires: Electronic, thermal, and thermoelectric properties. Phys. Rev. B 2014, 89, 115313. [Google Scholar] [CrossRef]
  8. Ryu, H. A multi-sub-band Monte Carlo study on dominance of scattering mechanisms over carrier transport in sub-10-nm Si nanowire FETs. Nanoscale Res. Lett. 2016, 11. [Google Scholar] [CrossRef] [PubMed]
  9. Ossig, G.; Schürrer, F. Simulation of non-equilibrium electron transport in silicon quantum wires. J. Comput. Electron. 2008, 7, 367–370. [Google Scholar] [CrossRef]
  10. Muscato, O.; Wagner, W.; Di Stefano, V. Numerical study of the systematic error in Monte Carlo schemes for semiconductors. ESAIM 2010, 44, 1049–1068. [Google Scholar] [CrossRef]
  11. Muscato, O.; Wagner, W.; Di Stefano, V. Properties of the steady state distribution of electrons in semiconductors. Kinet. Relat. Models 2011, 4, 809–829. [Google Scholar] [CrossRef]
  12. Muscato, O.; Di Stefano, V.; Wagner, W. A variance-reduced electrothermal Monte Carlo method for semicond. Comput. Math. Appl. 2013, 65, 520–527. [Google Scholar] [CrossRef]
  13. Lebon, G.; Jou, D.; Casas-Vázquez, J. Understanding Non-Equilibrium Thermodynamics; Springer: Berlin/Heidelberg, Germany, 2008. [Google Scholar]
  14. Mascali, G.; Romano, V. Hydrodynamic sub-band model for semiconductors based on the maximum entropy principle. Nuovo Cimento. C 2010, 33, 155–163. [Google Scholar]
  15. Mascali, G.; Romano, V. A non parabolic hydrodynamic sub-band model for semiconductors based on the maximum entropy principle. Math. Comput. Model. 2012, 55, 1003–1020. [Google Scholar] [CrossRef]
  16. Muscato, O.; Di Stefano, V. Hydrodynamic modeling of silicon quantum wires. J. Comput. Electron. 2012, 11, 45–55. [Google Scholar] [CrossRef]
  17. Muscato, O.; Di Stefano, V. Hydrodynamic simulation of a n+nn+ silicon nanowire. Contin. Mech. Thermodyn. 2014, 26, 197–205. [Google Scholar] [CrossRef]
  18. Castiglione, T.; Muscato, O. Non-parabolic band hydrodynamic model for silicon quantum wires. J. Comput. Theor. Transp 2016, in press. [Google Scholar]
  19. Muscato, O.; Castiglione, T. Electron transport in silicon nanowires having different cross-sections. Commun. Appl. Ind. Math. 2016, 7, 8–25. [Google Scholar] [CrossRef]
  20. Jacoboni, C.; Reggiani, L. The Monte Carlo method for the solution of charge transport in semiconductors with applications to covalent materials. Rev. Mod. Phys. 1983, 55, 645–705. [Google Scholar] [CrossRef]
  21. Zheng, Y.; Rivas, C.; Lake, R.; Alam, K.; Boykin, T.B.; Klimeck, G. Electronic Properties of Silicon Nanowires. IEEE Trans. Electron Devices 2005, 52, 1097–1103. [Google Scholar] [CrossRef]
  22. Nehari, K.; Cavassilas, N.; Autran, J.L.; Bescond, M.; Munteanu, D.; Lannoo, M. Influence of band structure on electron ballistic transport in silicon nanowire MOSFET’s: An atomistic study. Solid-State Electron. 2006, 50, 716–721. [Google Scholar] [CrossRef]
  23. Neophytou, N.; Paul, A.; Lundstrom, M.S.; Klimeck, G. Bandstructure Effects in Silicon Nanowire Electron Transport. IEEE Trans. Electron Devices 2008, 55, 1286–1297. [Google Scholar] [CrossRef]
  24. Neophytou, N.; Kosina, H. Atomistic simulations of low-field mobility in Si nanowires: Influence of confinement and orientation. Phys. Rev. B 2011, 84, 085313. [Google Scholar] [CrossRef]
  25. Gnani, E.; Reggiani, S.; Gnudi, A.; Parruccini, P.; Colle, R.; Rudan, M.; Baccarani, G. Band-Structure Effects in Ultrascaled Silicon Nanowires. IEEE Trans. Electron Devices 2007, 54, 2243–2254. [Google Scholar] [CrossRef]
  26. Jin, S.; Fischetti, M.V.; Tang, T.-W. Modeling of electron mobility in gated silicon nanowires at room temperature: Surface roughness scattering, dielectric screening, and band nonparabolicity. J. Appl. Phys. 2007, 102, 083715. [Google Scholar] [CrossRef]
  27. Wang, J.; Rahman, A.; Ghosh, A.; Klimeck, G.; Lundstrom, M. On the Validity of the Parabolic Effective-Mass Approximation for the I-V Calculation of Silicon Nanowire Transistors. IEEE Trans. Electron Devices 2005, 52, 1589–1595. [Google Scholar] [CrossRef]
  28. Trellakis, A.; Galik, A.T.; Pacelli, A.; Ravaioli, U. Iteration scheme for the solution of the two-dimensional Schrödinger–Poisson equations in quantum structures. J. Appl. Phys. 1997, 81, 7880–7884. [Google Scholar] [CrossRef]
  29. Lundstrom, M. Fundamentals of Carrier Transport; Cambridge University Press: Cambridge, UK, 2000. [Google Scholar]
  30. Ferry, D.K.; Goodnick, S.M.; Bird, J. Transport in Nanostructures; Cambridge University Press: Cambridge, UK, 2009. [Google Scholar]
  31. Jin, S.; Fischetti, M.V.; Tang, T.-W. Modeling of surface-roughness scattering in ultrathin-body SOI MOSFETs. IEEE Trans. Electron Devices 2007, 54, 2191–2203. [Google Scholar]
  32. Muscato, O.; Wagner, W. A class of stochastic algorithms for the Wigner equation. SIAM J. Sci. Comput. 2016, 38, A1483–A1507. [Google Scholar] [CrossRef]
  33. Camiola, V.D.; Mascali, G.; Romano, V. Numerical simulation of a double-gate MOSFET with a sub-band model for semiconductors based on the maximum entropy principle. Contin. Mech. Therm. 2012, 24, 417–436. [Google Scholar] [CrossRef]
  34. Camiola, V.D.; Mascali, G.; Romano, V. Simulation of a double-gate MOSFET by a non-parabolic energy-transport sub-band model for semiconductors based on the maximum entropy principle. Math. Comput. Mod. 2013, 58, 321–343. [Google Scholar] [CrossRef]
  35. Muscato, O.; Pidatella, R.M.; Fischetti, M.V. Monte Carlo and hydrodynamic simulation of a one dimensional n+nn+ silicon diode. VLSI Des. 1998, 6, 247–250. [Google Scholar] [CrossRef]
  36. Anile, A.M.; Muscato, O.; Romano, V. Moment equations with Maximum Entropy closure for carrier transport in semiconductor devices: Validation in bulk silicon. VLSI Des. 2000, 10, 335–354. [Google Scholar] [CrossRef]
  37. Muscato, O.; Di Stefano, V. Local equilibrium and off-equilibrium thermoelectric effects in silicon semiconductors. J. Appl. Phys. 2011, 110, 093706. [Google Scholar] [CrossRef]
  38. Di Stefano, V.; Muscato, O. Seebeck Effect in Silicon Semiconductors. Acta Appl. Math. 2012, 122, 225–238. [Google Scholar] [CrossRef]
  39. Muscato, O.; Di Stefano, V. An Energy Transport Model describing heat generation and conduction in silicon semiconductors. J. Stat. Phys. 2011, 144, 171–197. [Google Scholar] [CrossRef]
  40. Muscato, O.; Di Stefano, V. Hydrodynamic modeling of the electro-thermal transport in silicon semiconductors. J. Phys. A Math. Theor. 2011, 44, 105501. [Google Scholar] [CrossRef]
  41. Muscato, O.; Di Stefano, V. Heat generation and transport in nanoscale semiconductor devices via Monte Carlo and hydrodynamic simulations. COMPEL 2011, 30, 519–537. [Google Scholar] [CrossRef]
  42. Muscato, O.; Di Stefano, V. Electro-thermal behaviour of a sub-micron silicon diode. Semicond. Sci. Technol. 2013, 28, 025021. [Google Scholar] [CrossRef]
  43. Muscato, O.; Wagner, W.; Di Stefano, V. Heat generation in silicon nanometric semiconductor devices. COMPEL 2014, 33, 1198–1207. [Google Scholar] [CrossRef]
  44. Muscato, O.; Di Stefano, V. Electrothermal transport in silicon carbide semiconductors via a hydrodynamic model. SIAM J. Appl. Math. 2015, 75, 1941–1964. [Google Scholar] [CrossRef]
  45. Mascali, G. A hydrodynamic model for silicon semiconductors including crystal heating. Eur. J. Appl. Math. 2015, 26, 477–496. [Google Scholar] [CrossRef]
Figure 1. Schematic view of a SiNW. Electron transport is assumed to be one-dimensional in the x-direction.
Figure 1. Schematic view of a SiNW. Electron transport is assumed to be one-dimensional in the x-direction.
Entropy 18 00368 g001
Figure 2. SiNW band structure.
Figure 2. SiNW band structure.
Entropy 18 00368 g002
Figure 3. Cross-sections of a gate-all-around SiNW transistor.
Figure 3. Cross-sections of a gate-all-around SiNW transistor.
Entropy 18 00368 g003
Figure 4. Charge density and total potential along the cross-section at x = 48 nm and z = 0 nm, under a 1 V gate bias.
Figure 4. Charge density and total potential along the cross-section at x = 48 nm and z = 0 nm, under a 1 V gate bias.
Entropy 18 00368 g004
Figure 5. Linear density for the A and B valleys versus the simulation time.
Figure 5. Linear density for the A and B valleys versus the simulation time.
Entropy 18 00368 g005
Figure 6. The mean velocity (40) versus the simulation time, obtained with and without the surface roughness scattering mechanism.
Figure 6. The mean velocity (40) versus the simulation time, obtained with and without the surface roughness scattering mechanism.
Entropy 18 00368 g006
Figure 7. The mean energy (41) versus the simulation time, obtained with and without the surface roughness scattering mechanism.
Figure 7. The mean energy (41) versus the simulation time, obtained with and without the surface roughness scattering mechanism.
Entropy 18 00368 g007
Figure 8. The mean energy flux (42) versus the simulation time, obtained with and without the surface roughness scattering mechanism.
Figure 8. The mean energy flux (42) versus the simulation time, obtained with and without the surface roughness scattering mechanism.
Entropy 18 00368 g008
Table 1. Silicon nanowire constants.
Table 1. Silicon nanowire constants.
SymbolPhysical ConstantValue
m e electron rest mass9.1095 × 10 28 g
m A * effective mass A = Δ 4 valley [21]0.27 m e
m B * effective mass B = Δ 2 valley [21]0.94 m e
T L lattice temperature300 K
ρ mass density2.33 g/ cm 3
v s average sound speed9 × 10 5 cm/s
D a c acoustic-phonon deformation potential9 eV
D o intra-valley deformation potential g-scat [6]1.1 × 10 9 eV/cm
ω o intra-valley phonon energy [6]63.3 meV
Z o number equivalent valleys [6]1
D i v inter-valley deformation potential f-scat [6]2 × 10 8 eV/cm
ω i v inter-valley phonon energy [6]47.48 meV
Z i v number equivalent valleys [6]2
ε A 0 A valley energy minimum [21]0
ε B 0 B valley energy minimum [21]117 meV
Δ s r rms height [6]0.3 nm
λ s r correlation length [6]1.5 nm
Table 2. Low-field mobility (in c m 2 /Vs). SRS, surface roughness scattering.
Table 2. Low-field mobility (in c m 2 /Vs). SRS, surface roughness scattering.
μ A μ B μ low
no SRS55662430
with SRS13142108

Share and Cite

MDPI and ACS Style

Muscato, O.; Castiglione, T. A Hydrodynamic Model for Silicon Nanowires Based on the Maximum Entropy Principle. Entropy 2016, 18, 368. https://doi.org/10.3390/e18100368

AMA Style

Muscato O, Castiglione T. A Hydrodynamic Model for Silicon Nanowires Based on the Maximum Entropy Principle. Entropy. 2016; 18(10):368. https://doi.org/10.3390/e18100368

Chicago/Turabian Style

Muscato, Orazio, and Tina Castiglione. 2016. "A Hydrodynamic Model for Silicon Nanowires Based on the Maximum Entropy Principle" Entropy 18, no. 10: 368. https://doi.org/10.3390/e18100368

Note that from the first issue of 2016, this journal uses article numbers instead of page numbers. See further details here.

Article Metrics

Back to TopTop