Subsections

• 2.1.1 The Drift-Diffusion Transport Model
• 2.1.2 Types of Partial Differential Equations
• 2.1.3 Systematic Derivation of the Macroscopic Transport Models
• 2.1.4 The Six Moments Transport Model
• 2.1.5 The Energy-Transport Model
• 2.1.6 The Drift-Diffusion Transport Model

2.1 The Analytical Problem

The analytical problem is basically formed by nonlinear partial differential equations. These equations are the Poisson equation (2.1) and the current continuity equations for the two carrier types in semiconducting materials, electrons and holes,

$$\ensuremath{\ensuremath{\ensuremath{\mathbf{\nabla}}}\cdot (}\var... ...on \ensuremath{\ensuremath{\ensuremath{\mathbf{\nabla}}}\, \psi}) = -\varrho\ ,$$ (2.1)

$$\ensuremath{\ensuremath{\ensuremath{\mathbf{\nabla}}}\cdot \ensur... ...splaystyle \frac{\displaystyle \partial n}{\displaystyle \partial t} \right)\ ,$$ (2.2)

$$\ensuremath{\ensuremath{\ensuremath{\mathbf{\nabla}}}\cdot \ensur... ...splaystyle \frac{\displaystyle \partial p}{\displaystyle \partial t} \right)\ .$$ (2.3)

The unknown quantities of this equation system are the electrostatic potential $\psi$, and the electron and hole concentrations $n$ and $p$, respectively. $\mathrm{q}$ is the elementary charge, $\varrho$ denotes the space charge density and $R$ stands for the net recombination rate, which is defined as

$$R = R_n - G_n = R_p - G_p \ .$$ (2.4)

These three equations can be derived from the four Maxwell equations [100]:

