Semiconductor Device Problems

Mathematics Subject Classification

82D37; 35Q20; 35Q40; 35Q79; 76Y05

Definition

Highly integrated electric circuits in computer processors mainly consist of semiconductor transistors which amplify and switch electronic signals. Roughly speaking, a semiconductor is a crystalline solid whose conductivity is intermediate between an insulator and a conductor. The modeling and simulation of semiconductor transistors and other devices is of paramount importance in the microelectronics industry to reduce the development cost and time. A semiconductor device problem is defined by the process of deriving physically accurate but computationally feasible model equations and of constructing efficient numerical algorithms for the solution of these equations. Depending on the device structure, size, and operating conditions, the main transport phenomena may be very different, caused by diffusion, drift, scattering, or quantum effects. This leads to a variety of model equations designed for a particular situation or a particular device. Furthermore, often not all available physical information is necessary, and simpler models are needed, helping to reduce the computational cost in the numerical simulation. One may distinguish four model classes: microscopic/mesoscopic and macroscopic semiclassical models and microscopic/mesoscopic and macroscopic quantum models (see Fig. 1).

Semiconductor Device Problems, Fig. 1

Hierarchy of some semiconductor models mentioned in the text

Description

In the following, we detail only some models from the four model classes since the field of semiconductor device problems became extremely large in recent years. For instance, we ignore compact models, hybrid model approaches, lattice heat equations, transport in subbands and magnetic fields, spintronics, and models for carbon nanotube, graphene, and polymer thin-film materials. For technological aspects, we refer to [9].

Microscopic Semiclassical Models

We are interested in the evolution of charge carriers moving in an electric field. Their motion can be modeled by Newton’s law. However, in view of the huge number of electrons involved, the solution of the Newton equations is computationally too expensive. Moreover, we are not interested in the trajectory of each single particle. Hence, a statistical approach seems to be sufficient, introducing the distribution function (or “probability density”) f(x, v, t) of an electron ensemble, depending on the position $x \in \mathbb{R}^{3}$, velocity $v =\dot{ x} = dx/dt \in \mathbb{R}^{3}$, and time t > 0. By Liouville’s theorem, the trajectory of f(x(t), v(t), t) does not change during time, in the absence of collisions, and hence,

$$\displaystyle{ 0 = \frac{df} {dt} = \partial _{t}f +\dot{ x} \cdot \nabla _{x}f +\dot{ v} \cdot \nabla _{v}f\quad \mbox{ along trajectories}, }$$

(1)

where $\partial _{t}f = \partial f/\partial t$ and $\nabla _{x}f$, $\nabla _{v}f$ are gradients with respect to x, v, respectively.

Since electrons are quantum particles (and position and velocity cannot be determined with arbitrary accuracy), we need to incorporate some quantum mechanics. As the solution of the many-particle Schrödinger equation in the whole space is out of reach, we need an approximate approach. First, by Bloch’s theorem, it is sufficient to solve the Schrödinger equation in a semiconductor lattice cell. Furthermore, the many-particle interactions are described by an effective Coulomb force. Finally, the properties of the semiconductor crystal are incorporated by the semiclassical Newton equations.

More precisely, let $p = \hslash k$ denote the crystal momentum, where $\hslash$ is the reduced Planck constant and k is the wave vector. For electrons with low energy, the velocity is proportional to the wave vector, $\dot{x} = \hslash k/m$, where m is the electron mass at rest. In the general case, we have to take into account the energy band structure of the semiconductor crystal (see [4, 7, 8] for details). Newton’s third law is formulated as $\dot{p} = q\nabla _{x}V$, where q is the elementary charge and V (x, t) is the electric potential. Then, using $\dot{v} =\dot{ p}/m$ and $\nabla _{k} = (m/\hslash )\nabla _{v}$, (1) becomes the (mesoscopic) Boltzmann transport equation:

