(image) (image) [Previous] [Next]

Predictive and Efficient Modeling of Hot Carrier Degradation with Drift-Diffusion Based Carrier Transport Models

4.2 Semi-Classical Transport

Let us consider a semiconductor device in thermal equilibrium when no external electric field is applied and all the transient processes have relaxed. In this situation, the average energy of the carriers is equal to the thermal energy, while their average velocity is zero. This corresponds to the absence of a macroscopic current in the device. When a bias is applied, the carriers experience a force and accelerate until opposed by other forces such as scattering. Thus, the motion of particles inside a semiconductor device is a result of forces due to the applied bias, built-in potential and scattering potential. Their resultant force causes the momentum and position of the charge carriers to change. Thus, within the classical regime, to find the effect of fields on charge carriers, their momentum and position need to be evaluated.

4.2.1 Monte-Carlo Approach

In the Monte-Carlo method for solving the carrier transport problem, the Newton equation is evaluated for the particles. For a single ( \( i^{\mathrm {th}} \)) particle, the system of equations would read:

(4.1) \begin{equation} \frac {\mathrm {\partial }p_{i}}{\mathrm {\partial }t}=F(p,r,t)+R(p,r,t),\label {eq:MC_1} \end{equation}

(4.2) \begin{equation} \frac {\mathrm {\partial }r_{i}}{\mathrm {\partial }t}=u_{i}(t),\label {eq:MC_2} \end{equation}

where \( F(p,r,t) \) is due to external forces and \( R(p,r,t) \) is due to local forces. In the semi-classical approach, the effect of external forces is treated classically, but the local forces such as scattering are treated quantum mechanically. The above equations, 4.1 and 4.2 describe a classical carrier transport and are not fully applicable to modern devices. The above is equivalent to the BTE, band structure enters through the dispersion relation and scattering goes to R. Moreover, it is very computationally demanding to obtain a solution for all particles in a typical semiconductor device, containing \( \sim 10^{20} \) particles.

4.2.2 Boltzmann Transport Equation

Another approach to describe carrier transport semi-classically is to consider the statistical properties of the charge carriers by defining a distribution function \( f(p,r,t) \). The carrier distribution function \( f(p,r,t) \) gives the probability of finding a carrier with momentum in the range of \( (p,p+dp) \) and volume of \( (r,r+dr) \) and satisfies the Boltzmann transport equation:

(4.3) \begin{equation} \partial _{t}f+u\nabla _{\mathrm {r}}f+F\nabla _{\mathrm {p}}f=Q(f),\label {eq:BTE} \end{equation}

where \( Q(f) \) is the scattering operator. The semi-classical nature is incorporated in the relation of the carrier momentum \( p \) with the wave number \( k \) as \( p=\hbar k \), \( \hbar \) being the reduced Planck constant. The group velocity \( u \) in Equation 4.3 is given by:

(4.4) \begin{equation} u(p,r)=\nabla _{\mathrm {p}}E_{\mathrm {K}}(p,r), \end{equation}

where \( E_{\mathrm {K}} \) represents the kinetic energy of the carriers. The force \( F(p,r) \) acting on the carriers due to the electric \( E(r) \) and magnetic fields \( B(r) \) and inhomogeneous material properties is [174]:

(4.5) \begin{equation} F(p,r)=-\nabla _{\mathrm {r}}\left (E_{\mathrm {c,0}}(r)+\varepsilon (p,r)\right )+q\left (E(r)+u(p,r)\times B(r)\right ), \end{equation}

where \( E_{\mathrm {c,0}}(r) \) is the energetic position of the bottom of the conduction band and \( \varepsilon (p,r) \) represents the band structure. Scattering affects the distribution function via in-scattering and out-scattering:

