Subsections

• 4.6.1 Analytical High Field Velocity Model
• 4.6.2 Parallel Velocity Model
• 4.6.3 Perpendicular Velocity Model
• 4.6.4 Total Velocity for Arbitrary Field Direction
• 4.6.5 Simplification to Two-Dimensional Simulation Domains
• 4.6.6 Implementation Issues
• 4.6.7 Discussion

4.6 High-Field Mobility

Figure 4.8: Velocity saturation in Si.

With an increase in the applied lateral field, carriers gain energies above the ambient thermal energy and are able to transfer energy to the lattice by optical phonon emission. This effect leads to a saturation of the carrier velocity as shown in Fig. 4.8. To account for this velocity saturation, the mobility has to be reduced accordingly. Popular models including this basic feature are the Caughey and Thomas expression [Caughey67],

$$\mu = \frac{\mu_0}{\left[1 + \left(\displaystyle \frac{\mu_0 E}{v_s}\right)^\beta \right]^{1/\beta} }$$ (4.104)

and its modification [IuE04],

$$\displaystyle \mu_{E} = \frac{2 \mu_0 }{1+\left[1+\left(\displaystyle \frac{2\mu_0 E }{ \displaystyle v_{s} }\right)^{\beta}\right]^{1/\beta} }.$$ (4.105)

Here $\mu_0$ denotes the low field mobility and $v_{s}$ the saturation velocity. The parameter $\beta$ influences the transition from low to high fields. Although (4.104) and (4.105) are applicable for a uniform field, they can be used for device simulations where the fields are non-uniform by making an appropriate choice of the lateral field. Different alternatives are available such as (a) the magnitude of the gradient of the quasi-fermi potential, (b) the component of electric field in the direction of current flow, and (c) the component of electric field parallel to the interface. Any of the definitions can be considered for device simulation purposes.

4.6.1 Analytical High Field Velocity Model

In general, the high field mobility is modeled differently for the drift-diffusion and the hydrodynamic transport model. In the former case mobility is modeled as a function of the driving force, whereas in the latter case a dependence on the carrier temperature is usually assumed. To describe nonlocal transport effects occurring in aggressively scaled devices a mobility model for the hydrodynamic framework would be desirable. Such a model could include a three-valley band structure and deal with arbitrary strain conditions. It would capture the essential physics of multi-valley transport under a spatially rapidly varying electric field profile. However, one problem is complexity. A nonlinear system of nine unknowns, namely the valley populations, valley velocities and valley temperatures, has to be solved numerically. The strain and carrier temperature dependence for each valley would require careful modeling in order to obtain realistic carrier temperatures and, consequently, realistic valley population. In the past, multi-valley transport models have also been devised for compound semiconductors [Hänsch91]. As a matter of fact, it seems that such multi-valley transport models with separate carrier gases for each valley have never found application in commercial or academic device simulators.

To find a trade-off between physical rigor and an acceptable level of model complexity we abstained from the multi-valley approach and pursued a more empirical approach, where analytical expressions for the velocity-field characteristics are directly fitted to bulk Monte Carlo data. The model is restricted to such strain conditions where only one pair of X-valleys is shifted and four valleys remain degenerate. These conditions include biaxial stress and uniaxial stress applied along the $\langle$100$\rangle$ axes of Si. Another condition resulting in a separation of the $\Delta _2$ and $\Delta_4$ valleys is uniaxial stress in the [110] direction. Equation (3.42) shows that the valley splitting depends only on the diagonal elements of the strain tensor. The proposed mobility model is thus applicable, if two diagonal elements are equal, $\epsilon_{11} = \epsilon_{22} \ne \epsilon_{33}$.

To develop a clear understanding of the model, we have to consider three different coordinate systems.