$$\displaystyle\begin{array}{rcl} & & \partial _{t}f + \frac{\hslash } {m}k \cdot \nabla _{x}f + \frac{q} {\hslash }\nabla _{x}V \cdot \nabla _{k}f \\ & & = Q(f),\quad (x,k) \in \mathbb{R}^{3} \times \mathbb{R}^{3},\ t> 0,{}\end{array}$$

(2)

where Q(f) models collisions of electrons with phonons, impurities, or other particles. The moments of f are interpreted as the particle density n(x, t), current density J(x, t), and energy density (ne)(x, t):

$$\displaystyle\begin{array}{rcl} n =\int _{\mathbb{R}^{3}}fdk,\quad J& =& \frac{\hslash } {m}\int _{\mathbb{R}^{3}}kfdk, \\ ne = \frac{\hslash ^{2}} {2m}\int _{\mathbb{R}^{3}}\vert k\vert ^{2}fdk.& & {}\end{array}$$

(3)

In the self-consistent setting, the electric potential V is computed from the Poisson equation $\varepsilon _{s}\varDelta V = q(n - C(x))$, where ɛ _s is the semiconductor permittivity and C(x) models charged background ions (doping profile). Since n depends on the distribution function f, the Boltzmann–Poisson system is nonlinear.

The Boltzmann transport equation is defined over the six-dimensional phase space (plus time) whose high dimensionality makes its numerical solution a very challenging task. One approach is to employ the Monte Carlo method which consists in simulating a stochastic process. Drawbacks of the method are the stochastic nature and the huge computational cost. An alternative is the use of deterministic solvers, for example, expanding the distribution function with spherical harmonics [6].

Macroscopic Semiclassical Models

When collisions become dominant in the semiconductor domain, that is, the mean free path (the length which a particle travels between two consecutive collision events) is much smaller than the device size, a fluid dynamical approach may be appropriate. Macroscopic models are derived from (2) by multiplying the equation by certain weight functions, that is 1, k, and | k | ²∕2, and integrating over the wave-vector space. Setting all physical constants to one in the following, for notational simplicity, we obtain, using the definitions (3), the balance equations:

$$\displaystyle\begin{array}{rcl} \partial _{t}n +\mathop{ \mathrm{div}}\nolimits _{x}J =\int _{\mathbb{R}^{3}}Q(f)dk,\quad x \in \mathbb{R}^{3},\ t> 0,& &{}\end{array}$$

(4)

$$\displaystyle\begin{array}{rcl} \partial _{t}J +\mathop{ \mathrm{div}}\nolimits _{x}\int _{\mathbb{R}^{3}}k \otimes kfdk - n\nabla _{x}V =\int _{\mathbb{R}^{3}}kQ(f)dk,& &{}\end{array}$$

(5)

$$\displaystyle\begin{array}{rcl} & & \partial _{t}(ne) + \frac{1} {2}\mathop{\mathrm{div}}\nolimits _{x}\int _{\mathbb{R}^{3}}k\vert k\vert ^{2}fdk -\nabla _{ x}V \cdot \\ & &\qquad J = \frac{1} {2}\int _{\mathbb{R}^{3}}\vert k\vert ^{2}Q(f)dk. {}\end{array}$$

(6)

The higher-order integrals cannot be expressed in terms of the moments (3), which is called the closure problem. It can be solved by approximating f by the equilibrium distribution f₀, which can be justified by a scaling argument and asymptotic analysis. The equilibrium f₀ can be determined by maximizing the Boltzmann entropy under the constraints of given moments n, nu, and ne [4]. Inserting f₀ in (4)–(6) gives explicit expressions for the higher-order moments, yielding the so-called hydrodynamic model. Formally, there is some similarity with the Euler equations of fluid dynamics, and there has been an extensive discussion in the literature whether electron shock waves in semiconductors are realistic or not [10].