$$\ensuremath{\ensuremath{\ensuremath{\mathbf{\nabla}}}\ensuremath{... ...c{\displaystyle \partial \ensuremath{\mathbf{B}}}{\displaystyle \partial t} \ ,$$ (2.5)

$$\ensuremath{\ensuremath{\ensuremath{\mathbf{\nabla}}}\cdot \ensuremath{\mathbf{B}}} = 0 \ ,$$ (2.6)

$$\ensuremath{\ensuremath{\ensuremath{\mathbf{\nabla}}}\ensuremath{... ...c{\displaystyle \partial \ensuremath{\mathbf{D}}}{\displaystyle \partial t} \ ,$$ (2.7)

$$\ensuremath{\ensuremath{\ensuremath{\mathbf{\nabla}}}\cdot \ensuremath{\mathbf{D}}} = \varrho \ .$$ (2.8)

$\ensuremath{\mathbf{E}}$ and $\ensuremath{\mathbf{D}}$ are the electric field and displacement vectors, respectively, and $\ensuremath{\mathbf{H}}$ and $\ensuremath{\mathbf{B}}$ are the magnetic field and induction vectors, respectively. $\ensuremath{\mathbf{J}}$ is the conduction current density.

$\ensuremath{\mathbf{D}}$ is related to the electric field $\ensuremath{\mathbf{E}}$ by the permittivity tensor $\ensuremath{{\mathrm{\hat{\varepsilon}}}}$ (assumed to be a scalar $\varepsilon$ hereafter):

$$\ensuremath{\mathbf{D}} = \ensuremath{{\mathrm{\hat{\varepsilon}}}}\ensuremath{\mathbf{E}} \ ,$$ (2.9)

which is valid for all materials which have a frequency-independent permittivity and do not have piezoelectric or ferroelectric effects. The Poisson equation can be derived by introducing a vector potential $\ensuremath{\mathbf{A}}$ defined by $\ensuremath{\mathbf{B}} = \ensuremath{\ensuremath{\ensuremath{\mathbf{\nabla}}}\ensuremath{\times}\ensuremath{\mathbf{A}}}$, which is inserted in (2.5):

$$\ensuremath{\ensuremath{\ensuremath{\mathbf{\nabla}}}\ensuremath{... ...h{\ensuremath{\mathbf{\nabla}}}\ensuremath{\times}\ensuremath{\mathbf{A}}}) \ .$$ (2.10)

After interchanging the order of the time derivative and the curl operator and using the associative property of the curl operator, the argument of the curl operator can be substituted by the gradient of a scalar potential, because $\ensuremath{\ensuremath{\ensuremath{\mathbf{\nabla}}}\ensuremath{\times}\ensuremath{\ensuremath{\ensuremath{\mathbf{\nabla}}}\, \psi}}$ yields zero:

$$\ensuremath{\mathbf{E}} + \displaystyle \frac{\displaystyle \part... ...artial t} = - \ensuremath{\ensuremath{\ensuremath{\mathbf{\nabla}}}\, \psi} \ .$$ (2.11)

Since the wavelength is much larger than the typical semiconductor device dimensions, the quasi-stationary approximation can be assumed. With an operating frequency $\ensuremath{f_\mathrm{op}}= 200\,$GHz, $c_0$ as the speed of light in vacuum and $\varepsilon_\mathrm{r} = 11.9$ [217] and $\mu_\mathrm{r} = 1$ as the relative permittivity and permeability in silicon, respectively, $\lambda = c_0 / (\sqrt{\varepsilon_\mathrm{r} \,\mu_\mathrm{r}} {\ensuremath{f_\mathrm{op}}})$ equals $435\, \mu$m [91]. The typical gate length of industrially produced devices is already below $100\,$nm. In the quasi-stationary approximation, the time derivative of the vector potential can be neglected. By inserting $\ensuremath{\mathbf{E}} = - \ensuremath{\ensuremath{\ensuremath{\mathbf{\nabla}}}\, \psi}$ into (2.9) and by inserting the result into (2.8), the Poisson equation (2.1) is finally obtained. The space charge density $\varrho$ in semiconductors is composed of the charges of electrons, holes, the ionized dopant atoms and other charged defects:

$$\varrho = \mathrm{q}(p - n + C)\ . %\left(p + N_D^+\right) - \q \, \left(n + N_A^-\right) \ ,$$ (2.12)

The fixed charge $C$ is commonly modeled as

$$C = N_\mathrm{D}^+ - N_\mathrm{A}^- \ ,$$ (2.13)

where $N_\mathrm{D}^+$ is the concentration of the positively charged donor atoms and $N_\mathrm{A}^-$ the concentration of the negatively charged acceptor atoms. Dynamic recombination processes such as the Shockley-Read-Hall recombination require an additional dynamically changing trap charge $N_t$. It is important to note that not all dopants are electrically active and that not all electrically active dopants are always ionized. Particularly at low temperatures this fact has to be included in the model. However, in silicon at room temperature $N_D^+ = N_D$ and $N_A^- = N_A$ is assumed.

The Poisson equation without considering a magnetic field is finally obtained by

$$\ensuremath{\ensuremath{\ensuremath{\mathbf{\nabla}}}\cdot (}\var... ...}\, \psi}) = \mathrm{q}\left(n - p + N_\mathrm{A}^- - N_\mathrm{D}^+\right) \ .$$ (2.14)

The continuity equation for the conduction current density is derived by applying the divergence operator to (2.7) and using (2.8):

$$\ensuremath{\ensuremath{\ensuremath{\mathbf{\nabla}}}\cdot \ensur... ...ystyle \frac{\displaystyle \partial \varrho}{\displaystyle \partial t}} = 0 \ .$$ (2.15)

The equation can be interpreted, that all sources and sinks of the conduction current density are compensated by the time variation of the space charge density.

Whereas the Maxwell equations can be used to derive the Poisson equation and the current continuity equations, the current relations cannot be derived from them. The carrier transport equations are therefore discussed in the next sections.

2.1.1 The Drift-Diffusion Transport Model

Many causes of the current flows can be identified, for example gradients of the carrier concentrations, temperatures or material properties or a contribution determined by Ohm's law [137,138]. The latter component is called drift current and can be formulated as

$$\ensuremath{\mathbf{J}}^\mathrm{drift} = \sigma \ensuremath{\mathbf{E}} \ .$$ (2.16)

With $\sigma_n = \mathrm{q}n \mu_n$ as the electrical conductivity, $\mu_n$ as the mobility of electrons, and $\ensuremath{\mathbf{V}}_n = - \mu_n \ensuremath{\mathbf{E}}$ as the mean velocities of electrons, the current density is obtained by

$$\ensuremath{\mathbf{J}}^\mathrm{drift}_n = \mathrm{q}n \mu_n \ensuremath{\mathbf{E}} = - \mathrm{q}n \ensuremath{\mathbf{V}}_n \ .$$ (2.17)

Whereas the electric field accelerates the carriers, the velocity of the carriers is limited by various scattering mechanisms such as lattice vibrations (phonon scattering) [31]. The saturation velocity $V_\mathrm{sat}$ is generally modeled independently of the doping, whereas below $V_\mathrm{sat}$ an indirect relation is encountered. The diffusion current is caused by the thermal motion of the carriers described by the gradient of the carrier concentration. Derived from the law of diffusion, the following relation can be given:

$${ \ensuremath{\mathbf{J}}^\mathrm{diffusion}_n = \mathrm{q}D_n \ensuremath{\ensuremath{\ensuremath{\mathbf{\nabla}}}\, n} \ . }$$ (2.18)

$D_n$ is the diffusion coefficient for electrons which is often modeled using the mobility by the Einstein relation

$$D_n = \mu_n \frac{\ensuremath{\mathrm{k_B}}T_n}{\mathrm{q}}\ ,$$ (2.19)

where $\ensuremath{\mathrm{k_B}}$ is the Boltzmann constant. This relation is valid for conditions close to thermal equilibrium and for non-degenerate carrier systems (Boltzmann statistics). Both current components are added to obtain the isothermal drift-diffusion transport model

$$\ensuremath{\mathbf{J}}_n = \mathrm{q}D_n \ensuremath{\ensuremath... ...} - \mu_ n n \ensuremath{\ensuremath{\ensuremath{\mathbf{\nabla}}}\, \psi}) \ ,$$ (2.20)

which is applied in (2.2) to obtain the current continuity equation. Together with the Poisson equation (2.1), they form the basic semiconductor equations as given first by VanRoosbroeck [225]:

$${ \ensuremath{\ensuremath{\ensuremath{\mathbf{\nabla}}}\cdot (\va... ...ilon \ensuremath{\ensuremath{\ensuremath{\mathbf{\nabla}}}\, \psi})} = -\rho\ ,$$ (2.21)

$$\ensuremath{\ensuremath{\ensuremath{\mathbf{\nabla}}}\cdot \left(... ... R + \displaystyle \frac{\displaystyle \partial n}{\displaystyle \partial t}\ ,$$ (2.22)

$$\ensuremath{\ensuremath{\ensuremath{\mathbf{\nabla}}}\cdot \left(... ...R + \displaystyle \frac{\displaystyle \partial p}{\displaystyle \partial t}\ .{$$ (2.23)

To this set the lattice heat-flow equation (2.24) is frequently added to account for thermal effects such as self-heating in the device:

$$\ensuremath{\ensuremath{\ensuremath{\mathbf{\nabla}}}\cdot (}\kap... ... \frac{\displaystyle \partial T_{\mathrm{L}}}{\displaystyle \partial t} - H \ .$$ (2.24)

This equation requires modeling of the thermal conductivity $\kappa_{\mathrm{L}}$, the mass density $\rho_{\mathrm{L}}$, the heat capacity $c_{\mathrm{L}}$, and the locally generated heat $H$. As the parameters of equations (2.21) to (2.23) depend also on the lattice temperature and a gradient in the lattice temperature causes a current flow [227], the thermal drift-diffusion carrier continuity equations for electrons and holes are given as [194]

$$\ensuremath{\mathbf{J}}_n = \ensuremath{\ensuremath{\ensuremath{\... ... R + \displaystyle \frac{\displaystyle \partial n}{\displaystyle \partial t}\ ,$$ (2.25)

$$\ensuremath{\mathbf{J}}_p = \ensuremath{\ensuremath{\ensuremath{\... ... R + \displaystyle \frac{\displaystyle \partial p}{\displaystyle \partial t}\ ,$$ (2.26)

forming a coupled system with (2.24).

2.1.2 Types of Partial Differential Equations

Nonlinear partial differential equations of second-order can appear in three variants: elliptic, parabolic, and hyperbolic. The Poisson equation as well as the steady-state continuity equations form a system of elliptic partial differential equations, whereas the lattice heat-flow equation is parabolic.

This can be shown for the Poisson equation posed in a bounded domain $D \in R^{2}$ by rewriting the differential operators of (2.21) in (2.27) and for Cartesian coordinates in (2.28), which is then compared with a general partial differential equation of second order (2.29):

$${ \nabla^2 \psi = - \varrho/\varepsilon\ ,$$ (2.27)

$$\partial_x^2u + \partial_y^2u = - \varrho/\varepsilon\ ,$$ (2.28)

$$Au_{xx} + 2Bu_{xy} + Cu_{yy} + Du_{x} + Eu_{y} + F_u = G \ .{$$ (2.29)

The coefficients $A$ to $G$ are piecewise continuous functions of $x$ and $y$. In analogy to quadratic forms, one can calculate the determinant

$$d = AC - B^2 \ ,$$ (2.30)

and classify the given equation: in case $d > 0$ the equation is elliptic (for the Poisson equation: $A = C = 1, B = 0$), in case $d = 0$ parabolic and in case $d < 0$ hyperbolic.

To completely determine the solution of an elliptic partial differential equation boundary conditions have to be specified. Since parabolic and hyperbolic partial differential equations describe evolutionary processes, time normally is an independent variable and an initial condition is additionally required. As the transient continuity equations are parabolic, they require such an initial condition.

2.1.3 Systematic Derivation of the Macroscopic Transport Models

Macroscopic transport models can be systematically derived from the Boltzmann transport equation [134]

$$\displaystyle \frac{\displaystyle \partial f}{\displaystyle \part... ..._{\negmedspace \ensuremath{\mathbf{p}}}} f} = Q[f(\ensuremath{\mathbf{k}})] \ .$$ (2.31)

This equation is a time-dependent partial integro-differential equation in the six-dimensional phase space $(\ensuremath{\mathbf{k}}, \ensuremath{\mathbf{r}})$, which assumes that the carrier motion is determined by Newton's laws. The solution is the distribution function $f$, $\ensuremath{\mathbf{u}}$ denotes the group velocity, $\ensuremath{\mathbf{F}}$ stands for the external force. The scattering operator $Q[\cdot]$, which is in general nonlinear in $f$, represents the rate of change of $f$ due to collisions and is modeled via Fermi's Golden Rule.

The distribution function $f$ as solution of (2.31) is used to obtain the probability of finding a carrier inside a phase-space volume $\mathrm{d}\ensuremath{\mathbf{k}}\mathrm{d}\ensuremath{\mathbf{r}}$. The equilibrium distribution function is the Fermi-Dirac distribution. If the Pauli exclusive principle is neglected, the Maxwell distribution is obtained. In the drift-diffusion model the cold and in the energy-transport the heated Maxwell distribution is used. The six moments model uses an analytical distribution function in the approach pursued here [84].

This can be performed either directly by Monte Carlo simulations [112,115,123,122] or by the methods based on the expansion of the distribution function in momentum space into a series of spherical harmonics [97].

Monte Carlo simulations are particularly useful to obtain a physically accurate picture including various effects with as few approximations as possible. Since the simulations require much computational resources in terms of time and memory, it is rarely used on an engineering level. The extraction of small-signal parameters is often prohibitively costly, because the evaluation of a small-signal perturbation requires a low variance in the results.

Whereas for extremely small devices with gate lengths below 10$\,$nm systems based on the Wigner-Boltzmann equation [151,149,150,70] have been applied, the ballistic limit has been investigated by using transport models based on the Schrödinger [124] equation.

As its direct solution requires considerable computational resources, simpler solutions are obtained by investigating lower moments of the distribution function only. These include the electron concentration $n$, the electron temperature $\ensuremath{T_n}$, and the electron kurtosis $\ensuremath{\beta_n}$.

By multiplying (2.31) with a weight function $\phi (\ensuremath{\mathbf{k}})$ and integrating the result over the $\ensuremath{\mathbf{k}}$ space, the macroscopic transport equations are obtained. It is assumed that the Brillouin zone extends towards infinity, which is justified due to the exponential decline of the distribution function [240]. This procedure results in the following partial differential equations in ( $\ensuremath{\mathbf{r}}, t$)

$$\displaystyle \frac{\displaystyle \partial n \ensuremath{\langle ... ...emath{\mathbf{k}})] \,\, \mathrm{d}^3 \ensuremath{\mathbf{k}}} \ ,%\bestkiiint$$ (2.32)

while the coordinates of the $\ensuremath{\mathbf{k}}$ space are saturated. As macroscopic quantities describing the average behavior of the microscopic quantity $\phi (\ensuremath{\mathbf{k}})$, moments neither depend on $\ensuremath{\mathbf{k}}$ nor contain the respective information any more.

For the balance and flux equations, the weight functions ${\mathcal{E}}^i$ and $\ensuremath{\mathbf{p}}\,{\mathcal{E}}^i$ with the momentum $\ensuremath{\mathbf{p}}= \hbar \, \ensuremath{\mathbf{k}}$ are often used [84]. For the drift-diffusion transport model, the expansion is truncated at $i=1$, for the energy-transport transport model at $i=2$, and eventually for the six moments transport model at $i=3$. The balance and flux relations read as follows [83]:

$$\phi _0 = {\mathcal{E}}^0 = 1 \qquad \Rightarrow \qquad \displaystyle \frac{\displaystyle \partial n }{\displaystyle \partial t} \ +\ \ensuremath{\ensuremath{\mathbf{\nabla}}}\cdot (n \ensuremath{\mathbf{V}}\!_0 ) = \ n q_0 \$$ (2.33)

$$\phi _1 = \ensuremath{\mathbf{p}}{\mathcal{E}}^0 = \ensuremath{\mathbf{p}} \qquad \Rightarrow \qquad \displaystyle \frac{\displaystyle \partial n \ensuremath{\mathbf{P}}\!_0}{\displaystyle \partial t} \ +\ \ensuremath{\ensuremath{\mathbf{\nabla}}}\cdot (n \ensuremath{{\mathrm{\hat{U}}}}_{1}) \ - \ n \ensuremath{\mathbf{F}} = \ n \ensuremath{\mathbf{Q}}_0 \$$ (2.34)

$$\phi _2 = {\mathcal{E}}^1 \qquad \Rightarrow \qquad \displaystyle \frac{\displaystyle \partial n w_1}{\displaystyle \partial t} \ +\ \ensuremath{\ensuremath{\mathbf{\nabla}}}\cdot (n \ensuremath{\mathbf{V}}\!_1 ) \ - \ 1\ n \ensuremath{\mathbf{F}}\cdot \ensuremath{\mathbf{V}}\!_{0} = \ n q_1 \$$ (2.35)

$$\phi _3 = \ensuremath{\mathbf{p}}{\mathcal{E}}^1 \qquad \Rightarrow \qquad \displaystyle \frac{\displaystyle \partial n \ensuremath{\mathbf{P}}\!_1}{\displaystyle \partial t} \ +\ \ensuremath{\ensuremath{\mathbf{\nabla}}}\cdot (n \ensuremath{{\mathrm{\hat{U}}}}_{2}) \ - \ n \ensuremath{\mathbf{F}}\cdot ( w_1 \, \ensuremath{{\mathrm{\hat{I}}}}+ \, \ensuremath{{\mathrm{\hat{U}}}}_{1} ) = \ n \ensuremath{\mathbf{Q}}_1 \$$ (2.36)

$$\phi _4 = {\mathcal{E}}^2 \qquad \Rightarrow \qquad \displaystyle \frac{\displaystyle \partial n w_2}{\displaystyle \partial t} \ +\ \ensuremath{\ensuremath{\mathbf{\nabla}}}\cdot (n \ensuremath{\mathbf{V}}\!_2 ) \ - \ 2\ n \ensuremath{\mathbf{F}}\cdot \ensuremath{\mathbf{V}}\!_{1} = \ n q_2 \$$ (2.37)

$$\phi _5 = \ensuremath{\mathbf{p}}{\mathcal{E}}^2 \qquad \Rightarrow \qquad \displaystyle \frac{\displaystyle \partial n \ensuremath{\mathbf{P}}\!_2}{\displaystyle \partial t} \ +\ \ensuremath{\ensuremath{\mathbf{\nabla}}}\cdot (n \ensuremath{{\mathrm{\hat{U}}}}_{3}) \ - \ n \ensuremath{\mathbf{F}}\cdot ( w_2 \, \ensuremath{{\mathrm{\hat{I}}}}+ 2 \, \ensuremath{{\mathrm{\hat{U}}}}_{2} ) = \ n \ensuremath{\mathbf{Q}}_2 \$$ (2.38)

$$\hspace*{4.5mm}\vdots \hspace*{2.3mm}\vdots \$$

with

$$\ensuremath{\mathbf{P}}\!_i = \ensuremath{\langle \ensuremath{\ma... ...{u}}\ensuremath{\otimes}\ensuremath{\mathbf{p}}{\mathcal{E}}^{i-1} \rangle} \ .$$ (2.39)

The moments of the scattering integral are

$$n q_i = \ensuremath{\int {\mathcal{E}}^i \, Q[f(\ensuremath{\mathbf{k}})] \,\, \mathrm{d}^3 \ensuremath{\mathbf{k}}} \ ,$$ (2.40)

$$n \ensuremath{\mathbf{Q}}_i = \ensuremath{\int \ensuremath{\mathb... ...\, Q[f(\ensuremath{\mathbf{k}})] \,\, \mathrm{d}^3 \ensuremath{\mathbf{k}}} \ .$$ (2.41)

$\ensuremath{\mathbf{F}}$ stands for the external force, which is given as

$$\ensuremath{\mathbf{F}}(\ensuremath{\mathbf{r}}) = - \mathrm{q}\, \ensuremath{\mathbf{E}}(\ensuremath{\mathbf{r}})\ ,$$ (2.42)

for homogeneous band structures and neglected magnetic fields. In addition, the system is claimed to be diffusion dominated [174]. As a consequence, the convective terms and the time derivatives in the flux relations are neglected resulting in parabolic partial differential equation systems [89].

The moment $n$ is the carrier concentration, $n \ensuremath{\mathbf{V}}\!_0$ the average carrier flux, $w_1$ the average energy, $n \ensuremath{\mathbf{V}}\!_1$ the average energy flux, $w_2$ the average energy squared, and $n \ensuremath{\mathbf{V}}\!_2$ the average kurtosis flux. $\ensuremath{\mathbf{V}}\!_0$ is the average carrier velocity, $\ensuremath{\mathbf{P}}\!_0$ is the average momentum, $\ensuremath{{\mathrm{\hat{U}}}}_1$ the energy tensor, $\ensuremath{{\mathrm{\hat{U}}}}_2$ the second energy tensor, and $\ensuremath{{\mathrm{\hat{U}}}}_3$ is eventually the third energy tensor. Furthermore it can be noted that (2.33) - (2.38) are the continuity equation and current relation, the energy balance and energy-flux relation, and the kurtosis balance and kurtosis-flux relation, respectively.

As this equation system contains more unknowns than equations and each equation is coupled to the next higher one, the so-called closure problem comes into play: additional relations have to be introduced in order to close the equation system, thus making it solvable. There are several approaches to overcome this problem [84], for example the maximum entropy method [129], Grad's method [80], the cumulant closure method [176,235,192], or a Maxwell closure [85]. The choice of variables which are expressed by other ones is critical. The current density is proportional to the average velocity $\ensuremath{\mathbf{V}}\!_0$. Furthermore, the average energy $w_1$ is used to model many physical effects. For that reason, $\ensuremath{\mathbf{V}}\!_i$ and $w_i$ are chosen to remain as solution variables and the additional moments are expressed as functions of the solution variables. According to [82], the energy-like tensor $\ensuremath{{\mathrm{\hat{U}}}}_i$ are

$$\ensuremath{{\mathrm{\hat{U}}}}_i = \ensuremath{\langle \ensurema... ...\ensuremath{{\mathrm{\hat{I}}}}+ U_i^{(2)}(w_i, \ensuremath{\mathbf{V}}\!_i)\ ,$$ (2.43)

where the anisotropic second term is neglected. By applying the macroscopic relaxation time approximation [87,82] the equation system in the form of (2.44) and (2.45) is obtained.

$${ \displaystyle \frac{\displaystyle \partial n w_i}{\displaystyle... ...) - i \ensuremath{\mathbf{F}}\cdot n \ensuremath{\mathbf{V}}\!_{i-1} = n q_i\ ,$$ (2.44)

$$\ensuremath{\ensuremath{\mathbf{\nabla}}}\cdot (n \ensuremath{{\m... ...+ i \, \ensuremath{{\mathrm{\hat{U}}}}_{i} ) = n \ensuremath{\mathbf{Q}}_i \ ,{$$ (2.45)

with

$$q_i = - \frac{w_i - w_{i,\mathrm{eq}}}{\tau_i} \ ,$$ (2.46)

$$\ensuremath{\mathbf{Q}}_i = - \frac{\ensuremath{\mathbf{V}}\!_i}{\mu_i}\ .$$ (2.47)

The fluxes can be explicitly written as

$$n \ensuremath{\mathbf{V}}\!_i = - \frac{2 \mu_i H_{i+1}}{3 \mathr... ...) - n \, \ensuremath{\mathbf{F}}w_i \frac{ 3 + 2 i H_{i}}{2 H_{i+1}} \Bigl) \ .$$ (2.48)

The non-parabolicity correction factors $H_i$ equal unity in the case of parabolic bands and are modeled as either energy-dependent using a simple analytical expression [28], by the incorporation of bulk Monte Carlo data [84,219] or via analytic models for the distribution function [82].

Note that (2.45) still contains the tensors $\ensuremath{{\mathrm{\hat{U}}}}_i$ which have to be approximated by the unknowns of the equation system $w_i$ and $\ensuremath{\mathbf{V}}\!_i$. Finally, the following variable transformation according to [86] is introduced:

$$w_1 = \frac{3}{2} \ensuremath{\mathrm{k_B}}T \ , \qquad w_2 = \underbrace{(\frac{15}{4} {\mathrm{k_B^2}} T^2)}_ \beta \ , \qquad w_3 = \underbrace{(\frac{105}{8} {\mathrm{k_B^3}} T^3)}_ \gamma .$$ (2.49)

$\beta$ and $\gamma$ are the contributions of the non-Maxwell closure. Thus, in order to obtain the Maxwell closure, $\beta = 1$ and $\gamma = 1$ is used. The quantity $\beta$ is the kurtosis of the distribution function and indicates the deviation from the Maxwell shape. The quantity $\gamma$ depends on the applied closure relation and equals $\beta^{c}$ for the generalized Maxwell closure, which will be used in the following. Then, the following form of the flux equations can be obtained:

$${ n \ensuremath{\mathbf{V}}\!_0 = - A_0 \Bigl( \ensuremath{\ensur... ...nabla}}}(n \ensuremath{\mathrm{k_B}}T) - n \ensuremath{\mathbf{F}}h_0 \Bigr)\ ,$$ (2.50)

$$n \ensuremath{\mathbf{V}}\!_1 = - A_1 \Bigl( \ensuremath{\ensurem... ... T^2 \beta ) - n \ensuremath{\mathbf{F}}h_1 \ensuremath{\mathrm{k_B}}T\Bigr)\ ,$$ (2.51)

$$n \ensuremath{\mathbf{V}}\!_2 = - A_2 \Bigl( \ensuremath{\ensurem... ...3 \gamma) - n \ensuremath{\mathbf{F}}h_2 {\mathrm{k_B^2}} T^2 \beta \Bigr)\ , {$$ (2.52)

with

$$A_i = \frac{\mu_i H_{i+1}}{2^i \mathrm{q}} \frac{1}{3} \prod_{j=0}^{i+1}(1+2j) \qquad h_i = \frac{3 + 2 i H_i}{(3 + 2 i)H_{i+1}} \ .$$ (2.53)

The balance equations are derived as

$${ \displaystyle \frac{\displaystyle \partial n}{\displaystyle \pa... ...th{\ensuremath{\mathbf{\nabla}}}\cdot (n \ensuremath{\mathbf{V}}\!_0) \ = 0 \ ,$$ (2.54)

$$C_1 \displaystyle \frac{\displaystyle \partial n T}{\displaystyle... ...\ensuremath{\mathbf{V}}\!_{0} = - n \, C_1 \frac{T - T_\mathrm{eq}}{\tau_1} \ ,$$ (2.55)

$$C_2 \displaystyle \frac{\displaystyle \partial n T^2 \beta}{\disp... ...- n \, C_2 \frac{T^2 \beta - T^2_\mathrm{eq} \, \beta_\mathrm{eq}}{\tau_2} \ ,{$$ (2.56)

with $C_1 = 3 \ensuremath{\mathrm{k_B}}/2$ and $C_2 = 15 {\mathrm{k_B^2}} / 4$. By using the equilibrium solution of the Boltzmann transport equation, the equilibrium values $T_\mathrm{eq}$ and $\beta_\mathrm{eq}$ can be calculated for the Maxwell distribution.

The equilibrium carrier temperature defined by (2.49) is different from the lattice temperature, because a non-parabolic band structure is used. One obtains $T_\mathrm{eq} \approx 309.452\,\mathrm{K}$ and $\beta_\mathrm{eq}\approx 1$ [84]. These are also the values used for the Dirichlet boundary conditions for $T$ and $\beta$. For the carrier concentration a standard Ohmic contact model is used, corresponding to a cold Maxwell distribution at the contacts.

In the following sections, the transport models for the parabolic case are used. Thus, $h_i = H_i = 1$, $T_\mathrm{eq} = T_{\mathrm{L}}$, and $\beta_\mathrm{eq} = 1$, resulting in

$$A_0 = \frac{\mu_0}{\mathrm{q}} \ , \qquad A_1 = \frac{5\,\mu_1}{2\,\mathrm{q}} \ , \qquad A_2 = \frac{35\,\mu_2}{4\,\mathrm{q}} \ .$$ (2.57)

2.1.4 The Six Moments Transport Model

The closure of the six moments transport model at $\phi _6$ was determined as $\gamma = \beta^{c}$. Thus, the following flux equations can be derived:

$${ \ensuremath{\mathbf{J}}_n = -\mathrm{q}(n \ensuremath{\mathbf{V... ...f{\nabla}}}(n \ensuremath{\mathrm{k_B}}T) - n \ensuremath{\mathbf{F}}\Bigr) \ ,$$ (2.58)

$$\ensuremath{\mathbf{S}}_n = (n \ensuremath{\mathbf{V}}\!_1) = \ph... ...2}} T^2 \beta ) - n \ensuremath{\mathbf{F}}\ensuremath{\mathrm{k_B}}T\Bigr) \ ,$$ (2.59)

$$\ensuremath{\mathbf{K}}_n = (n \ensuremath{\mathbf{V}}\!_2) = \ph... ...^3 \beta^{c}) - n \ensuremath{\mathbf{F}}{\mathrm{k_B^2}} T^2 \beta \Bigr) \ .{$$ (2.60)

The balance equations are

$${ \displaystyle \frac{\displaystyle \partial n}{\displaystyle \pa... ...th{\ensuremath{\mathbf{\nabla}}}\cdot (n \ensuremath{\mathbf{V}}\!_0) \ = 0 \ ,$$ (2.61)

$$\frac{3 \, \ensuremath{\mathrm{k_B}}}{2} \displaystyle \frac{\dis... ...rac{3 \, \ensuremath{\mathrm{k_B}}}{2} \, \frac{T - T_{\mathrm{L}}}{\tau_1} \ ,$$ (2.62)

$$\frac{15 \, {\mathrm{k_B^2}}}{4} \displaystyle \frac{\displaystyl... ...{15 \, {\mathrm{k_B^2}}}{4} \, \frac{T^2 \beta - T_{\mathrm{L}}^2}{\tau_2} \ .{$$ (2.63)

With $c = 2.7$ the best match for $w_3$ was obtained in comparison to Monte Carlo simulations [83]. The boundary conditions are derived by applying a cold Maxwell distribution:

$${ f_{\mathrm{eq}} = A \exp{\left(-\frac{{\mathcal{E}}}{\ensuremath{\mathrm{k_B}}T_{\mathrm{L}}}\right)}\ ,$$ (2.64)

$$n = \int f \, g \, \mathrm{d}{{\mathcal{E}}}\ , {$$ (2.65)

resulting in the boundary conditions for the higher moments

$$w_1 = {{\mathcal{E}}}_\mathrm{eq} = \frac{1}{n} \int {{\mathcal{E}}} f\,g \, \mathrm{d}{{\mathcal{E}}} = \frac{3}{2} \ensuremath{\mathrm{k_B}}T\ ,$$ (2.66)

$$w_2 = \frac{1}{n} \int {{\mathcal{E}}^2} f\,g \, \mathrm{d}{{\mathcal{E}}} = \frac{3 \cdot 5}{2 \cdot 2} \, {\mathrm{k_B^2}} T^2\ .$$ (2.67)

2.1.5 The Energy-Transport Model

By using the closure $\beta = 1$ at $\phi _4$, the flux equations of the energy-transport transport model can be obtained:

$${ \ensuremath{\mathbf{J}}_n = -\mathrm{q}(n \ensuremath{\mathbf{V... ...f{\nabla}}}(n \ensuremath{\mathrm{k_B}}T) - n \ensuremath{\mathbf{F}}\Bigr) \ ,$$ (2.68)

$$\ensuremath{\mathbf{S}}_n = (n \ensuremath{\mathbf{V}}\!_1) = \ph... ...m{k_B^2}} T^2)- n \ensuremath{\mathbf{F}}\ensuremath{\mathrm{k_B}}T\Bigr) \ . {$$ (2.69)

The balance equations are (2.61) and (2.62).

2.1.6 The Drift-Diffusion Transport Model

For the drift-diffusion transport model, the closure at $\phi _2$ was found to be $\ensuremath{T_n}=T_{\mathrm{L}}$, yielding the following balance and flux equation:

$${ \ensuremath{\mathbf{J}}_n = -\mathrm{q}(n \ensuremath{\mathbf{V... ... \ensuremath{\mathrm{k_B}}T_{\mathrm{L}}) - n \ensuremath{\mathbf{F}}\Bigr) \ ,$$ (2.70)

$$\displaystyle \frac{\displaystyle \partial n}{\displaystyle \part... ...nsuremath{\ensuremath{\mathbf{\nabla}}}\cdot \ensuremath{\mathbf{J}}_n = 0 \ .{$$ (2.71)

These equations look very familiar, as they were the result of the derivation in Section 2.1.1. In the flux equation (2.70), the external force is substituted by $\ensuremath{\mathbf{F}}= -\mathrm{q}\, \ensuremath{\mathbf{E}} = \mathrm{q}\, \ensuremath{\ensuremath{\ensuremath{\mathbf{\nabla}}}\, \psi}$.

Eventually, (2.70) must be inserted into (2.71), which is the balance equation corresponding to (2.2).

$$\ensuremath{\ensuremath{\mathbf{\nabla}}}\cdot \ensuremath{\mathb... ...rm{q}\, \displaystyle \frac{\displaystyle \partial n}{\displaystyle \partial t}$$ (2.72)

$$\ensuremath{\ensuremath{\mathbf{\nabla}}}\cdot \Bigl[-\mathrm{q}\... ...rm{q}\, \displaystyle \frac{\displaystyle \partial n}{\displaystyle \partial t}$$ (2.73)

When $T_{\mathrm{L}}$ is assumed to be constant, the isothermal drift-diffusion model as already shown in (2.22) can be finally obtained:

$$\ensuremath{\ensuremath{\mathbf{\nabla}}}\cdot \Bigl(\underbrace{... ...\displaystyle \frac{\displaystyle \partial n}{\displaystyle \partial t} + R \ ,$$ (2.74)

with the net recombination rate $R = R_0$ according to

$$n q_0 = n \frac{\ensuremath{\langle 1 \rangle} - \ensuremath{\langle 1 \rangle_\mathrm{eq}}}{\tau_0} + R_0 = R_0$$ (2.75)