1. The principal coordinate system has to be oriented such that the unit vectors ${\mathbf{e}}_1$, ${\mathbf{e}}_2$, and ${\mathbf{e}}_3$ correspond to the [100], [010], and [001] crystallographic directions, respectively. In this system the $\Delta_4$-valleys are aligned along the [100] and [010] directions, whereas the $\Delta _2$- valleys are aligned along the [001] direction.

2. The unit vectors ${\mathbf{e}}_x$, ${\mathbf{e}}_y$, and ${\mathbf{e}}_z$ constitute the device coordinate system. In this system the device geometry is defined. For performing device simulations it is essential to transform all transport parameters into this coordinate system.

3. A polar coordinate system is employed, comprising a unit vector along the field direction, ${\mathbf{e}}_E = {\mathbf{E}}/\vert{E}\vert$, and two orthogonal vectors ${\mathbf{e}}_\theta$ and ${\mathbf{e}}_{\varphi}$. The polar axis is aligned with the [001] direction. In terms of the polar angle $\theta$ and the in-plane (azimuth) angle $\varphi$, the unit vectors are defined as follows.

   $${\mathbf{e}}_E = \begin{pmatrix}\displaystyle \sin(\theta)\cos(\v... ...phi = \frac{1}{\sin(\theta)}\frac{\partial{{\mathbf{e}}_E}} {\partial{\varphi}}$$ (4.106)

   
   

4.6.2 Parallel Velocity Model

Although (4.105) can describe the high field behavior in unstrained Si, it cannot produce the peculiarities of the velocity field relation in the strained case. Here the following expression has been suggested which can handle all types of velocity-field characteristics resulting from the Monte Carlo simulations performed [Dhar05a,Dhar06].

$$\displaystyle v_{E} = \frac{2 \mu_{EE} E}{1+\left[1+\left(\displa... ...1/\beta} } + \displaystyle v_{s}\xi\frac{(E/ \eta)^\gamma}{1+(E / \eta)^\gamma}$$ (4.107)

Here $\mu_{EE}$ denotes the low-field mobility in the field direction, obtained by projection of the low-field mobility tensor as $\mu_{EE} = {\mathbf{e}}_E^T\cdot \mu_0\cdot{\mathbf{e}}_E$. The last term incorporated in (4.107) models the velocity kink shown in Fig. 5.12. The relevance of the parameter $\xi$ is twofold: It accounts for the velocity plateau occurring approximately at $v_s(1-\xi)$ and also signifies the small negative differential mobility occurring in strained Si for higher strain levels. The parameters $\eta$ and $\gamma$ are fitting parameters.

All parameters depend on the strain-induced valley splitting, $\Delta\varepsilon=\varepsilon(\Delta_2)-\varepsilon(\Delta_4)$. The following empirical expressions were assumed.

$$v_s = v_{s1} + v_{s2} \cdot\Delta{\varepsilon}$$ (4.108)

$$\beta = \beta_{1} + \beta_{2} \cdot \Delta{\varepsilon}$$ (4.109)

$$\eta = \eta_{1} + \eta_{2}\cdot\Delta{\varepsilon}$$ (4.110)

$$\gamma = \gamma_{1} + \gamma_{2}\cdot\Delta{\varepsilon}$$ (4.111)

$$\xi = \frac{(\Delta{\varepsilon}/\xi_{1})} {1+(\Delta{\varepsilon}\cdot \xi_{2}/\xi_{1})^2}$$ (4.112)

For all parameters except $\xi$, a linear dependence was found to be sufficient. The parameter $\xi$ was modeled by the rational expression in (4.112). The parameters $v_{si}$, $\beta_{i}$, $\eta_{i}$, $\gamma_{i}$, $\xi_{i}$ where $i=1,2$, are constants for a particular field direction. We have chosen the three high symmetry directions [100], [110], and [001] and two additional directions [101] and [11$\sqrt{2}$]. These five sample directions form a spherical triangle on a unit sphere as shown in Fig. 4.9.

Figure 4.9: Spherical triangle with the five field directions chosen.

4.6.3 Perpendicular Velocity Model

For the cases where the field is not oriented along a high symmetry direction, it is observed that an electron velocity component perpendicular to the field direction develops. The component ${\mathbf{v}}_\theta$, although small for low stress levels, has a significant magnitude for intermediate field regimes and results in a total velocity different from the parallel velocity. For symmetry reasons, the velocity component in the ${\mathbf{e}}_\varphi$ direction vanishes for all five sample directions.

The perpendicular velocity component vanishes for fields along the [100], [110] and [001] directions. For the field directions [101] and [11$\sqrt{2}$], the normal velocity can be expressed in terms of $v_E$ and $v_3$.

$$v_{\perp} = v_{E} - \sqrt{2}v_{3}$$ (4.113)

After fitting $v_E$, the component $v_3$ is fitted using an expression similar to (4.107)

$$\displaystyle v_{3} = \frac{\sqrt{2} \mu_{33} E}{1+\left[1+\left(... ...laystyle \frac{v_{s}}{\sqrt{2}}\xi\frac{(E / \eta)^\gamma}{1+(E / \eta)^\gamma}$$ (4.114)

To ensure the correct low-field behavior, $v_3 = \mu_{33}E_3$, the magnitude of the electric field $E$ in the first term in (4.107) has to be replaced by $E_3 = E/\sqrt{2}$ to obtain (4.114). The correct high-field limit is introduced by replacing $v_s$ in (4.107) by $v_{3,s} = v_s/\sqrt{2}$.

The total electron velocity vector for a fixed field direction is obtained by addition of the two components.

$${\mathbf{v}}_t = v_E{\mathbf{e}}_E + v_{\theta}{\mathbf{e}}_{\theta}$$ (4.115)

where ${\mathbf{e}}_{E}$ and ${\mathbf{e}}_{\theta}$ are the unit vectors parallel and perpendicular to the field direction.

4.6.4 Total Velocity for Arbitrary Field Direction

The velocity-field characteristics can be extended to other field directions using a spherical harmonics interpolation.

$$\Omega(\theta,\varphi) = \sum_{l=0}^{\infty}\sum_{m=0}^{l}a_{lm}P_l^m[\cos(\theta)]\cos(m\varphi)$$ (4.116)

Here, $\Omega$ is the function to be interpolated, $a_{lm}$ denote the expansion coefficients and $P^m_l$ are the associated Legendre polynomials. From the symmetry properties $\Omega(\theta,\varphi+\pi/2) = \Omega(\theta,\varphi)$ and $\Omega(\theta+\pi,\varphi) = \Omega(\theta,\varphi)$ it follows that $l$ must be even and $m=4n$. Truncating (4.116) after the 4$^{th}$ order yields

$$\Omega(\theta,\varphi) = a_{00}P_0^0(\chi) + a_{20}P_2^0(\chi)+ a_{40}P_4^0(\chi) + a_{44}P_4^4(\chi)\cos(4\varphi)$$ (4.117)

where $\chi=\cos(\theta)$. Evaluating (4.117) for the field directions [100], [110], [001], [101] and [11$\sqrt{2}$] gives

$$\Omega_{100} =a_{00} - \frac{1}{2}a_{20} + \frac{3}{8}a_{40} + 105a_{44}$$ (4.118)

$$\Omega_{110} =a_{00} - \frac{1}{2}a_{20} + \frac{3}{8}a_{40} - 105a_{44}$$ (4.119)

$$\Omega_{001} =a_{00} + a_{20} + a_{40}$$ (4.120)

$$\Omega_{011} =a_{00} + \frac{1}{4}a_{20} - \frac{13}{32}a_{40} + \frac{105}{4}a_{44}$$ (4.121)

$$\Omega_{11\sqrt{2}} =a_{00} + \frac{1}{4}a_{20} - \frac{13}{32}a_{40} - \frac{105}{4}a_{44}$$ (4.122)

To determine the four coefficients from the overdetermined system (4.118) to (4.122), we solve (4.118) to (4.120) exactly and minimize the error in (4.121) and (4.122), as shown in Appendix B.

$$a_{00} = \frac{1}{3}(\Omega_{100}+\Omega_{110}+\Omega_{001})-\frac{7}{12}a_{40},$$ (4.123)

$$a_{20} =\frac{1}{3}(2\Omega_{001}-\Omega_{100}-\Omega_{110})-\frac{5}{12}a_{40},$$ (4.124)

$$a_{40} =\frac{8}{35}(\Omega_{100}+\Omega_{110}+2\Omega_{001}-2\Omega_{011}-2\Omega_{11\sqrt{2}}),$$ (4.125)

$$a_{44} =\frac{1}{210}(\Omega_{100}-\Omega_{110}).$$ (4.126)

4.6.5 Simplification to Two-Dimensional Simulation Domains

For two-dimensional simulation domains a simpler interpolation method can be used. In the following we consider the case that ${\mathbf{e}}_3$ lies in the simulation domain, such that $\varphi$ is fixed and only the variation of $\theta$ has to be considered. An alternative case would be a uniaxial stress orthogonal to the simulation domain, for example, in the width direction of a MOSFET. Then $\theta=\pi/2$ and $\varphi$ varies.

In the first case the quantity $\Omega$ can be interpolated using the following polynomial.

$$\Omega(\theta) = b_0 + b_2\cos^2(\theta) + b_4\cos^4(\theta)$$ (4.127)

Considering the special case of transport in the ($010$) plane ($\varphi=0$), we can write the equation system

$$\Omega_{100} =b_0$$ (4.128)

$$\Omega_{101} =b_0+ \frac{b_2}{2} + \frac{b_4}{4}$$ (4.129)

$$\Omega_{001} =b_0+ b_2 + b_4$$ (4.130)

which gives the coefficients

$$b_2 =-3\Omega_{100} + 4\Omega_{101} - \Omega_{001}$$ (4.131)

$$b_4 =2\Omega_{100} - 4\Omega_{101} + 2\Omega_{001}$$ (4.132)

Similarly, for transport in a ( $\overline{1}10$) plane ( $\varphi=\frac{\pi}{4}$), we have

$$\Omega_{110} =b_0$$ (4.133)

$$\Omega_{11\sqrt{2}} =b_0+ \frac{b_2}{2} + \frac{b_4}{4}$$ (4.134)

$$\Omega_{001} =b_0+ b_2 + b_4$$ (4.135)

giving

$$b_2 =-3\Omega_{110} + 4\Omega_{11\sqrt{2}} - \Omega_{001}$$ (4.136)

$$b_4 =2\Omega_{110} - 4\Omega_{11\sqrt{2}} + 2\Omega_{001}$$ (4.137)

4.6.6 Implementation Issues

The present model has been derived for a uniform electric field ${\mathbf{E}}$. To apply it in a drift diffusion based device simulator, the electric field in the model has to be replaced by an appropriately defined driving force ${\mathbf{F}}_n$. Typical definitions of the driving force employed in practical devices simulators are the electric field component along the current density vector or the gradient of the quasi Fermi level.

For the two-dimensional cases described in Section 4.6.5, only one angle has to be determined. With the constant, two-dimensional vector ${\mathbf{e}}_3$ denoting the [001] direction, one obtains for the polar angle

$$\cos^2(\theta) = \frac{ ({\mathbf{F}}_n\cdot{\mathbf{e}}_3)^2 } { {\mathbf{F}}_n\cdot{\mathbf{F}}_n }.$$ (4.138)

The involved vectors are two-dimensional and specified in the device coordinate system.

The task of finding a tensor $\ensuremath{{\underline{\mu}}}$ that relates two given vectors ${\mathbf{F}}_n$ and ${\mathbf{v}}$ by ${\mathbf{v}} = \ensuremath{{\underline{\mu}}}{\mathbf{F}}_n$ has no unique solution. In two dimensions, this vector equation denotes two scalar equations, whereas the 2 by 2 matrix $\ensuremath{{\underline{\mu}}}$ has four unknown elements. The straightforward assumption of a diagonal matrix with diagonal elements $\mu_{xx} = v_x / F_{n,x}$ and $\mu_{yy} = v_y / F_{n,y}$ cannot be made, because this definition becomes singular whenever the driving force vector is parallel to the $x$ or $y$ axis.

To solve this problem one can probe the velocity vector for a second driving force vector. Using Cartesian coordinates, a vector orthogonal to ${\mathbf{F}}_n = (F_{n,x}, F_{n,y})$ can be easily found as ${\mathbf{G}}_n = (-F_{n,y}, F_{n,x})$. The four matrix elements can now be uniquely determined from the following two vector equations.

$${\mathbf{v}}_t({\mathbf{F}}_n) =\ensuremath{{\underline{\mu}}} \hspace*{1mm} {\mathbf{F}}_n$$ (4.139)

$${\mathbf{v}}_t({\mathbf{G}}_n) =\ensuremath{{\underline{\mu}}} \hspace*{1mm} {\mathbf{G}}_n$$ (4.140)

Implementation of a full tensorial mobility model in two dimensions would require the following quantities. We consider some edge $\alpha$ of an unstructured mesh, connecting the nodes $i$ and $j$. The unit vector along this edge is given by ${\mathbf{e}}_\alpha = ({\mathbf{r}}_j - {\mathbf{r}}_i) / \vert{\mathbf{r}}_j - {\mathbf{r}}_i\vert$. A vector ${\mathbf{e}}_\beta$ orthogonal to ${\mathbf{e}}_\alpha = (\alpha_x, \alpha_y)$ can be defined as ${\mathbf{e}}_\beta = (-\alpha_y, \alpha_x)$. The box integration method requires the projection of the current density onto the edge, $j_\alpha = {\mathbf{j}} \cdot {\mathbf{e}}_\alpha$. For a mobility tensor with non-zero off-diagonal elements this current component becomes

$$j_\alpha = q n_\alpha (\mu_{\alpha\alpha} F_\alpha + \mu_{\alpha\beta} F_\beta),$$ (4.141)

where $n_\alpha$ denotes the carrier concentration at the mid point edge $\alpha$. This equation states that the current component along the edge is driven not only by the driving force component along the edge, $F_\alpha$, but also by the perpendicular component $F_\beta$. The Scharfetter-Gummel scheme [Scharfetter69] gives a discrete representation of the component $F_\alpha$ as a function of the variables at nodes $i$ and $j$. The perpendicular component $F_\beta$, however, cannot be determined from the variables at the two nodes. It can only be estimated by some kind of interpolation of the parallel components at neighboring edges. One possible extension of the Scharfetter-Gummel discretization has been proposed in [Egley93].

4.6.7 Discussion

The high-field velocity model presented has been derived for bulk Si. The model can consistently be used with the drift-diffusion transport model whenever the spatial variations of the potential are sufficiently smooth. However, even nowadays where gate lengths are in the deca-nanometer regime, drift-diffusion based simulations are still widely used to estimate transistor performance, despite the fact that this model cannot capture the strong non-local transport effects. The latter are usually accommodated by changing parameters, in particular by considerably increasing the saturation velocity. Also the high-field velocity model presented cannot include non-local effects, as it is based on populations of the band minima being functions of the local electric field. One should be aware that extending the local drift-diffusion equation with a local high-field velocity model cannot extend the limitation with respect to strongly non-local transport occurring in relevant technologies. Up-scaling the velocity by some technology-dependent parameter may still add some physics to the model, as a rough estimate of the direction-dependence of high-field transport.