Diffusion models, which do not exhibit shock solutions, can be derived by a Chapman–Enskog expansion around the equilibrium distribution f₀ according to $f = f_{0} +\alpha f_{1}$, where α > 0 is the Knudsen number (the ratio of the mean free path and the device length) which is assumed to be small compared to one. The function f₁ turns out to be the solution of a certain operator equation involving the collision operator Q(f). Depending on the number of given moments, this leads to the drift-diffusion equations (particle density given):

$$\displaystyle\begin{array}{rcl} & & \partial _{t}n +\mathop{ \mathrm{div}}\nolimits _{x}J = 0,\quad J = -\nabla _{x}n + n\nabla _{x}V, \\ & & \quad \quad x \in \mathbb{R}^{3},\ t> 0, {}\end{array}$$

(7)

or the energy-transport equations (particle and energy densities given)

$$\displaystyle\begin{array}{rcl} & & \partial _{t}n +\mathop{ \mathrm{div}}\nolimits _{x}J = 0,\quad J = -\nabla _{x}n + \frac{n} {T}\nabla _{x}V, \\ & & \qquad x \in \mathbb{R}^{3},\ t> 0, {}\end{array}$$

(8)

$$\displaystyle\begin{array}{rcl} & & \partial _{t}(ne) +\mathop{ \mathrm{div}}\nolimits _{x}S + nu \cdot \nabla _{x}V = 0, \\ & & \qquad S = -\frac{3} {2}(\nabla _{x}(nT) - n\nabla _{x}V ), {}\end{array}$$

(9)

where $ne = \frac{3} {2}nT$, T being the electron temperature, and S is the heat flux. For the derivation of these models; we have assumed that the equilibrium distribution is given by Maxwell–Boltzmann statistics and that the elastic scattering rate is proportional to the wave vector. More general models can be derived too, see [4, Chap. 6].

The drift-diffusion model gives a good description of the transport in semiconductor devices close to equilibrium but it is not accurate enough for submicron devices due to, for example, temperature effects, which can be modeled by the energy-transport equations.

In the presence of high electric fields, the stationary equations corresponding to (7)–(9) are convection dominant. This can be handled by the Scharfetter–Gummel discretization technique. The key idea is to approximate the current density along each edge in a mesh by a constant, yielding an exponential approximation of the electric potential. This technique is related to mixed finite-element and finite-volume methods [2]. Another idea to eliminate the convective terms is to employ (dual) entropy variables. For instance, for the energy-transport equations, the dual entropy variables are $w = (w_{1},w_{2}) = ((\mu -V )/T,-1/T)$, where μ is the chemical potential, given by $n = T^{3/2}\exp (\mu /T)$. Then (8) and (9) can be formulated as the system:

$$\displaystyle{\partial _{t}b(w) -\mathop{\mathrm{div}}\nolimits (D(w,V )\nabla w) = 0,}$$

where $b(w) = (n, \frac{3} {2}nT)^{\top }$ and $D(w,V ) \in \mathbb{R}^{2\times 2}$ is a symmetric positive definite diffusion matrix [4] such that standard finite-element techniques are applicable.

Microscopic Quantum Models

The semiclassical approach is reasonable if the carriers can be treated as particles. The validity of this description is measured by the de Broglie wavelength $\lambda _{B}$ corresponding to a thermal average carrier. When the electric potential varies rapidly on the scale of $\lambda _{B}$ or when the mean free path is much larger than $\lambda _{B}$, quantum mechanical models are more appropriate. A general description is possible by the Liouville–von Neumann equation:

$$\displaystyle{i\varepsilon \partial _{t}\widehat{\rho } = [H,\widehat{\rho }] := H\widehat{\rho } -\widehat{\rho } H,\quad t> 0,}$$

