Subsections

• 4.1.1 Poisson's and Continuity Equation
• 4.1.2 Carrier Transport Equations
• 4.1.3 The Drift-Diffusion Model
• 4.1.4 Higher-Order Transport Models

4.1 Semiconductor Equations

4.1.1 Poisson's and Continuity Equation

In macroscopic semiconductor device modeling, Poisson's equation and the continuity equations play a fundamental role. Poisson's equation, one of the basic equations in electrostatics, is derived from the Maxwell's equation $\ensuremath{\mathbf{\nabla}}\cdot \ensuremath{\mathbf{D}}= \rho$ and the material relation $\ensuremath{\mathbf{D}}= \ensuremath{\ensuremath{{\mathrm{\hat{\epsilon}}}}}\ensuremath{\mathbf{E}}.$ $\ensuremath{\mathbf{D}}$ stands for the electric displacement field, $\ensuremath{\mathbf{E}}$ for the electric field, $\rho$ is the charge density, and $\ensuremath{\ensuremath{{\mathrm{\hat{\epsilon}}}}}$ is the permittivity tensor. Using the electrostatic potential $\psi$ with $\ensuremath{\mathbf{E}}= - \ensuremath{\mathbf{\nabla}}\psi$ leads to Poisson's equation

$$\ensuremath{\mathbf{\nabla}}\cdot ( \ensuremath{\ensuremath{{\mathrm{\hat{\epsilon}}}}}\ensuremath{\mathbf{\nabla}}\psi ) = - \rho .$$ (4.1)

Even for high frequencies, where the wavelengths are typically much smaller than the device dimension, the quasi stationary approximation used to justify Poisson's equation is still valid. Also the permittivity tensor $\ensuremath{\ensuremath{{\mathrm{\hat{\epsilon}}}}}$ is considered time invariant in this derivation. In isotropic materials like silicon the permittivity can be additionally approximated by the scalar value $\epsilon.$ Furthermore the permittivity is often considered to be constant within a material segment.

In semiconductors the charge density is commonly split into fixed charges which are in particular ionized acceptors $\ensuremath{N_{\mathrm{A}}^-}$ and donors $\ensuremath{N_{\mathrm{D}}^+}$ and into free charges which are electrons $n$ and holes $p.$ Including the acceptors, donors, electrons, and holes into (4.1), Poisson's equation can be written as

$$\ensuremath{\mathbf{\nabla}}^2 \psi = \frac{\ensuremath{\mathrm{q... ...silon} (n - p + \ensuremath{N_{\mathrm{A}}^-}- \ensuremath{N_{\mathrm{D}}^+}) .$$ (4.2)

The right hand side of (4.2) becomes zero in ideal oxides, and Poisson's equation reduces to the Laplace equation. However, the existence of oxide charges $N_\mathrm{ox}$ has to be considered, which yields

$$\ensuremath{\mathbf{\nabla}}^2 \psi = - \frac{\ensuremath{\mathrm{q}}}{\epsilon} N_\mathrm{ox}.$$ (4.3)

Additionally, at interfaces and surfaces, the interface trap concentration $N_\mathrm {it}$ has to be considered.

The continuity equation, can be also derived from Maxwell's equations and reads

$$\ensuremath{\mathbf{\nabla}}\cdot \ensuremath{\mathbf{J}}+ \ensuremath{\ensuremath{\frac{\partial \rho}{\partial t}}} = 0 .$$ (4.4)

The current density $\ensuremath{\mathbf{J}}$ is split into $\ensuremath{\mathbf{J}_{n}}$ and $\ensuremath{\mathbf{J}_{p}},$ for the contribution of electrons and holes, respectively. The electron and hole concentrations are written as $n$ and $p,$ respectively, and $\ensuremath{\mathrm{q}}$ is used to represent the elementary charge. By introducing the recombination rate $R,$ two separate continuity equations can be formulated as

$$\ensuremath{\mathbf{\nabla}}\cdot \ensuremath{\mathbf{J}_{n}}- \e... ...math{\ensuremath{\frac{\partial n}{\partial t}}} = + \ensuremath{\mathrm{q}}R ,$$ (4.5)

for electrons and

$$\ensuremath{\mathbf{\nabla}}\cdot \ensuremath{\mathbf{J}_{p}}+ \e... ...math{\ensuremath{\frac{\partial p}{\partial t}}} = - \ensuremath{\mathrm{q}}R ,$$ (4.6)

for holes. The separation into two equations allows independent transport modeling of the carrier types. The rate $R$ represents the net rate, which is zero in thermal equilibrium, following that generation and recombination are balanced. Generation and recombination rates of electrons and holes are expressed using physically or empirically based recombination models [11], which will be presented in Section 4.2. The equations (4.5) and (4.6) can alternatively be derived from Boltzmann's transport equation using the method of moments [15].

4.1.2 Carrier Transport Equations

A semiclassical description of carrier transport is given by Boltzmann's transport equation (BTE) which describes the evolution of the distribution function in the six-dimensional phase space ($x,$ $y,$ $z,$ $p_x,$ $p_y,$ $p_z$ ). Unfortunately, analytical solutions exist only for very simple configurations. One popular approach for solving the BTE in arbitrary structures is the Monte Carlo method [122] which gives highly accurate results. However, due to the statistical nature of the Monte Carlo method solutions are computationally very expensive. Due to the good agreement with experiments [123] results are often used as reference for the evaluation of simpler models. Alternatively, the spherical harmonics expansion (SHE) method as a deterministic numerical solution method of the BTE was presented already in the 1960s [124]. However, it has just been recently, that an efficiently use on real devices has been realized [125].

Device simulations on an engineering level require simpler transport equations which can be solved for complex structures within reasonable time. One method to perform this simplification is to consider only moments of the distribution function [126]. Depending on the number of moments considered in the model, different transport equations can be derived. The use of the first two moments, leads to the well known drift-diffusion model, a widely used approach for modeling carrier transport.

4.1.3 The Drift-Diffusion Model

Beside the derivation of the drift-diffusion by the method of moments [11], it is also possible using basic principles of irreversible thermodynamics [127]. The resulting electron and hole current relations contain at least two components caused by carrier drift and carrier diffusion. Additionally, the gradient of the lattice temperature ( $\ensuremath{\mathbf{\nabla}}T_{\mathrm{L}}$ ) acts as a driving force on the free carriers leading to [128]

$$\ensuremath{\mathbf{J}_{n}}= \ensuremath{\mathrm{q}}n \mu_n \ensu... ...}n + \ensuremath{\mathrm{q}}n D^T_n \ensuremath{\mathbf{\nabla}}T_{\mathrm{L}},$$ (4.7)

and

$$\ensuremath{\mathbf{J}_{p}}= \ensuremath{\mathrm{q}}p \mu_p \ensu... ...}p - \ensuremath{\mathrm{q}}p D^T_p \ensuremath{\mathbf{\nabla}}T_{\mathrm{L}}.$$ (4.8)