\begin{eqnarray*} Q(f)=\sum f(p)(1-f(p'))S(p,p')-f(p')(1-f(p))S(p',p), \end{eqnarray*}

where \( S(p,p') \) represents the transition rate for a particle from a state denoted by the momentum \( p \) to a state denoted by the momentum \( p' \), while \( 1-f(p') \) represents the probability that the state at \( p' \) is empty. For non-degenerate semiconductors (low carrier concentration), the Pauli principle can be neglected and thus the scattering operator becomes linear and has an easier structure:

(4.6) \begin{equation} Q(f)=\sum f(p)S(p,p')-f(p')S(p',p), \end{equation}

The function \( S \) is assumed to account for all scattering mechanisms. The rates for the scattering mechanisms are calculated using the overlap integrals and deformation potentials following Fermi’s Golden Rule [175].

Since the BTE is a seven dimensional integro-differential equation, it is very challenging to obtain an exact solution. The Monte-Carlo method is widely used to numerically solve the BTE, but this method requires huge computational resources especially to resolve the high-energy tail of the carrier distribution function. In this subsection more computationally efficient approaches to solve the Boltzmann transport equation are discussed. Among these are the expansion of the distribution function in momentum space into a series of spherical harmonics, the method of moments, i.e. six moments model, the hydrodynamics or energy transport scheme, and the drift-diffusion scheme which is commonly used by device engineers.

4.2.2.1 Spherical Harmonics Expansion

The spherical harmonics expansion (SHE) method relies on the expansion of the carrier distribution function \( f(x,k,t) \) into a series of spherical harmonics [18]:

(4.7) \begin{equation} f(x,k,t)=\sum _{l=0}^{\infty }\sum _{m=-l}^{l}f_{l,m}(x,\epsilon ,t)\mathrm {Y^{\mathit {l,m}}(\theta ,\varphi )},\label {eq:SHE} \end{equation}

where \( k \) is the three-dimensional wave vector which is expressed in terms of spherical coordinates \( \epsilon \), \( \theta \), \( \varphi \). Owing to its deterministic nature, the computational effort in the SHE method is reduced as compared to the stochastic Monte-Carlo method. However, Equation 4.7 is still five dimensional and requires a fair amount of computational resources. These high computational demands of the SHE approach are also due to the system of coupled partial differential equations which needs to be solved for each expansion coefficient ( \( f_{l,m} \)).

4.2.2.2 The Moments Method

The macroscopic models to solve the BTE are common simplifications where only a few moments of the carrier distribution function are considered, like the carrier concentration and the carrier temperature [175, 176, 177, 178]. The moments are macroscopic quantities derived by integrating the distribution function over the three dimensional \( k \) space with a suitable weight function \( \phi \left (k\right ) \). The \( i^{\mathrm {th}} \) moment of the BTE is given as [174]:

(4.8) \begin{equation} \left \langle \phi _{i}\right \rangle =\int \phi _{i}\left (k\right )f\left (k\right )d^{3}k. \end{equation}

This reduces the complexity of the seven dimensional BTE by eliminating the three \( k \) coordinates. Some information is lost in this simplification which may not be important for device engineering applications or for obtaining an approximate analysis of charge transport. The weight functions \( \phi _{i}\left (k\right ) \) are chosen as increasing powers of the vector \( k \) with some scaling factor to obtain physically meaningful quantities. However, each equation for the moment \( i \) contains the next higher order moment \( i+1 \), resulting in the variables outnumbering the equations. In order to close this system of equations, usually an a priori assumption about the shape of the distribution function is made, for example using a Maxwellian distribution. The first few moments important to this work are listed in Table 4.1.

Order Moment
\( \phi _{0} \) \( n \)
\( \phi _{1} \) \( \frac {3}{2}nk_{\mathrm {B}}T_{\mathrm {n}} \)
\( \phi _{2} \) \( \frac {5*3}{2*2}n\left (k_{\mathrm {B}}T_{\mathrm {n}}\right )^{2} \)
\( \phi _{3} \) \( \frac {7*5*3}{2*2*2}n\left (k_{\mathrm {B}}T_{\mathrm {n}}\right )^{3} \)

Table 4.1: The first four moments of the Maxwell distribution.

The derivation of the moments of a Maxwell distribution for closure is given in Appendix B. Applying the moments in Table 4.1 to the BTE, and using the diffusion1 and relaxation time approximations2, the transport equations can be formulated [179]. For electrons, the first three transport equations can be written as:

(4.9–4.11) \{begin}{align} \phi _{0}: & \qquad \partial _{t}n-\frac {1}{q}\nabla J_{\mathrm {n}}=-R,\label {eq:order_0}\\ \phi _{1}: & \qquad J_{\mathrm {n}}=\frac {q\tau _{\mathrm {m}}}{m}\left (\frac {2}{3}\nabla \left \langle \phi _{2}\right \rangle +qE\right ),\label {eq:order_1}\\ \phi _{2}: & \qquad \frac {3}{2}k_{\mathrm {B}}\partial _{t}\left (nT_{\mathrm {n}}\right )+\nabla \left \langle \phi _{3}\right \rangle -EJ_{\mathrm {n}}=-\frac {3}{2}k_{\mathrm {B}}n\frac {T_{\mathrm {n}}-T_{\mathrm {L}}}{\tau \varepsilon }+G\varepsilon n,\label {eq:order_2} \{end}{align}

where \( \tau _{\mathrm {m}} \) is the momentum relaxation time, \( T_{\mathrm {n}} \) and \( T_{\mathrm {L}} \) are the electron and lattice temperature, respectively, \( R \) the net recombination rate, and \( E \) the electric field.

1 The diffusion approximation, implying that the diffusion of carriers dominates over drift, leads to the conclusion that the antisymmetric part of the distribution function is negligible and the symmetric part is isotropic.

2 The relaxation time approximation suggests that the perturbed distribution function will relax exponentially into the equilibrium distribution function.

4.2.2.3 Drift-diffusion

The drift-diffusion transport model has been the workhorse of device simulation for most device engineering applications. The DD model treats the transport either phenomenologically, i.e., carriers drift in electric fields and diffuse in concentration gradients or uses the first two moments of the BTE and hence the two transport Equations 4.9 and 4.10. The equation system has to be closed at \( \phi _{2} \) for which the relation \( \phi _{2}=\nicefrac {3}{2}\,nk_{\mathrm {B}}T_{\mathrm {L}} \) is used. It is worth mentioning that the DD equation system does not provide any information about the average carrier energy or carrier temperature \( T_{\mathrm {n}} \). Since, a thermal equilibrium condition is assumed treating \( T_{\mathrm {n}}=T_{\mathrm {L}} \). The DD model is then expressed as:

(4.12–4.13) \{begin}{align} \nabla J_{n}= & q\left (R+\partial _{t}n\right ),\\ J_{n}= & \mu _{n}k_{\mathrm {B}}\left (\nabla \left (nT_{\mathrm {L}}\right )+\frac {q}{k_{\mathrm {B}}}En\right ),\nonumber \\ = & qn\mu _{n}E+qD_{n}\nabla n,\label {eq:_DD_n} \{end}{align}

Similarly for holes, the DD equations read:

(4.14–4.15) \{begin}{align} \nabla J_{p}= & -q\left (R+\partial _{t}p\right ),\\ J_{p}= & qp\mu _{p}E-qD_{p}\nabla p,\label {eq:DD_p} \{end}{align}

where \( \mu _{n/p}=q\tau _{n/p}/m_{n/p} \) is the mobility of the respective charge carriers in the device with \( \tau _{n/p} \) being their momentum relaxation time and \( m_{n/p} \) their mass. The electric field \( E \) in Equations 4.13 and 4.15 is obtained from the Poisson equation A.3, while the diffusion coefficients ( \( D_{n/p} \)) are determined from the Einstein relation:

(4.16) \begin{equation} D_{n/p}=\frac {k_{\mathrm {B}}T_{\mathrm {L}}}{q}\mu _{n/p}=V_{\mathrm {T}}\mu _{n/p},\label {eq:Einstein's relation} \end{equation}

where \( V_{\mathrm {T}} \) is the thermal voltage which has a value of \( 26 \) \( \, \)mV ( \( k_{\mathrm {B}}=1.38\times 10^{-23} \) \( \, \)J/K, \( q=1.6\times 10^{-19} \) \( \, \)C) at room temperature ( \( 300 \) \( \, \)K). Equation 4.16 underestimates the diffusivity since the lattice temperature is used instead of the carrier temperature. Moreover, the numerous assumptions made during the derivation of the DD method limit its physical relevance for real device operation. However, many mobility models have been proposed to enhance the validity range of the DD approach by incorporating the electric field, doping dependencies and the role of the interface trap density to account for degradation effects [180, 181, 124]. Another major drawback of the DD method is that effects of rapidly changing electric fields, like the velocity overshoot, cannot be described [182]. In order to estimate the average carrier energy, the parameters obtained with the DD simulations are post-processed using the homogeneous energy balance equation as given in [28, 147]. However, in reality, the average energy lags behind the electric field if the latter changes rapidly. Since, to first order, the mobility depends (inversely) on the average energy, and not on the electric field, the mobility will not reduce instantly when the electric field increases. Thus, an overshoot in velocity will be observed until the carrier energy reached equilibrium with the electric field. Since the carrier energy depends directly on the electric field in the homogeneous energy balance equation, the velocity overshoot cannot be described. Nevertheless, due to its simplicity and low computational requirements, the DD model is the most widely used approach, especially for device engineering applications.

4.2.2.4 Hydrodynamic/Energy Transport Model

The full Hydrodynamic (HD) transport model involves numerically challenging equations due to their hyperbolic nature. Thus, the energy transport model, which is a simplified version of the HD transport model is often used. The energy transport employs the first three or four moments of the BTE, thereby, incorporating spatial dependencies of the average carrier energy. Using Equations 4.9, 4.10, and 4.11 and the corresponding closure as in [179], the HD transport model is given as:

(4.17–4.20) \{begin}{align} \nabla J_{n}= & q\left (R+\partial _{t}n\right ),\label {eq:HD_1}\\ J_{n}= & \mu _{\mathrm {n}}k_{\mathrm {B}}\left (\nabla \left (nT_{\mathrm {n}}\right )+\frac {q}{k_{\mathrm {B}}}En\right ),\label {eq:HD_2}\\ \nabla S_{n}= & -\frac {3}{2}k_{\mathrm {B}}\partial _{t}(nT_{\mathrm {n}})+EJ_{\mathrm {n}}-\frac {3}{2}k_{\mathrm {B}}n\frac {T_{\mathrm {n}}-T_{\mathrm {L}}}{\tau _{\mathrm {\varepsilon }}}+G\varepsilon n,\label {eq:HD_3}\\ S_{n}= & \frac {5}{2}\frac {k_{\mathrm {B}}}{q}T_{\mathrm {n}}J_{n}-\kappa _{\mathrm {n}}\nabla T_{\mathrm {n}}.\label {eq:HD_4} \{end}{align}

The electric field \( E \) is determined from the Poisson equation A.3, \( \tau _{\mathrm {\mathrm {\varepsilon }}} \) is the energy relaxation time, \( S\mathrm {_{n}} \) the heat flux density. The thermal conductivity \( \kappa _{\mathrm {n}} \) is obtained using the Wiedemann-Franz law, \( \kappa _{\mathrm {n}}T_{\mathrm {n}}=\left (\nicefrac {5}{2}-p'\right )\left (\nicefrac {k_{\mathrm {B}}}{q}\right )^{2}q\mu _{\mathrm {n}}nT_{\mathrm {n}} \) where \( p' \) is a correction factor. In contrast to the DD model, the energy transport model includes non-local effects. Thus, a phenomenon like the velocity overshoot can be approximately described since \( \mu _{\mathrm {n}} \) depends on \( T_{\mathrm {n}} \), which in turn depends on \( E \) in a non-local manner (Equations 4.19 and 4.20). However, the energy transport model predicts velocity overshoot also in rapidly decreasing fields, which is not observed in Monte-Carlo simulations. This artifact was attributed to the truncation of the model after the \( 4^{\mathrm {th}} \) moment [183, 184].

4.2.2.5 Six Moments Model

By considering two additional moments of the BTE along with those in the HD model, the six moments model can be derived [176]. An additional quantity, the kurtosis (skewness) of the distribution function is introduced as [174]:

(4.21) \begin{equation} \beta =\frac {3\left \langle \varepsilon ^{2}\right \rangle }{5\left \langle \varepsilon \right \rangle ^{2}}. \end{equation}

The quantity \( \beta \) represents the deviation of the shape of the distribution function from the Maxwellian distribution. This model better represents phenomena like impact ionization, velocity overshoot, hot-carrier effects, but it is quite complex and numerically challenging. Thus, the model has not been widely used in industry and has been restricted to academic interest.