for the density matrix operator $\widehat{\rho }$, where $i^{2} = -1$, ɛ > 0 is the scaled Planck constant, and H is the quantum mechanical Hamiltonian. The operator $\widehat{\rho }$ is assumed to possess a complete orthonormal set of eigenfunctions (ψ _j ) and eigenvalues $(\lambda _{j})$. The sequence of Schrödinger equations i ɛ ∂ _t ψ _j = H ψ _j ($j \in \mathbb{N}$), together with the numbers $\lambda _{j} \geq 0$, is called a mixed-state Schrödinger system with the particle density $n(x,t) =\sum _{ j=1}^{\infty }\lambda _{j}\vert \psi _{j}(x,t)\vert ^{2}$. In particular, $\lambda _{j}$ can be interpreted as the occupation probability of the state j.

The Schrödinger equation describes the evolution of a quantum state in an active region of a semiconductor device. It is used when inelastic scattering is sufficiently weak such that phase coherence can be assumed and effects such as resonant tunneling and quantum conductance can be observed. Typically, the device is connected to an exterior medium through access zones, which allows for the injection of charge carriers. Instead of solving the Schrödinger equation in the whole domain (self-consistently coupled to the Poisson equation), one wishes to solve the problem only in the active region and to prescribe transparent boundary conditions at the interfaces between the active and access zones. Such a situation is referred to as an open quantum system. The determination of transparent boundary conditions is a delicate issue since ad hoc approaches often lead to spurious oscillations which deteriorate the numerical solution [1].

Nonreversible interactions of the charge carriers with the environment can be modeled by the Lindblad equation:

$$\displaystyle{i\varepsilon \partial _{t}\widehat{\rho }=[H,\widehat{\rho }] + i\sum _{k}\Big(L_{k}\widehat{\rho }L_{k}^{{\ast}}-\frac{1} {2}(L_{k}^{{\ast}}L_{ k}\widehat{\rho } +\widehat{\rho } L_{k}^{{\ast}}L_{ k})\Big),}$$

where L _k are the so-called Lindblad operators and L _k ^∗ is the adjoint of L _k . In the Fourier picture, this equation can be formulated as a quantum kinetic equation, the (mesoscopic) Wigner–Boltzmann equation:

$$\displaystyle\begin{array}{rcl} \partial _{t}w + p \cdot \nabla _{x}w +\theta [V ]w = Q(w),& & \\ (x,p) \in \mathbb{R}^{3} \times \mathbb{R}^{3},\ t> 0,& &{}\end{array}$$

(10)

where p is the crystal momentum, $\theta [V ]w$ is the potential operator which is a nonlocal version of the drift term $\nabla _{x}V \cdot \nabla _{p}w$ [4, Chap. 11], and Q(w) is the collision operator. The Wigner function $w = W[\,\widehat{\rho }\,]$, where W denotes the Wigner transform, is essentially the Fourier transform of the density matrix. A nice feature of the Wigner equation is that it is a phase-space description, similar to the semiclassical Boltzmann equation. Its drawbacks are that the Wigner function cannot be interpreted as a probability density, as the Boltzmann distribution function, and that the Wigner equation has to be solved in the high dimensional phase space. A remedy is to derive macroscopic models which are discussed in the following section.

Macroscopic Quantum Models

Macroscopic models can be derived from the Wigner–Boltzmann equation (10) in a similar manner as from the Boltzmann equation (2). The main difference to the semiclassical approach is the definition of the equilibrium. Maximizing the von Neumann entropy under the constraints of given moments of a Wigner function w, the formal solution (if it exists) is given by the so-called quantum Maxwellian M[w], which is a nonlocal version of the semiclassical equilibrium. It was first suggested by Degond and Ringhofer and is related to the (unconstrained) quantum equilibrium given by Wigner in 1932 [3, 5]. We wish to derive moment equations from the Wigner–Boltzmann equation (10) for the particle density n, current density J, and energy density ne, defined by:

$$\displaystyle\begin{array}{rcl} n =\int _{\mathbb{R}^{3}}M[w]dp,\quad J =\int _{\mathbb{R}^{3}}pM[w]dp,& & {}\\ ne = \frac{1} {2}\int _{\mathbb{R}^{3}}\vert p\vert ^{2}M[w]dp.& & {}\\ \end{array}$$