$\mu_n$ and $\mu_p$ represent the carrier electron and hole mobility, $D^T_n$ and $D^T_p$ the thermal diffusion coefficients, and $D_n$ and $D_p$ the diffusivity, respectively. The diffusivities are commonly expressed via the mobility invoking the Einstein relation

$$D_\nu = \mu_\nu \frac{\ensuremath{\mathrm{k_B}}T_{\mathrm{L}}}{\ensuremath{\mathrm{q}}} .$$ (4.9)

The index $\nu$ stands for $n$ and $p,$ respectively, while $\ensuremath{\mathrm{k_B}}$ is the Boltzmann constant. The Einstein relation is strictly valid only in equilibrium [16].

The equations (4.7) and (4.8) together with (4.5), (4.6), and (4.2) form the drift-diffusion model which was first presented by Van Roosbroeck in the year 1950 [129]. Rigorous derivations from the BTE show that many simplifications are required to obtain the drift-diffusion equations as will be shown. Simplifications include, for example, the assumption of a single parabolic band structure and the cold Maxwellian carrier distribution function. This assumes the carrier temperature equal to the lattice temperature. Nevertheless, due to its simplicity and its excellent numerical properties, the drift-diffusion equations have become the workhorse for most TCAD applications.

4.1.4 Higher-Order Transport Models

To obtain a better approximation of the BTE, higher-order transport models can be derived using more than just the first two moments [130]. The most prominent models beside the drift-diffusion model are the energy-transport/hydrodynamic models which use three or four moments. These models are based on the work of Stratton [131] and Bløtekjær [132]. A detailed review is given in [15]. In addition to the quantities used in the drift-diffusion model, the energy flux and the carrier temperatures are introduced as independent variables. Also new balance and flux equations are required, which introduce additional transport parameters. The carrier energy distribution function is here modeled using the heated Maxwellian distribution. Modeling of carrier mobility and impact-ionization benefit from more accurate models based on the carrier temperatures rather than the electric field. This advantage is caused by the non-local behavior of the average energy with respect to the electric field and becomes especially relevant for small device structures. This context is explained in Fig. 4.2.

Figure 4.2: Electron temperature along the channel of comparable n$^+$ -n-n$^+$ structures with varying channel lengths $L_\mathrm {C}$ (left).The spatial coordinates have been normalized to get an overlapping electric field. Therefore, the local electric field approach yields the same results for all device lengths. However, it can be seen that the local electric field approach is sufficient for larger structures only, but gives poor results for small structures. The distribution functions (right) are shown for $L_\mathrm {C} = 0.2$ µm at the positions A to D marked in the left figure. Note that the average energies at the points A and C are the same, whereas the distribution function looks completely different. Also note the high-energy tail at point D where the carrier temperature is already close to the lattice temperature. (Data taken from [133])

    

More detailed examinations in the far sub-micron area show that describing the energy distribution function using only the carrier concentration and the carrier temperature is still not sufficient for specific problems which depend on high energy tails (see Fig. 4.2). Hot carrier modeling in small structures, for example, which is based on accurate modeling of the high-energy part of the distribution function would require more complex models. Also the reliability modeling benefits of the detailed knowledge of the distribution function (more on this is highlighted in Chapter 6). One method for improving the approximation of the distribution function is the six moments method [134]. It introduces the kurtosis, which is the deviation of the distribution function from the heated Maxwellian.

It is evident that higher-order transport models give a closer solution of the BTE and therefore lead to a better agreement between simulation results and measurements of real devices. This effect is especially relevant for small structures, where non-local effects gain importance (see Fig. 4.2). Simulation results with drift-diffusion in deep sub-micrometer structures therefore seem to be very questionable [135].

The high-voltage devices considered in this work are relatively large. Hence, a rather high number of mesh points is required for a proper discretization. In combination with more elaborative transport equations, this leads to a higher computational time. Additionally, the convergence properties degrade significantly for higher moments models [136]. However, the relatively large dimensions of the high-voltage devices justify the use of the drift-diffusion model in this work.

In case advanced transport models have to be solved in complex devices, it is also possible to combine different transport equations in one simulation. To accomplish this, the semiconductor domain must be partitioned into separate segments. One segment must contain the critical areas, e.g. the channel area. Hence, the higher-order transport equations are solved for this segment, and the drift-diffusion model is used for the remaining ones. Here, it is essential to define proper boundary conditions between the segments [137].