1. Introduction

For a gas of charge carriers of q-type, where q stands for (n) for electrons or (p) for holes, the semiclassical Boltzmann transport equation (BTE) may be written as follows, when charge carriers are not strongly correlated [1]

(1.a)

where f_q = f_q(k, x, t) is the carrier distribution function of the q-type gas in the phase space, u_q is the group velocity of carriers, F = ± e(ζ + u_q x B₎ is the applied force due to electric field ζ and magnetic field B and (∂f_q/∂t)_col is the rate of change of the distribution function due to the internal collisions between the charge carriers and the crystal lattice (phonons and other defects).

Since the validity of the semiclassical transport approaches for simulating extremely scaled devices (sub ) and ultra fast devices (working up to the terahertz range) has been lately a subject of debate, we just remind here the main assumptions on which the BTE in semiconductors is based [2].

● Collisions are localized in space, so that the mean free path of charge carriers is much greater than the de-Brouglie wavelength,

● Collisions are instantaneous in time such that the collision time is much smaller than the time between subsequent collisions,

● External electric field force does not vary greatly over small distances in the order of wave packet length,

● Particles are weakly correlated (Bogoliubov assumption [3]) and the many body effect is negligible.

Mathematically speaking, the BTE is valid as long as the following Barker condition is realized [4]:

(1.b)

where is the maximum applied electric field, m* is the carrier effective mass, is collision time and is mean free time between collisions and q = k–k’ is the wave vector of (emitted or absorbed) phonons during collisions.

For strongly ionized gasses and degenerate semiconductors, where particle-particle collisions are not negligible, the Boltzmann transport equation may be replaced by the Fokker-Planck equation [5-6]. For ultra small devices, whose active layers dimensions (at least one of them) are in the order of De Broglie wavelength, the BTE should be replaced by one of the quantum transport formulations [7-11], like the Wigner-Boltzmann transport equation (WBTE) [9-11]. However, there exist several approaches to include the degeneracy effects and quantum corrections to the semiclassical transport models (e.g., [12-13]).

The solution of semiclassical BTE is carried out by a variety of methods [14-26], among which the iterative methods [15], the matrix methods [16], the Monte Carlo method [17-21], the cellular automata method [22], truncated expansion (usually in Legendre polynomials and spherical harmonics) methods [23-26]. In addition, the so-called energy transport model (ETM) and moment methods [27], which are sometimes called the hydrodynamic model (HDM), have been also used to solve the BTE in the hot carrier regime. In the ETM approach the distribution function is decomposed into an even and odd parts and the collision term is expressed using a microscopic relaxation time approximation [28].

The hydrodynamic description of the gas of charge carrier in semiconductors consists in the study of the evolution of certain macroscopic quantities, whose physical meaning is significant and whose quantity is eventually measurable. So, the HDM is constructed by multiplying the semiclassical BTE with a number of weight functions and integrating these equations over the phase or momentum space, so that a set of differential equations in physical space and time is obtained. This method is hence an averaging process that leads to the loss of some microscopic information in the distribution function. However, in many practical cases the resultant average (or macroscopic) parameters retained by the equations in the physical space are sufficient to capture the essential features of the carrier transport phenomena. The HDM has been widely used in the study of aerodynamic flow, and was suggested to study the hot electron transport in multivalley semiconductor devices by Bløtekjaer [29]. Since then, several contributions have been introduced to improve the hot carrier transport models in semiconductors and the HDM in particular [30-47]. For instance, the nonparabolicity of energy bands has been included in the HDM for the first time by Thoma et al [38]. The energy transport model ETM was reformulated by Chen, Kan and Ravaioli, who assumed non-Maxwellian and nonparabolicity correction factors, to solve the contradiction involved in evaluating a non-zero heat flux term for a Maxwellian distribution [39]. Also, the number of considered moments in the HDM has been increased in order to improve the HDM accuracy [44]. Unfortunately, all the added improvements to the HDM were associated with the addition of a huge number of new physical parameters [45-46]. On the other hand, some authors neglected the role of some important parameters in order to reduce the complexity of the model. This led to several paradoxal results and the failure of the HDM, in certain situations [47]. We try here to generalize the previously introduced models while reducing the number of physical parameters, by correlating them all to the energy relaxation time. We start by introducing the system of hydrodynamic equations in section II. In section III we discuss the system closure hypotheses. We present the physical definition of the carrier temperature. In section IV we present an almost complete description of the physical parameters in the bulk of semiconductor devices. Finally we present our conclusions in section VIII.

2. System of Hydrodynamic Equations

The multiplying weight functions in the HDM for semiconductors are usually chosen as powers of increasing order of the carrier vector k_q, or carrier group velocity u_q, or even the carrier momentum p_q. The powers of these parameters are usually formatted using some appropriate scaling factors to get physically meaningful quantities. For instance, if we adopt the powers of the carrier group velocity as weight functions (e.g., 1, u_q and u_q²), then we use the carrier energy E_q = ½m_q^*u_q², instead of u_q² as the second-order multiplier. So, the set of conservation equations of carrier concentration, momentum and energy may be obtained by multiplying the BTE by the first powers of the carrier group velocity u_q (e.g., 1, u_q and E_q), and then integrating the two sides of the resulting equations over the whole k-space (see Appendix I). For the sake of simplicity and brevity we write here only the electron equations in a single-valley semiconductor. In multi-valley semiconductors, we may assume a single lumped valley or write similar equations for each valley of the semiconductor [37]. The electron density-, momentum- and energy- conservation equations, in a semiconductor, can be put in the following form:

(2.a)

(2.b)

(2.c)

where n, v_n, P_n and ω_n are the electron concentration, mean velocity, electrokinetic pressure tensor and mean energy, respectively. Also, S_n is the electron energy flow and m_n^-1 is the mean inverse mass of electrons. We dully note that we don’t adopt the effective mass approximation here, and the second derivative of electron energy just appears as a result of the averaging process in the k-space (Appendix I). Therefore, the mean inverse mass of carriers is generally a tensor whose components m_n^-1_(i,j) = (1/ħ²) <∂²E_n/∂k_i∂k_j> are supposed functions of the mean carrier energy. When the energy band structure of the semiconductor material is assumed parabolic, m_n^-1 is a constant and equal to the conventional electron inverse effective masse (m_n*)^-1 at the bottom of conduction band. However, in this work we consider m_n as a scalar quantity averaged over the semiconductor valleys At the bottom of conduction band the average mass will be denoted by m_no and its value in Si , according to optical measurements is about 0.275 m_o. The subscript ‘col.’ In the above equations means the rate of change of the quantity between brackets due to electron gas collisions with the semiconductor crystal lattice and its defects. The sign stands for the tensorial product of two terms such that:

(3.a)

The electrokinetic pressure tensor in the above equations is defined as follows:

(3.b)

We discuss the significance of this definition in the following section. The carrier energy flow (or energy flux) is given by the following expression, for electrons:

(4.a)

where Q_n is called the electron heat flux. In the effective mass approximation, where we express the electron energy by E_n= ½ m^*_n u_n², the electron heat flux vector Q_n is defined as the third central moment of the electron distribution function f_n:

(4.b)

So, the heat flux component Q_n(i) is the flux of random thermal energy across a surface with its normal oriented in the i^th direction. However, as the electron distribution function is not known we don’t make use of this equation in our model. We dully note that the formal definition of the gas pressure tensor P_n,(ij) is the flux of the i^th component of the momentum across a surface element moving with the mean velocity oriented with its normal in the j^th direction [48], such that:

(5.a)

where is the mass density of the gas. In the ideal gaz theory, the gaz mass density is a macroscopic quantity and is taken outside the averaging brackets. So, the gas mass density is given by:

(5.b)

such that P_n ≡ P_n. However, in clasical transport models, the electron gas mass density is usually related to the electron effective mass such that Hence, our definition of the gas pressure, P_n, is different from the classic definition of electron gas pressure tensor, P_n:

(5.c)

Therefore, the electrokinetic pressure tensor P_n may be related to the classic definition P_n via a correction factor, g_n, such that:

(5.d)

We call this factor, the electron gas pressure g-factor. This factor can be calculated from MC simulation in the bulk of semiconductor with arbitrary band structure.

In absence of magnetic field, the applied force on electrons is given by F = -eζ, where the electric field ζis equal to the negative gradient of the electrostatic potential . In heterogeneous devices, the applied force on electrons is given by F= -∇E_c, where is the edge of the conduction band. Similarly, the applied force on holes is F=∇E_v, where and E_vo being the edge of the valence band. The electrostatic potential distribution can be found by solving the Poisson equation:

(6)

where p is the hole concentration, εis the dielectric constant, P is the polarization vector, N_t is the density of traps and Dop = |N_d⁺-N_a^-| is the net ionized concentration in the semiconductor.

3. Closing the System of Hydrodynamic Equations

In order to close the system of conservation equations for charge carrier density, average momentum and average energy (2-a,b,c) we usually adopt some assumptions to approximate the terms which depend on the knowledge of the distribution function. The main assumptions are usually concerned with the following items:

1. definition of the carrier gas pressure and temperature,

2. approximation of collision terms,

3. modeling the heat flux term (the third center moment of the distribution function), and

4. modeling any other higher moment terms, if the number of moments is greater than 3.

When the kinetic pressure of the gas of electrons is considered as a scalar quantity, one can define the electron gas temperature, T_n, by the law of perfect gases P_n = n k_B This assumption does not imply a spherically-symmetric distribution function about its mean value in the momentum space (like the displaced Maxwellian distribution). Equation (5-a) can be merged with the hypothesis of the perfect gas law, in one relation defining the electron gas temperature as follows:

(7)

This coincides with the physical definition of the temperature as a measure of the average kinetic energy of the gas particles [49]. Using the above definitions of the electron gas pressure and temperature, we can easily prove that the mean electron energy ω_n is related to the scalar electron temperature T_n and scalar average mass m_n as follows:

(8)

Although we make use of the expression E_n=½ m_n* u_n² in order to derive the above relation between the eventual error due the effective mass approximation is absorbed in the gas pressure correction factor, which can be calculated by MC simulation. However, it is well known that the convection part of the the average electron energy is generally small compared the thermal part, except for device regions where the velocity overshoot is pronounced.

As for the collision terms, we make use of the macroscopic relaxation time approximation [31]. In this approximation, the collision term in the momentum conservation equation is expressed in terms of a macroscopic momentum relaxation time as follows:

(9.a)

Also, the collision term in the energy conservation equation can be expressed in terms of a macroscopic energy relaxation time as follows:

(9.b)

where is the electron energy at thermal equilibrium and T_L is the lattice temperature. Note that the collision term in the above two equations is decomposed into two terms. So, the momentum and energy relaxation times, include the effect of collisions over optical while the effect of impact ionization collisions is included in the [∂n/∂t]_col. The collision term [∂n/∂t]_col., which appears explicitly in the carrier density conservation equation can be expressed as follows:

(9.c)

where (U = R - G) is the net recombination-generation rate in the semiconductor. Note that the recombination rate, R, regroups the relevant recombination mechanisms in the semiconductor device (e.g., radiative, Shockley-Reed-Hall and Auger recombination mechanisms [50]). Also, the generation rate, G, regroups the applicable generation mechanisms in the semiconductor device (e.g., thermal and optical generation as well as the generation by impact ionization mechanism).

Using the relaxation-time approximation and our definition of the carrier temperature, the set of hydrodynamic equations (HDE’s) can be formulated as follows for a single valley semiconductor:

(10.a)

(10.b)

(10.c)

The electron drift mobility is defined as follows:

(11.a)

Similar equations can be written for holes. We dully note that the average electron energy is measured above the conduction band edge Ec, while the average hole energy is measured under the valence band edge Ev.

In steady state (up to the terahertz frequency range), the carrier momentum conservation equations can be further simplified by suppressing the partial time-derivatives so that:

(11.b)

(11.c)

The above conservation equations are similar to the current equations of the isothermal DDM, except for the acceleration and thermal diffusion terms. Also the carrier energy flow (or energy flux) can be put in the following forms, for electrons and holes:

(11.d)

(11.e)

As indicated by (4-b), the carrier heat flux terms is defined in terms of the carrier distribution function, which is unknown. Therefore, the carrier-heat flux is usually modeled by the Fourier relation [51]:

(12.a)

where k_n is the electron thermal conductivity. The carrier thermal conductivity is usually modeled by the Wiedemann-Franz law, which states:

(12.b)

where is the electrical conductivity and is the Lorenz number (for electrons). In the original Wiedeman-Franz formulation (using the microscopic relaxation time approximation) is given by:

(12.c)

where r is an exponent in the microscopic relaxation time expression, which is given by:

(12.d)

A more rigorous approach for modelling the carrier heat flux term consists in taking the third moment of the BTE into account. This procedure increases the number of differential equations to be solved and the number of parameters of the HDM. However, we show in section 5.4 how to derive a more accurate model of the carrier heat flux term and the Lorentz number, on the basis of the third moment of the BTE. We do this without increasing the number of differential equations of the HDM.

4. About the Definition of the Carrier Temperature

In this section we compare the different approaches for modeling the carrier temperature. Formally, temperature is a property which governs the transfer of thermal energy, or heat, between one system and another. The electronic temperature concept was initially introduced in 1947 by Fröhlich [52], and since then it has been exploited by several authors to model the high-field transport in semiconductors. In the isothermal drift-diffusion model (DDM), the carriers temperature is assumed equal to the lattice temperature T_L. However, the so-called static electron temperature T_n may be roughly estimated from the electric field distribution in several ways [53]. For instance, the static temperature may be defined as follows:

(13)

where is the electron drift mobility, v_n^sat is the electron saturation velocity andis the low-field drift mobility. Also, ζ^// is the parallel component of the electric field to the electron current density such that ζ ^// = ζ .J_n/|J_n|. The above equation is based on the assumption that the parallel component of the diffusion coefficient is almost constant and independt of the electric field [33]. Such static forms of the carrier temperature exclude the lag effect between the electric field and the carrier temperature, and hence the velocity overshoot phenomenon will not be observed. However, the definition of the carrier temperature in the simulation is given by:

(14.a)

Also, the carrier temperature in the conventional hydrodynamic models, where the effective mass approximation is assumed and the anisotropy of the distribution function is ignored, is defined by:

(14.b)

It is well known, that the above two definitions result in different temperatures [38]. In fact, the conventional hydrodynamic definition (14-b) is implicitly based on the displaced Maxwellian distribution at the electronic temperature.

(15.a)

where A is a constant. Assuming that the collisions are fully randomized, such that the drift energy may be neglected with respect to the thermal energy, then the electron distribution may be approximated to the so-called heated Maxwellian as follows:

(15.b)

where E_n = ½ m_n^*u_n² is the electron energy and B is constant. According to Staratton [28], the electronic temperature has no physical meaning, unless the distribution function is a displaced Maxwellian as described by equation (15). However, the Maxwellian distribution is only valid near equilibrium, and in highly doped regions, where carrier-carrier collisions are dominant. In fact, it has been confirmed by several authors (e.g., [21]) that the energy distribution function is not Maxwellian at high electron energies. In the so-called improved energy transport model [39], the nonparabolicity of the band structure is taken into account and the distribution function f_n(E) is appended by a non-Maxwellian factor such that:

(16.a)

with f_m(E) is the heated Maxwellian distribution function at an elevated temperature T_n. Then, a modified (or equivalent) carrier temperature T_m is defined as follows:

(16.b)

However, referring to our system of HDE’s (8), we make use of the formal definition for the electron temperature tensor:

(17.a)

When the carrier temperature tensor is reduced to a scalar quantity (such that T_n = T_n I where I is the unity tensor), then the above definition may be simply written as follows:

(17.b)

In this case the g-factor may be calculated from the following relation:

(17.c)

Figure 1. Electron temperature according to MC [29], classic HDM and our model

Figure 1 depicts the variation of T_n according to MC (with a nonparabolicity factor 0.5 eV^-1) and the classic hydrodynamic model (14-b) and to our definition (17-b). Note that assuming g_n ≈1 brings an error less than 3% in T_n, for electric fields up to 300 kV/cm. This error is less drastic than the error encountered when neglecting the convective part in the calculation of the electron energy. In fact, most of the commercial simulators neglect totally the convection part of the carrier mean energy. If we adopt the assumption of unity g_n, we can drop the g_n factor from the above system of hydrodynamic equations while keeping the non-parabolicity effects, via the average mass term. In this case, the scalar temperature definition reduces to:

(17.d)

This approximation means that we assume the electrokinetic pressure tensor P_n = nk_BT_n, and may be considered as a closure condition for the set of HDEs. However, if later full-band MC (FBMC) showed a significant deviation of the g-factor from unity at very high fields, it could be easily retained in the proposed HDM to get more accurate results.

5. Energy-Dependent Physical Parameters

The accurate modelling of transport parameters, as the carrier drift mobility and impact ionization rate at high energies, is vital in any transport model. In the DDM, these parameters were usually expressed as functions of the local electric field [54-55]. Nevertheless, we know that these parameters express physical mechanisms involving energy exchange and should be expressed as functions of the carrier energy, rather than the local electric field. Fortunately, the solution of the HDM produces such information about the carrier energy distribution and permits to express these parameters in a physical manner, along the semiconductor device. In the following subsections, we discuss the modelling of the energy-dependent physical parameters in the bulk of silicon devices.

5.1. Carrier Drift Mobility

5.1.1. Carrier Drift Mobility in the Bulk:

In the drift-diffusion model, the high-field drift mobility in the bulk of a semiconductor is usually expressed by the Caughey-Thomas model, in its scaled form [56]:

(18.a)

where is the low-field mobility of electrons, which is related to doping, temperature and injection level, and may be expressed by a suitable model, e.g., Dorkel-Leturcq [57] or Klassen [58] models. Also v_n^sat is the electron saturation velocity and may be expressed by the model of Canali et al. [59]:

(18.b)

However, in the hydrodynamic model where the carrier average energy can be determined, the high-field mobility can be expressed in terms of the carrier mean energy The energy-dependent drift mobility relation may be written in the following form [60]:

(18.c)

This model resembles the model of Hänsch and Mattausch [61], but we consider an energy-dependent relaxation time and we make no a priori assumption about the distribution function in its derivation [62]. Figure 2(a) shows the variation of electron drift mobility as a function of electron temperature, according to our model and according to the model of Hänsch-Mattausch. In our model is taken from MC simulation, as shown in Figure 4(a). In order to compare our model with the field-dependent phenomenological relations, one can substitute equation (18-b) into the set of HDE’s in steady state homogeneous case to obtain:

(19.a)

The solution of this quadratic equation (in μ_n ) is given by the following relation:

(19.b)

where we replaced ζ by the parallel field component ζ_//, which is responsible of carrier heating, to satisfy the Thornber scaling theory [63]. As shown in Figure 2(b), this model coincides with the conventional field-dependent Caughey-Thomas mobility model [56] as well as the measured drift mobility, according to Canali et al [59]. The change of saturation velocity with lattice temperature may be also modelled as follows:

Figure 2(a). Electron drift mobility in Si at 300K versus electron temperature, according to the model of Hänsch-Mattausch with different values of [43] and according to our model

Figure 2(b). Electron drift mobility in Si at 300K, according to Canali et al [56], Caughey-Thomas model [52] and our model

5.1.2. Carrier Drift Mobility near the Surface

The drift mobility near the surface of a semiconductor is greatly reduced because of the additional scattering mechanisms due to the rupture of crystal periodicity, the surface phonons, the surface roughness and the eventual presence interface states. In addition, the presence of confined electrons into a thin layer near the surface of MOSFETs (the inversion layer) gives rise to quantization effects (subbands). As suggested by Sabins and Clemens [64], the reduction in the effective surface mobility may be expressed as a unique function of the effective vertical field ζ_┴. The most famous model of the surface mobility in MOSFET devices has been, for long time, the Yamaguchi model [65]:

(20)

where ζ_┴= |ζxJn|/|Jn|, ζ_c is a constant (1.49x10⁴ V/cm for electrons and 1.87x10⁴ V/cm for holes in Si), and is the bulk drift mobility. According to Hänsch-Mattausch model [61], the drift mobility along the semiconductor surface is given by the following relation:

(21.a)

which resembles equation (19-b) but with the low-field mobility replaced with the reduced surface mobility which is given by:

(21.b)

where are constants. In fact, the channel mobility can be modeled by the Schwarz empirical formula, when we neglect the effect of generated surface states [66]:

(21.c)

Another model, due to Schwarz and Russek [67], considers the channel broadening quantum effects and expresses the channel drift mobility as follows:

(21.d)

where are constants. The effective vertical field ζ_┴ is sometimes related to the surface potential and may be approximated by the following relation, for MOSFETs wih n+ polysilicon gate under strong inversion conditions [68]:

(22)

where are the dielectric constants of silicon and gate oxide, and x_ox are the MOSFET gate voltage, threshold voltage and gate oxide thickness, respectively.

A more elaborate and universal field-dependent mobility model, which reproduces the measured mobility data of Takagi et al. [69] at the surface of semiconductor devices is given by Darwish et al. [70]. According to Darwish’s model, the total carrier mobility is expressed by an equation similar to the scaled Hansch-Murrata model (21.a), with is given by:

(23.a)

Here is again the bulk mobility, is the carrier mobility limited by surface roughness and is the carrier mobility limited by surface acoustic phonons. In the later model, the bulk mobility, is expressed by the klassen model [58], the surface acoustic, , follows the Lombardi model [71]:

(23.b)

where B and C are constants (B=3.61x10⁷ cm/s for electrons and 1.51x10⁷cm/s for holes, C=1.7x10⁴ cm^14/3/V^2/3s.K for electrons and 4.18x10³ cm^14/3/V^2/3s.K for holes) and the surface roughness limited mobility, is given by:

(23.c)

where are constants, depending on the technology for electrons and 4.1x10¹⁵ V/cm for holes,

5.2. Impact Ionization Rate

The generation by impact ionization is one of the most important mechanisms limiting the performance of semiconductor devices in general and the submicron MOSFETs in particular. In fact, the electron-hole pair generation by impact ionization, G_ii , influences so many characteristics of semiconductor devices such as the breakdown voltage and hot carrier injection currents. The impact ionization mechanism was identified, about 50 years ago, by McKay [72]. Since then, the literature of the impact ionization mechanism has been extensive, e.g. [73-108], but only a few reviews are available. The generation rate by impact ionization G_ii has been always expressed in terms of Townsend coefficients (or ionization rates of electrons and holes), as follows:

(24.a)

In the DDM, the impact ionization rates are conventionally expressed as functions of the local electric field. For instance, Wolff [73] used the first two terms of the Legendre expansion of the momentum distribution function, and derived the following electron impact ionization rate:

(24.b)

where and C are constants. According to Baraff [74], this procedure is only justifiable at very high electric fields where electron-phonon collisions are frequent and the momentum distribution function is isotropic. Also, Keldish [75] solved analytically the BTE at high energy assuming parabolic band structure and calculated the impact ionization rate using the symmetric part of the distribution function. The Keldish solution attracted the attention because it tends to the Wolff relation (exp[-C/ζ²]) at high fields and tends to a simple exponential exp[-C’/|ζ|] (where C’ is another constant) at moderate fields, which coincides with the experiment. In fact, the measurements of Chynoweth [74], Lee et al [77], Van Overstraten and Deman [78], Kotani and Kawazu [79], Maes [80] and Takayangi [81] showed that the impact ionization rate at moderate fields follow a simple exponential function. According to Grant’s model [82], the electron impact ionization rate follows the following exponential relation:

(24.c)

where ζ_cn is the critical field at the onset of impact ionization. Shockley tried to find out an explanation of this experimental exponential relation in his study of p-n junctions breakdown [83]. According to Shockley, the impact ionization mechanism is mainly due to streaming of “lucky” electrons which can escape from collisions, and hence the ionization rate of such electrons can be expressed as follows:

(24.d)

where d = (E_n^th/eζ) is the minimum path an electron should travel to gain the threshold energy and is the mean free path between collisions with optical phonons. This model has been utilized to describe the breakdown phenomenon as well as the hot carrier injection and leakage currents in MOSFET devices in DDM-based simulation for long time [84-85]. Unfortunately, the parameters of the lucky-electron exponential model are dispersed among the different measurements of different authors in the literature [84]. In fact, the impact ionization phenomenon happens when drifting electrons have enough energy to ionize the valence electrons of lattice atoms. As the electron energy depends on the whole path of the electron under the effect of both the electric field and collisions, the local electric field value cannot directly reflect the real value of its energy. Thus, the local electric field dependent models of the impact ionization are not physically adequate, and their success in certain classic cases is due to their adjustable parameters. However, with the advent of the hydrodynamic model for semiconductors, where the carrier energy distribution is available, several authors tried to express the impact ionization rates in terms of the carrier mean energy or temperature.

5.2.1. Schöll-Quade Impact Ionization Model

The model of Schöll and Quade [87-88] for imact ionization rate is based on the heated Maxwellian approximation for the symmetric part of the distribution function. According to this model, the impact ionization rate of electrons is given by:

(25)

where is the electron ionization threshold energy, is a constant and erfc(x) is the complementary error function. Note that this model considers a parabolic energy band (constant effective mass m_n^*). In order to get adequate results at moderate fields, both E_n^th and are usually considered as adjustable parameters, with values which are extremely far from their physical meaning. For instance, Souissi et al [89] considered E_n^th ≈ 4eV and in order to fit this model with experiment.

5.2.2. Baccarani & Stork Impact Ionization Model

The model of Baccarani and Stork belongs to the category of semi-empirical models, which are based on the Shockley lucky-electron model. According to Baccarani and Stork [90], the rate of impact ionization is given by:

(26)

where are constants. According to Baccarani and Stork, the ionization threshold energy E_n^th =1.1eV, the energy relaxation length the mean free path between opical phono collisions and the adjustable parameter c is equal to 2x10¹² cm^-2. However, Souissi et al. [89] used the model of Baccarani and Stork in their investigation about impact ionization in bipolar transistors and found that the best fit with measured values implies E_n^th=5.077eV. It has been shown in the latter work that both the Schöll-Quade and Baccarani-Stork models may be used to predict the multiplication factor across the base-collector zone with reasonable error, only over a limited range of applied bias.

5.2.3. Non-Maxwellian Distribution Function–based Impact Ionization Models

In order to get rid of the disadvantages of the above mentioned models, some authors tried to calculate the impact ionization rate as a function of carrier energy, starting from its microscopic definition, on the basis of non-Maxwellian distribution functions [91-98]. For instance, Matsuzawa, Komahara and Wada [90] suggested the following form for the hot carrier tail of electrons distribution function:

(27.a)

where C₁ is a constant. Also, Sonoda et al. [96] and later Grasser et al [97] assumed the following two-components distribution, to account for both cold and hot carriers populations in nanoscale MOSFET devices:

(27.b)

where C₁, C₂, E_ref and b are adjustable parameters, which are used to tune the hot carrier distribution component, and the cold carrier (Maxwellian) component, in order to obtain reasonable results over a wide range of carrier energies. Therefore, when C₁=0, the distribution function becomes a heated maxwellian and when C₂=0 it tends to the tail distribution function. Using this four-parameters distribution function and assuming parabolic energy bands, Morris, Pass and Abebe introduced a simple formulation for the impact ionization rate on the basis of Gauss-Laguerre approximation [98].

It should be noted that several efforts have been devoted for the development a nonlocal impact ionization model on the basis of the hot-electron subpopulation, e.g., [99]-[100]. T-W.Tang and Nam have introduced a simplification of such hot-carrier subpopulation models [101]. The impact ionization coefficient is then fitted to homogeneous MC simulation in the bulk of Si as follows:

(28)

5.2.4. Advanced Impact Ionization Model

Unfortunately, the success of the above mentioned energy-dependent models of impact ionization is limited to specific cases at small and moderate fields. In addition, these models make use of a number of adjustable parameters, whose fitting values have sometimes irrelevant physical meaning. We think that one source of errors is their lake to a link with the realistic band structure of the semiconductor material. The following phenomenological model accounts for the energy band structure and can be used, with no need to curve fitting of its physical parameters [102].

(29.a)

where are the electron and hole average ionisation rates. Starting from the scaled form of the impact ionization rate, according to the Thornber theory [103]:

(29.b)

where are the critical fields at which the electrons surmount acoustic phonons, optical phonons and impact ionization scattering thresholds. Assuming that the ionization rate can be fully described by the kinetic moments one can write the energy conservation equation (10-c) in the following form in static field conditions:

(29.c)

where the first term represent the generation by impact ionization, the second term represents the loss of electron energy due to recombination ( being the electrons lifetime and n_o is the equilibrium electron density) and the third term represents the loss of electron energy due to collisions with optical phonons. The second term may be dropped from the above equation because. Also the the first and third terms can be combined, with an effective energy relaxation time, which includes the effect of impact ionization at high energies. In fact, this effect is not significant, according to Schöll and Quade [88]. Substituting the product eζ.v_n from (27-c) into (27-b), we get the following expression for the electron ionization rate at moderate and high fields:

(30.a)

where is the electron mean free path between optical phonon collisions and E_p is the mean optical phonon energy. At low and moderate fields, this model can be approximated to the following simple form:

(30.b)

where is the asymptotic value of impact ionization rate and is the mean ionization threshold energy for electrons. However, at high fields the impact ionization rate is proportional to Note that we consider and So, even at low fields the exponential realation does not imply a Maxwellian distribution because E_n^I is energy dependent. The energy dependence of the energy relaxation time, which implicitly includes the energy band structure effects, can be obtained from MC simulation or by measurement as indicated in the next section IV-C. Unlike classical models, the hard threshold energy (E_n^th) is replaced with a soft (energy-dependent) mean threshold energy (E_n^I) in our model. This coincides with Bude and Hess conclusions about the impact ionization soft threshold [102]. The calculation of as a function of the electron energy may be carried out by MC simulation in the bulk, using the relation [105]:

(31.a)

where is the total scattering rate by optical phonons. It may be also estimated from the so-called lucky-electron exponential model LE-EM [86]. In our analysis we used the following phenomenological relation, which is based on Crowell and Sze model [106], at the carrier temperature:

(31.b)

where is mean free path between collisions at low fields (when T_n ≈T_L). Figs. 3(a) and (b) depict the variation of in Si, according to the above relations. In our calculations, we took E_p = 60meV and and exponentially extrapolated the available MC data of as shown in figure 4(a). For the matter of comparison, we plotted according to Baccarani-Stork and Scholl-Quade models, beside our model. The Scholl-Quad model produces exceesively high rates when the physical value of is taken into account. So, we set equal to 12.6fs, about 3 orders of magnitude greater than its physical value [86]. Note that the value of in our scaled model is greater than all other models and continues to increase at high energies. In fact, the experimental values of may exceed 10⁶ cm^-1. For instance, the measurement of Kotani and Kawazu of is usually fitted with cm^-1 [79]. Also, the homogeneous MC simulation [105], indicates that the impact ionization pre-exponential is almost 10⁶ cm^-1 at high energies, as depicted by equation (26). It worth noting that the impact ionization rate at high energies (in the saturation region) is greatly influenced by the roll-off of the enrgy relaxation time, while E_n^th and have a relative influence at low and medium energies. This explains the failure of the previously suggested models, to properly estimate at high energies, even with adjusting E_n^th and the pre-exponential factors with non physical values, while fixing both For the matter of verification, Figure 14 depicts the I-V characteristics of a reverse biased p-i-n diode near avalanche, with different transport approaches (the DDM, the classic HDM with Scholl-Quade model and our proposed HDM model). As shown, the I-V characteristics obtained by our model are softer and the breakdown voltage is smaller and more realistic than the classic HDM results.

Figure 3(a). Electron mean free path between optical phonon collisions in Si at 300K, according to the lucky-electron exponential model (LE-EM), after [69], and according to the Crowel-Sze model [85] at the electronic temperature

Figure 3(b). Electron impact ionization rate in Si at 300K, as a function of average electron energy, according to Baccarani-Stork model Scholl-Quade model In all models

Figure 3 (c). Average ionization energy in Si at 300K, according to our model model

Figure 4(a). Variation of the electron energy relaxation time with electron average energy in Si at 300K, according to MC simulation. In the above curve, the energy bands are assumed parabolic. The lower curve takes the non-parabolicity of the conduction band into account (non-parabolicty factor 0.5V-1)

Figure 4 (b). Experimental value of electron and holes conductivity mass, after Rife [95]

5.2.5. Energy Broadening of Impact Ionization Threshold

The determination of the impact ionization threshold has been a subject of extensive research, both theoretically and experimentally [107]. In fact the measured values of the ionization threshold and those used in the theoretical models vary widely. As we pointed out in the last section, the mean ionization energy E_n^I in our model of impact ionization has a soft distribution around E_n^th. In order to evaluate the broadening of mean ionization energy, we define as follows:

(32)

Figure 3(c) depicts the ionization energy broadening as a function of mean electron energy. This broadening is mainly due to the variation of the energy relaxation time and the mean free path at high energies. Curiously, one may wonder if there is any reason to relate this phenomenon with the broadening in energy levels due to finite lifetime of energy states between scattering events at high energy [108-109]. In fact, it has been shown experimentally in [110] that there exists a line broadening transition, representing the transition from phonon scattering dominated to impact ionization dominated transport not only in silicon but also in any material with a bandgap. In addition, some theoretical studies of the impact ionization on the basis of quantum transport showed that a fixed impact-ionization threshold does not exist, and the impact-ionization scattering rate is drastically enhanced around the semiclassical threshold by the intracollisional field effect [111]. Of course, the broadening factor may be better evaluated by precise MC simulation, taking into account the impact ionization mechanism. However, with our simple model of and the exponentially extrapolated we found that is centralized around a minimum (almost zero) at , which is interestingly equal to 3/2 E_n^th. It should be noted here that we chose E_n^th= Eg ≈ 1.2 and chose _o.= 1000Å, such that the distribution of is as close as possible from the experimental LE-EM curve. As shown in figure 3(a), the two curves almost intersect at =1.2eV.

5.3. Carrier Energy Relaxation Time

The choice of the energy relaxation time is crucial in the HD simulation because it influences both the static and dynamic characteristics of the device to be simulated. According to the hydrodynamic theory, the exact modelling of the carrier energy relaxation time (and all other relaxation times) should be expressed as function of all the carrier moments like [112]. However, as we usually truncate the series of moments and satisfy ourselves with the few first moments, it is believed that much of the information about the energy relaxation time is related to the carrier mean energy and density [30]. According to Jacoboni and Reggiani [112], the energy relaxation time should decrease at high fields in a phenomenological way representing the tendency of the isotropic part of the distribution function to decay towards its equilibrium value. So, the later authors used the following equation, to evaluate the carrier energy relaxation time for warm electrons:

(33.a)

This relation, which is sometimes called the Seeger formula [114], reflects the energy continuity at steady state:

(33.b)

According to Seeger’s static formula, the energy relaxation time can be determined by substituting the static field relations from MC simulation [115]. Gonzalez et al. employed this relation to calculate the fitting parameters of their energy relaxation time model [116]. On the other hand, several authors considered the energy relaxation time as a constant when the electron distribution function is warm enough and the electron velocity tends to saturate. This is a good approximation when the band structure of the material is parabolic. For non-parabolic bands, the energy relaxation time continues to decrease when electron energy is further increased, as shown in Figure 4(a). In fact, it has been shown that the value of the energy relaxation time greatly influences the hydrodynamic simulation results in deep submicron transistors [117]. This is due the fact that the intervalley phonons becomes more effective in dissipating the electron energy at high energies, which lead to decreasing the electron energy relaxation time, at higher energies. For accurate modelling of the energy relaxation time, we propose the following phenomenological relation between the electron average effective mass m_n and for energies above the equilibrium value

(34)

where m_no is the density of states effective mass at the bottom of the conduction band (near equilibrium). This model is based on the physical fact that heavy electrons relax faster than lighter ones. In fact, the increased electron-phonon coupling at high fields (where the electron is said to be dressed in a cloud of phonons) leads to a modification of electron states and changes their energy (by addition of the self energy) and this increases the mass of electrons and makes their motion slower. Therefore, at high energies where the electron average mass increases as electron energy increases, the electron dissipates energy more efficiently and the energy relaxation time is shorter. According to the proposed relation, we only need the knowledge of the low-field value of the electron energy relaxation time and the average mass as a function of electron energy to obtain The low-field value of the average electron mass, m_no, is equal to the conductivity effective mass at the conduction band edge, and can be measured using the optically-detected cyclotron resonance (ODCR). The ODCR method is more accurate than the classical CR and has been recently used to measure the effective mass of several semiconductor compounds [118]. According to the optical measurements of Rife [119] the conductivity effective mass at 300K m_no =0.275m_o and m_no =0.365m_o as shown in figure 4(b). The ultrafast laser based pump-probe techniques have been used to measure the energy relaxation as well as the momentum relaxation times in Si. Using this technique, a time constant of 32±5 fs associated with momentum relaxation and an electron-phonon energy relaxation time of 260±30 fs have been extracted from the coherent-transient variations in (001) Si at 300K [120].

In order to verify our model, we used the MC simulation to calculate the average electron mass as a function of the average electron energy, and compared it with the measured ratio of (m_no/m_n), as shown in Figure 5. As shown in figure, the normalized value of the average inverse mass decreases as the electron energy increases. Then we calculated the energy relaxation time using the above phenomenological relation (above when electrons start to warm up). As shown in Figure 5, the discrepancy between the the energy relaxation time, according to MC and the calculated values according to our model lies within the measurement error (~10%) for a wide range of carrier temperatures (up to 2000K and may be extrapolated to higher values). Alternatively, our model can be used to calculate the energy dependence of the average electron mass from the energy dependence of the electron energy relaxation time The energy dependence of of several semiconductors (e.g., AlGaAs/GaAs, AlGaN/AIN/GaN) has been measured using the microwave-noise method [121]. This method is based on the measurement of the noise temperature, which is equal to the electron temperature T_n in the semiconductor, according to the following relation:

(35)

where we neglect the convection part of the electron energy and consider However, we have no experimental data yet about the energy relaxation time in Si (as a function of carrier temperature or energy) using this method. So, the proposed model which correlates to the already measurable should be very useful in this case.

Figure 5. Comparison between the normalized inverse mass (mno/mn), according to optical measurement [52] and our model. The model is calculated from the normalized energy relaxation time according to MC data, by our phenomenological relation (25)

Figure 6 (a). Heat flux beta factor in Si at 300K, calculated from MC simulation, with non-parabolicty factor 0.5V-1

Figure 6(b). Relative energy flow mobility in Si at 300K, calculated from MC simulation, with non-parabolicty factor 0.5V-1

5.4. Advanced Heat Flux Model

As we pointed out earlier, the heat flux term in the majority of previous HDM’s was either neglected or modelled by the Fourier relation (12-a). It has been shown that this model is not accurate for modelling the heat flux across semiconductor interfaces in general and p-n junctions in particular [122]-[124]. In order to increase the accuracy of the hydrodynamic model, we propose the following model, which is based on the third moment of the BTE in steady state:

(36.a)

where we suppressed the time derivative (∂S_n/∂t) and used the relaxation time approximation ( is the energy flux relaxation time). The above conservation equation (of energy flux), may be further simplified by approximating the average terms, as follows [123]:

(36.b)

where the correction factor _’ expresses the anisotropy of the distribution function and the correction factors expresses how much the distribution function is deviated from the Maxwellian distribution. These parameters can be estimated by MC simulation in the bulk of semiconductor. Then we can write the energy flux as the sum of a conduction part Q_n and a convection part (where J_n is a multiplier):

(36.c)

Here and the electron heat flux Q_n () is given by:

(37.a)

with

(37.b)

(37.c)

Note that the convection energy flux is similar to our 3-moments formulation (11), except for the correction factor and that is replaced with 5/2(k_BT_n). However, we just called the third moment equation in order to derive a more accurate formula of the heat flux term., which contains the effect of both the carrier temperature and the carrier density gradients.

5.5. Modified Wiedman-Franz Relation

If the anisotropy of the distribution function is negligible, we can substitute Then the heat flux term can be put in the following simple form:

(38.a)

Figures 6(a) and (b) depict the electron heat flux beta factor and the relative energy flux mobility as a function of electron energy. We calculated these parameters by MC simulation in the bulk of Si at 300Kwith non-parabolicity factor 0.5V^-1. Alternatively, they may be evaluated, as functions of carrier energy, and doping level, by appropriate physical models. For instance, the energy flux mobility may be modelled in much the same way as the drift mobility as follows [124]:

(38.b)

where is the energy flux mobility at low field and Also, the heat flux beta factor may be approximately expressed as follows in Si at 300K:

(38.c)

where are constants. For Si at 300K with doping levels Dop ≤10¹⁹cm^-3 and electron energy , we have: , .

When the carrier-concentration gradient is negligible (with respect to the temperature gradient), the heat flux formula (35-a) can be further reduced to the conventional Fourier relation (12-a), with the following carrier-thermal conductivity:

(39.a)

The above equation is equivalent to the conventional Wiedemann-Franz relation with the following ‘modified’ energy- and doping-dependent Lorenz number:

(39.b)

In the extreme case of very low doping and very low electric field, where as compared to (5/2+r) in the conventional model (12-c). Figure 13 depicts the heat flux as obtained by the classic Fourrier relation and the new model. As shown in figure, the consideration of the carrier gradients in the new model broaden and smooth the heat fluc profile near junctions.

5.6. Tunneling & Hot Carrier Injection Currents

The tunnelling and hot carriers injection currents paly an important role in modern semiconductor devices in general and the ultra short MOSFET devices in particular. There are two tunnelling (through the barrier) mechanisms and four hot carrier injection (over the barrier) mechanisms which are commonly encountered in MOSFET devices, namely:

● Field-assisted (Fowler-Nordheim or FN) tunnelling,

● direct tunnelling (DT),

● drain avalanche hot carrier injection (DAHCI),

● channel hot electron injection (CHEI),

● substrate hot electron injection (SHI) and

● secondary generated hot electron injection (SGHEI)

These mechanisms are usually monitored by measuring the gate current, I_G, and the substrate current, I_sub. In the following subsections we describe briefly these parameters and how they are modeled within the framework of the hydrodynamic model.

5.6.1. Tunneling Currents

The tunnelling currents are either direct or of field-assisted (Fowler-Nordheim) type. The direct tunnelling can take place at low gate voltages via thin oxides layers (less than 50Å) and has a weak dependence on the gate field. The FN tunnelling is a field-assisted mechanism which is strongly dependent on the gate field and dominates in modern MOS structures. The gate tunnelling current density is given by the following integration [125]:

(40)

where T(E_n) is the tunneling probability, g_c(E_n) is the density of states in the conduction band, is the electron group velocity perpendicular to the Si surface . The integration in the above equation starts from the semiconductor conduction band edge, as a reference. For trapezoidal barrier of upper height and lower hight (), where V_ox is the potential voltage across the insulator gate, the tunnelling probability is usually expressed using the WKB approximation as follows:

(41.a)

In the upper triangular part of the barrier, where , and

(41.b)

in the lower rectangular part of the barrier where . Here are the electric field and effective mass of electrons at the insulator (for SiO₂, m_ox =0.42m₀). In the case of the field-assisted Fowler-Nordheim mechanism, we consider the triangular part of the barrier, which is described by (38-a). Assuming the energy distribution function Maxwellian and the energy bands parabolic, we get the following field-dependent gate current density:

(42)

where the are constants, related to the energy barrier height between the insulator and the injecting conductor as follows:

(43.a)

(43.b)

The tunnelling mechanism can be initiated by low energy carriers (near equilibrium). However, it is also possible for hot electrons to tunnel when their energy is not sufficient to surmount the barrier. In order to obtain more realistic results of the gate tunnelling current in short MOSFET devices, some authors assumed a non-Maxwellian distribution function and considered the a nonparaboicity of energy bands [126]. The total gate current may be partitioned between the source and drain, as follows:

(44.a)

where

(44.b)

where W and L are the MOSFET channel width and effective length, respectively. When the substrate bias is high enough, the thermoionic emission across the gate oxide dominates the tunnelling current. The thermoionic emission current can be calculated by a formula similar to equation (40), with replacing the the tunnelling probability, T , by the probability that an electron will not scatter in the image potential well in the oxide region, P.

(45.a)

with

(45.b)

where x_ox is the gate oxide thickness and is the mean free path of electrons inside the gate oxide (for SiO₂, ). In order to account for both tunnelling and thermoionic emission, the above equation should be modified as follows [127]:

(45.c)

Finally it should be noted that the Fowler-Nordheim tunnelling gate current may photogenerate e-hole pairs at the substrate and hence it is likely to be the main origin of the substrate leakage current, rather than the anode hole injection [128].

5.6.2. Hot Carrier Injection Currents

The channel hot electron injection (CHEI) occurs when the gate voltage is high enough and V_G ≈V_D.  Channel carriers which are accelerated by lateral fields and transported from the source to the drain may be driven towards the gate oxide before they reach the drain. The drain avalanche hot carrier injection (DAHCI) occurs when a high drain voltage is applied in non-saturation mode (V_D>V_G) which results in a high electric field near the drain, and hence accelerating channel carriers into the drain depletion region. The accelerated channel carriers collide with Si atoms valence electrons, creating electron-hole pairs by impact ionization. Some of the generated electron-hole pairs are again accelerated and may acquire sufficient energy to surmount the Si/SiO₂ barrier. The hot carriers that surmount the gate oxide barrier can inject into the gate oxide where they may be trapped and gives rise to a shift in the MOSFET threshold voltage, surface mobility and dynamic conductance. The substrate hot electron injection (SHEI) occurs at high substrate bias, |V_B| >> 0 and strong inversion, usually with both drain and source grounded. Under this condition, carriers of one type in the substrate are driven by the substrate field toward the Si-SiO₂ interface. Carriers can be generated by external optical or thermal excitation. As carriers move toward the substrate-oxide interface, they gain further kinetic energy from the high field in surface depletion region. They eventually overcome the surface energy barrier and get injected into the gate oxide, where some of them are trapped. The secondary generated hot electron injection (SGHEI) occurs under conditions similar to DAHCI, when V_D>V_G, which lead to impact ionization of hot carriers. However, SGHEI involves secondary carriers that are created by an earlier incident of impact ionization and driven under the influence of the substrate bias. This bias produces a field that drives hot carriers toward the surface region, where they gain more energy to overcome the gate oxide barrier [124]. The so-called “channel-initiated secondary electron (CHISEL) is a variant of the SGHEI, which relies on ionization feedback and is activated by a negative substrate bias V_B<0. CHISEL is used as a reliable programming technology in recent low-voltage Flash EEPROM devices, with floating gate lengths down to 0.2 [131]. Figure 18 illustrates the different hot injection and tunnelling mechanisms in MOSFET devices.

The injection of hot-electrons into the gate-oxide has been usually modelled by the lucky-electron model [83]. According to Hu et al., the probability of a channel electron reaching the gate terminal as a combination of the probabilities of the following events:

● the electron gains sufficient energy in the lateral electric field to overcome the interfacial

● potential barrier and retains this energy after a collision directs its momentum towards the interface,

● the electron reaches the interface without suffering any more collisions, and

● the electron is not scattered back into the semiconductor in the image-force potential well present close the interface

The lucky-electron model, despite of its drastic assumptions, remains the most useful approach to estimate the injection currents with reasonable agreement with experimental results. On the basis of the lucky-electron model, it can be shown that the fraction of the channel carrier supply which is injected into the gate oxide is given by [129] :

(46.a)

where, is the interfacial potential barrier, is the electron mean free path and is the effective lateral electric field. It can be shown that the effective barrier potential is given by:

(46.b)

where a and b are constants, expressing the main barrier height and lowering effects due to image force and quantum tunnelling. For Si/SiO₂ interface, , a = and b = 2.59x10^-4 (Vcm)^½. As the substrate current is generated due to impact ionization in the lateral electric field close to the drain, it is linearly proportional to the drain current. On the basis of the lucky-electron concept, it can be shown that the substrate currents is given by :

(47)

Where I_ds is the drain current, L_d is the length of the pinch-off region and is the electron impact ionization coefficient in the velocity saturated part of the channel. So, the substrate current can be given by:

(48)

where, E_n^th is the impact ionization threshold energy for electrons. As seen from this equation, the ratio of the substrate and drain currents provides a direct measure of the average lateral electric field. In addition, this multiplication ratio, is often used as a monitor for the device degradation.

6. Lattice Heat Conservation Equation

Some recent studies demonstrated the important role of lattice heat and thermal activation of the low voltage impact ionization in deep-submicrometer MOSFETs [132]. In addition, the recent studies of silicon on insulator (SOI) devices by the classical HDM, showed certain anomalies due to lack of a heat evacuation mechanism at the insulator interface [133]. Some authors suggested to cure this problem by considering a tensorial temperature and modifying the closure condition for energy flux term [134]. Othors authors suggested coupling the lattice heat equation to the set of HD equations [135].

The classical heat diffusion equation, which resolves the temporal and spatial distribution of lattice temperature along the semiconductor, reads:

(49.a)

where is the lattice density, c_p is the lattice specific heat and their product is equal to the heat capacitance () of the semiconductor material. Also k^th is the thermal conductivity of the semiconductor and the heat source term H_s is given by [136]:

(49.b)

where is the absolute average electron energy and is the absolute average hole energy in the crystal lattice. In Si at 300K, gm/cm³ and c_p =7.03x10⁴ cm²/Ks² and their temperature dependence is negligible [137]. We dully note that the carrier thermal conductivity is a small fraction of the Si material thermal conductivity k^th because heat conduction in semiconductors is dominated by phonon transport. It should be also noted that for T_L > 100K, the thermal conductivity in Si is practically independent of doping concentration [ 138]. The temperature dependence of k^th in the bulk of crystalline Si may be modelled according to Glassbenner and Slack [139] as follows:

(49.c)

where and k_o ≈ 1.45 for Si. However, in thin semiconductor layers the thermal conductivity is reduced due to the phonon scattering process at layer boundaries (the size effect), as described in [140]. An analytical expression describing the thermal conductivity as a function of distance, into thin Si layers can be found in [141]. A more rigorous approach to solve the problem of heat evacuation in SOI devices, consists in solving the BTE for phonons. In fact it is well known that the LO phonons have small group velocities (~10⁵ cm/s) and anharmonically decay into faster acoustic phonons which in turn transport the energy out off the semiconductor device. So, the phonon BTE may be simplified by assuming a two-fluid phonon transport model involving a stationary optical reservoir mode and a propagating acoustic mode. The latter approach may consider not only the heat evacuation by optical phonons, but also their anharmonic decay to ‘fast propagating’ acoustic phonons. As a compromise between the above mentioned approaches, the a lattice heat continuity equation may be added to the set of HDE’s with phonon energy relaxation term. The transformed phonon BTE in the relaxation time approximation reads [142]

(50.a)

where is the phonon energy per unit volume, υ the phonon velocity and is the phonon energy relaxation time and H_s is the volumetric heat generation rate absorbed from hot electrons, which can be extracted from the solution of the HDE’s. According to the kinetic theory, the phonon energy relaxation time is related to the phonon velocity v_ph and mean-free path by the following relation [143]:

(50.b)

So, the lattice energy conservation equation may be written as follows:

(50.c)

where W_L = ∫C_v dT_L is the lattice energy density, W_o is the lattice energy density at equilibrium (at T_L=T_o) and S_L = -k^th∇ T_L is the lattice energy flux (or heat flux). We note that the overall lattice energy conservation dictates the following relation:

(51.a)

where is equal to the RHS of the energy conservation equation (10-c) for electrons and for holes. When the net generation-recombination rate is negligible, represents the local energy transfer rate from crystal lattice to the electron and hole gaz subsystems. Therefore we can write the following correlation function between the charge carriers relaxation times and phonons relaxation time:

(51.b)

Note that our formulation of the lattice heat conservation equation is tightly coupled with the set of HDEs, via the average electron and hole average energies (which makes part of the heat source term H_s) and the energy relaxation correlation function.

7. Quantum Corrections

The so-called quantum hydrodynamic model (QHDM) can be derived by a variety of methods [144-152], such as the first few moments of the WBTE or using the Madelung transform [147] or the principle of minimum entropy [148]. The formulation of the QHDM usually retains the lowest order term in the series expansion of a complicated quantum operator acting on the Wigner distribution function and the system potential energy. From the mathematical point of view, the most significant difference between the semiclassical HDM and the QHDM is the quantum potential term, V_q, which appears in the average energy expression [150]:

(52.a)

where is the average momentum of electrons and V_q is given by:

(52.b)

The quantum potential term is sometimes merged with the electron rature, to form the so-called quantum rature T_qn:

(52.c)

The quantum correction appears also in the current equation (which is a reduced form of the momentum conservation equation) in the form of an additional current J_q. The following electron current equation depicts this correction term [151]:

(53.a)

with

(53.b)

Note that the quantum correction term is significant when the electron concentration n varies so rapidly over small distances in the order of de-Broglie wavelength (h/p). Adding the quantum correction term in this manner is sometimes referred to as the density gradient (DG) method. The DG theory has been used to model the quantum confinement as well as the tunnelling currents in MOS devices [152].

8. Quasi-Ballistic Transport Considerations

It has been shown that the operation of sub 0.1 micron MOSFET transistors is quasi-ballistic, in a small portion of the channel [153-158]. Some authors shows that quasi-ballistic transport in nano-transistors takes place near the source region and others showed that quasi-ballistic transport and back scattering take also place near drain region as well [157]. In all cases, one should be aware that even when the transport is full ballistic, the contact resistance of the device is still there limiting the transport process. However, the quasi-ballistic transport can be also considered within the framework of the semiclassical theory [156]. In fact, it has been shown that semiclassical transport models along the MOSFET channel may remain reliable down to 10nm [157]. In addition, it has been shown in the latter article, that the semiclassical model is more reliable than the quantum simulation, even when the ratio of coherent transport along the MOSFET channel is not negligeable (down to 10nm). Figure 7 depicts the electron velocity distribution along a 50nm n-i-n structure, in the ballistic mode and using different transport approaches.

9. Sample of Simulation Results

Figure 1 depicts the electron temperature T_n, according to Monte Carlo simulation, the classic HDM and according to our definition. As shown in figure, the electronic temperature is reasonably estimated using our model when taking g_n≈1, for electric field up to 200 kV/cm. For the matter of demonstration of the proposed HDM and comparison with other DDM-based simulation, we simulated the carrier transport across a reverse-biased silicon p-i-n diode. The simulated diode is 100m long and has a p⁺ side with Gaussian distribution and peak doping of 10¹⁹cm^-3., as shown in Figure 8. The diode base n^--layer is homogenously doped with 10¹⁴cm^-3, and the metallurgical junction is almost at 11 apart from the p⁺ contact. Figures 9 and 10 depict the two-dimensional distributions of electron and hole density and temperatures inside the p-i-n diode at 75V reverse bias. Figures 11 and 12 depict the one-dimensional distributions of electron and hole density and temperature along the p-i-n diode at 500V reverse bias. Also figures 13 and 14 depict the one-dimensional distributions of electron and hole energy flux and heat flux along the p-i-n diode at 500V reverse bias. In addition. Figure 15 shows the I-V characteristics of the simulated diode near breakdown. We discussed in [159] the effect of the hydrodynamic simulation on the distribution of charge carriers and showed that discrepancy from the DDM-based results is due to the carrier thermal back diffusion. In order to show the important role of the energy relaxation time model in the HD simulation, we traced the electron distribution across the p⁺-n junction at various values of the Figure 16 depicts the electron density distribution along a reverse biased p-i-n diode at three values of the energy relaxation time. Note that the electron distribution, which is obtained using the HDM with lower tends to the DDM classical distribution. This means that for shorter energy relaxation times, the electrons are more efficient in dissipating their energy. For the limiting case when tends to zero, the electrons will tend to dissipate their gained energy at once and their average energy is equal the lattice energy This is exactly the case of the DD simulation, where we consider the electron temperature equal to the lattice temperature (T_n=T_L) at all fields.

Figure 7. Effect of quasi-ballistic transport on the electron velocity across an n-i-n diode, according to several models [118]

Figure 8. (a) Two-dimensional doping profile of the simulated p-i-n diode

Figure 8. (b) Lateral doping profile along the x-axis of the p-i-n diode

Figure 9. (a) Electron and (b) hole density distribution in the p-i-n diode (at –75V bias)

Figure 10. (a) Electron and (b) hole temperature distribution in the p-i-n diode (at -75V bias, 300K)

Figure 11. Electron density distribution along the p-i-n diode (at 500V reverse bias, 300K), on the basis of DDM and HDM

Figure 12. Electron and hole temperature distribution along a reverse biased p-i-n diode (at –500V), as calculated from the DDM electric field distribution (static) and by HD simulation

Figure 13. Electron heat flux distribution along the p-i-n diode at 500V reverse bias, 300K

Figure 14. Electron energy flux distribution along the p-i-n diode at 500V reverse bias, 300K

Figure 15. The I-V characteristics of the simulated the p-i-n diode near avalanche breakdown, as obtained by 2-dimmensional DD and HD simulation at 300K. The surface of the p-i-n diode is properly passivated and protected by field plates, so that the breakdown occurs at the junction curvature inside the p-i-n diode

Figure 16. Electron density distribution along a reverse biased p-n junction, with different values of energy relaxation time

Figure 17. Illustartion of different hot carrier injection mechanisms in MOSFET Devices

Figure 18. Illustration of hot carrier injection currents in MOSFET Devices

(54)

We may conclude from this discussion that the DDM simulation is a hypothetical case of the HDM where we assume and hence and this case means a virtual semiconductor whose collisions are so efficient that electrons relax to their equilibrium energy in almost no time.

(54)

The electrons in such a virtual semiconductor (with zero energy relaxation time) will not allowed to heat up even at high electric fields such that the mean electron energy will be almost equal to the lattice energy This result has been already reported for majority carriers in MOSFET devices [97], but has er been reported for minority carriers in bipolar devices, as far as the author knowledge.

10. Conclusions

In conjunction with sophisticated Monte Carlo methods that incorporate both realistic band structure and carrier scattering processes, the HDM is a basic tool to understand hot carrier behaviour in ever smaller and ultra fast semiconductor devices. We presented in this article a complete hydrodynamic model based on the first three moments of the semiclassical BTE for electrons, holes and phonons, with a few number of parameters. The model includes the band structure non-parabolicity effects and advanced energy-dependent models of the carrier mobility, impact ionization rate and heat flux. As the accurate estimation of the carrier temperature in semiconductor devices, influences the estimation of so many important parameters, such as the breakdown voltage and leakage currents, we proposed a simple definition of the carrier temperature, which takes the band structure into account. The definition results in temperatures which are closer to the simulation results than those obtained by the classical hydrodynamic models. We also proposed a phenomenological correlation between the electron average mass and energy relaxation time which cuts down the number of parameters of the HDM. In conclusion, our proposed system of HDE’s has a few energy-dependent parameters and they are all correlated to the energy relaxation time The later parameter can be obtained from the MC simulation [161] or by measurement. The proposed model has been applied to investigate the carrier transport in certain semiconductor devices and led to adequate results. From our simulation, we conclude that the undeliberated choice of the energy relaxation time value in the HD simulation leads to drastic results and may be to the failure of the HDM in certain cases, as reported in [47]. For instance, we show, by simulation, that reducing the energy relaxation time to unreasonable short values will mask the thermal diffusion mechanism and gives erroneous results, in both MOSFET and bipolar devices.

APPENDIX: Moments of the BTE

Consider the distribution function f(x,k,t), which obeys the BTE. Consider A(k) a microscopic quantity, which is a function of the wave vector k. Now, let’s multiply the BTE by A(k) and integrate the two sides of the equation over the entire k-space:

(55)

We have also, the mean value of , which we note by is given by:

(56)

Therefore, the 1^st term of the L.H.S of equation (55) becomes:

(57-a)

The 2^nd term of the L.H.S of equation (55) becomes:

(57-n)

Also, the 3^rd term of the L.H.S of equation (55) becomes:

(57-d)

where we effectuated an integration by parts on the k-variable. The index j of the force component F_j appears for any orthonormal basis (e.g. for Cartesian coordinates, j = x, y, z). The totally integrated term in equation (I-2c) disappears because . This means that all the moments of the BTE are finite. On the other hand, the second term of (I-2c) may be written as follows:

(57)

Here, we put , which is true if F is independent of k or when F is normal to k. Fortunately, this is the case for electric field force () and Lorentz’s force (). Now, by regrouping all (55) components we get:

(58)

By replacing A(k) by 1, and we obtain the 3 moment equations of the BTE. Table 1 summarizes the results.

Table 1. The integration of different terms of the LHS of the BTE, after being multiplied by A(k)