Such a program was carried out by Degond et al. [3], using the simple relaxation-type operator $Q(w) = M[w] - w$. This leads to a hierarchy of quantum hydrodynamic and diffusion models which are, in contrast to their semiclassical counterparts, nonlocal.

When we employ only one moment (the particle density) and expand the resulting moment model in powers of ɛ up to order O(ɛ⁴) (to obtain local equations), we arrive at the quantum drift-diffusion (or density-gradient) equations:

$$\displaystyle\begin{array}{rcl} & & \partial _{t}n +\mathop{ \mathrm{div}}\nolimits _{x}J = 0,\quad J = -\nabla _{x}n + n\nabla _{x}V + \frac{\varepsilon ^{2}} {6}n\nabla _{x} {}\\ & & \quad \quad \quad \quad \quad \quad \quad \quad \quad \Big(\frac{\varDelta _{x}\sqrt{n}} {\sqrt{n}} \Big),\quad x \in \mathbb{R}^{3},\quad \ t>0. {}\\ \end{array}$$

This model is employed to simulate the carrier inversion layer near the oxide of a MOSFET (metal-oxide-semiconductor field-effect transistor). The main difficulty of the numerical discretization is the treatment of the highly nonlinear fourth-order quantum correction. However, there exist efficient exponentially fitted finite-element approximations, see the references of Pinnau in [4, Chap. 12].

Formally, the moment equations for the charge carriers and energy density give the quantum energy-transport model. Since its mathematical structure is less clear, we do not discuss this model [4, Chap. 13.2].

Employing all three moments n, nu, ne, the moment equations, expanded up to terms of order O(ɛ⁴), become the quantum hydrodynamic equations:

$$\displaystyle\begin{array}{rcl} & & \partial _{t}n +\mathop{ \mathrm{div}}\nolimits J = 0,\quad \partial _{t}J +\mathop{ \mathrm{div}}\nolimits _{x}\Big(\frac{J \otimes J} {n} + P\Big) {}\\ & & \qquad + n\nabla _{x}V = -\int _{\mathbb{R}^{3}}pQ(M[w])dp, {}\\ & & \partial _{t}(ne) -\mathop{\mathrm{div}}\nolimits _{x}((P + ne\mathbb{I})u - q) + \nabla _{x}V \cdot {}\\ & & J = \frac{1} {2}\int _{\mathbb{R}^{3}}\vert p\vert ^{2}Q(M[w])dp,\quad x \in \mathbb{R}^{3},\ t> 0, {}\\ \end{array}$$

where $\mathbb{I}$ is the identity matrix in $\mathbb{R}^{3\times 3}$, the quantum stress tensor P and the energy density ne are given by:

$$\displaystyle\begin{array}{rcl} P = nT\mathbb{I} - \frac{\varepsilon ^{2}} {12}n\nabla _{x}^{2}\log n,\quad ne = \frac{3} {2}nT& & {}\\ +\frac{1} {2}n\vert u\vert ^{2} - \frac{\varepsilon ^{2}} {24}n\varDelta _{x}\log n,& & {}\\ \end{array}$$

$u = J/n$ is the mean velocity, and $q = -(\varepsilon ^{2}/24)n(\varDelta _{x}u + 2\nabla _{x}\mathop{ \mathrm{div}}\nolimits _{x}u)$ is the quantum heat flux. When applying a Chapman–Enskog expansion around the quantum equilibrium, viscous effects are added, leading to quantum Navier–Stokes equations [5, Chap. 5]. These models are very interesting from a theoretical viewpoint since they exhibit a surprising nonlinear structure. Simulations of resonant tunneling diodes using these models give qualitatively reasonable results. However, as expected, quantum phenomena are easily destroyed by the occurring diffusive or viscous effects.