Subsections

• 7.2.1 Derivation
• 7.2.2 Discretization of Edges
• 7.2.3 Discretization of the Right Hand Side
• 7.2.4 Limitations of the Box Discretization Method

7.2 Box Discretization

The discretization of the partial differential semiconductor equations in space and time is required to obtain difference equations which can be solved using numerical methods. A common approach for the discretization of the differential equations is the box discretization or box integration method [259,12], also known as the finite volume method.

7.2.1 Derivation

To apply the box discretization method, the mesh has to fulfill the Delaunay criterion and can therefore be split into boxes using a Voronoi tessellation [256]. By applying the Voronoi tessellation a dual mesh is created so that every point in the domain is assigned to its closest vertex. All points which are closest to the same vertex form together an element of the dual mesh which is called box volume throughout this work (see Fig. 7.4). The transformation is bijective.

The box discretization method reformulates the equations based on this Voronoi tessellation and can be equally used in two- and three-dimensional simulation domains. Since the two-dimensional devices have a specified width which acts as a multiplier, all boxes are considered as volumes and the box boundary elements are considered as surfaces. The basic method of the box discretization concerns the divergence operator in the form

$$\ensuremath{\mathbf{\nabla}}\cdot \ensuremath{\mathbf{F}}= g ,$$ (7.1)

where $\ensuremath{\mathbf{F}}$ and $g$ are an arbitrary flux and generation term, respectively. The equation can be transformed by integration over a volume $V$

$$\int_V \ensuremath{\mathbf{\nabla}}\cdot \ensuremath{\mathbf{F}}\ensuremath{ \mathrm{d}}V = \int_V g \ensuremath{ \mathrm{d}}V,$$ (7.2)

and by applying Gauss's theorem,

$$\oint\limits_{\partial V} \ensuremath{\mathbf{F}}\ensuremath{ \mathrm{d}}\ensuremath{\mathbf{A}}= \int_V g \ensuremath{ \mathrm{d}}V,$$ (7.3)

where ${\partial V}$ is the surface of V and $\ensuremath{\mathbf{A}}$ is the outwardly oriented surface area vector. In the next step equation (7.3) is applied to each Voronoi box in the mesh (see Fig. 7.5). With this spatial quantization in boxes the two-dimensional boundary $\partial V$ is transferred to a polygon (in 2D) or polyhedron (in 3D) and can be split into planar surface elements perpendicular to the edges leading to the neighboring vertex. Every vertex $i$ has an associated box volume $V_i$ and a set of $N_i$ neighboring boxes. Equation (7.3) can be approximated, i.e. discretized, for each vertex as

$$\sum_{j \in N_i} F_{ij} \ensuremath{A_{ij}}= \int_{V_i} g \ensuremath{ \mathrm{d}}V .$$ (7.4)

$F_{ij}$ is the flux through the boundary area $\ensuremath {A_{ij}}$ between the boxes $i$ and $j$ (compare Fig. 7.5) and can be interpreted as the projection of $\ensuremath{\mathbf{F}}$ on the edge between the vertices $i$ and $j$ as

$$F_{ij} = \bm{\xi}_{ij} \cdot \ensuremath{\mathbf{F}}.$$ (7.5)

$\bm{\xi}_{ij}$ is the unit vector pointing along edge $ij.$ From the discretization point of view, $F_{ij}$ must be the value of the quantity at the mid-point of the edge which is exactly the boundary between the boxes. As will be shown in the next section, the proper discretization of $F_{ij}$ is relevant for numerically stable results and depends on the type of differential equation system used.

Figure 7.5: Voronoi box $i$ (blue, dashed) of vertex $i$ with connections to its neighboring vertices forming six adjacent triangular elements. The flux $F_{ij}$ from box $i$ to box $j$ through the area $\ensuremath {A_{ij}}$ is depicted. The edge length between $i$ and $j$ is $d_{ij}.$

The physical equations that describe semiconductor devices are laws of conservation. The derivation of the box discretization method is inherently conservative [260]. The box discretization scheme is, therefore, widely used in semiconductor device simulations [12]. In drift-diffusion simulations, the differential equations that have to be solved are Poisson's equation (4.1) and the current continuity equations (4.5) and (4.6). The flux quantities $\ensuremath{\mathbf{F}}$ introduced in (7.1) are the dielectric flux density $\ensuremath{\mathbf{D}}= \epsilon \ensuremath{\mathbf{E}}$ for Poisson's equation and the current densities $\ensuremath{\mathbf{J}}_\nu$ for the continuity equation, respectively. The generation terms $g$ in (7.1) are the charge density $\rho$ and the carrier generation/recombination rates $R$ /$G,$ respectively.

A considerable advantage of the box discretization method is that the only information required is the unstructured neighborhood information. This neighborhood information consists basically of two lists. First, a list of all boxes in the simulation domain, together with their center points ($i$ ), i.e. the vertices, and the box volumes ($V_i$ ). Second, a list of all edges connecting the vertices, together with the edge length ($d_{ij}$ ), i.e. the distance between the vertices, and the common surface between the neighboring boxes ( $\ensuremath {A_{ij}}$ ). It has to be noted that no more information about the elements is necessary for the evaluation of (7.4), which makes the box discretization independent of the box shape, including its dimensionality. This makes this scheme very straightforward to implement and a dimensional independent set-up of the simulator can be realized.

7.2.2 Discretization of Edges

The calculation of the one-dimensional projection of the flux on the edge depends on the type of the differential equation. For the dielectric flux density $\ensuremath{\mathbf{D}}$ in Poisson's equation the flux from box $i$ to box $j$ along the connecting edge through the common boundary area $\ensuremath {A_{ij}}$ is $D_{ij}.$ This density is commonly approximated using the directional derivative of the electrostatic potential:

$$D_{ij} = \bm{\xi}_{ij} \cdot \ensuremath{\mathbf{D}}= - \epsilon ... ...\partial \bm{\xi}_{ij}}} \approx - \epsilon \frac{\psi_j - \psi_i}{d_{ij}} .$$ (7.6)

$d_{ij}$ is the distance between the two vertices and the permittivity $\epsilon$ is considered as a constant in each element. In MINIMOS-NT, for example, the permittivity is taken as the average between the permittivities associated with the vertices $i$ and $j$ as $\epsilon _{ij} = \left.\left( \epsilon_i + \epsilon_j \right)\middle/2\right..$

For the current densities $\ensuremath{\mathbf{J}}_{\nu}$ of the carrier type $\nu,$ the discretization of $J_{\nu,ij}$ is not as straightforward. Insertion of the drift-diffusion current relations from equations (4.7) and (4.8) in the continuity equations (4.5) and (4.6) results in a second order parabolic partial differential equation. Using a simple finite difference approach like in (7.6) leads to numerical oscillations if the drift term dominates over the diffusion term [261]. Very fine meshes would be necessary to stabilize the system. A stable discretization can be obtained using the Scharfetter-Gummel method [39] instead. Here, the drift-diffusion current equations (4.7) and (4.8) are used to solve the one-dimensional carrier concentrations, $n$ and $p,$ respectively, along the edge. The boundary conditions of the carrier concentrations are given using the values at the corresponding vertices $i$ and $j.$ The values of $J_{n,ij}$ and $J_{p,ij},$ $E_{ij},$ and $\mu_n$ and $\mu_p$ are considered constant along the edge. Solving this one-dimensional differential equation results in

$$J_{n,ij} = \frac{\ensuremath{\mathrm{q}}\mu_n V_\mathrm{T}}{d_{ij... ...g( n_j {\mathscr{B}}(\triangle_{ij}) - n_i {\mathscr{B}}(-\triangle_{ij}) \Big)$$ (7.7)

for electrons and

$$J_{p,ij} = - \frac{\ensuremath{\mathrm{q}}\mu_p V_\mathrm{T}}{d_{... ...g( p_j {\mathscr{B}}(-\triangle_{ij}) - p_i {\mathscr{B}}(\triangle_{ij}) \Big)$$ (7.8)

for holes, where the temperature voltage is defined as $V_\mathrm{T}= \left.\ensuremath{\mathrm{k_B}} T\middle/\ensuremath{\mathrm{q}}\right.,$ $\triangle_{ij} = \left.\left( \psi_j-\psi_i \right)\middle/V_\mathrm{T}\right.$ and ${\mathscr{B}}$ is the Bernoulli function defined as

$${\mathscr{B}}(x) = \frac{x}{\exp (x) - 1} .$$ (7.9)

7.2.3 Discretization of the Right Hand Side

Up to now, the generation term $g$ on the right hand side of (7.4), which is the charge density $\rho$ in Poisson's equation and the carrier generation and recombination rate in the carrier equations, was only represented as an integral in continuum space. Most implementations calculate this integral by partitioning the box into pieces and adding the contributions to the integral of the box. Three different methods are considered in this work, each partitioning the Voronoi box differently. For the first assumption, the box is not split and the generation is considered constant over the whole box volume, as it is depicted in Fig. 7.6(a). This approach reduces the integral term of (7.4) to a simple product [12,262] which reads

$$\int\limits_{V_i} g \ensuremath{ \mathrm{d}}V \approx g_i V_i .$$ (7.10)

Figure 7.6: The Voronoi box $i$ (blue) and the different approaches to decompose the generation for discretization. Areas of constant generation are surrounded in bold red (dashed plus solid), parts contributing to the box are shown in solid red, for (a) the box based methods, (b) the element based methods, and (c) the method proposed by Laux et al. [263].

This formulation perfectly fits to the box discretization method, because it only requires quantities which are stored on vertices and the only geometry information needed is the unstructured neighborhood information. This makes the assembling procedure of (7.10) very simple. MINIMOS-NT [120], for example, is one simulation tool using this technique. As will be discussed in Section 7.3, the calculation of vector quantities is somewhat more involved in this approach than for the other methods.

Another approach for assembling the generation integral is to assume the generation rate to be constant within one mesh element. In this case, the integral for one box volume is assembled using contributions from all adjacent elements $e.$ Using the naming conventions from Fig. 7.6(b), this reads

$$\int\limits_{V_i} g \ensuremath{ \mathrm{d}}V \approx \sum_{e \in \mathrm{E}_i} g_{e,i} V_{e,i} ,$$ (7.11)

where $\mathrm{E}_i$ is the set of all elements that contribute to the volume $V_i.$ In this scheme the generation term is estimated for each element (compare dashed red line in Fig. 7.6(b)), which might be a triangle, tetrahedra, box, or other mesh element. Its contributions can be used for all boxes the element is part of [264]. This element-wise assembly results in an implicitly increased resolution and a more accurate physical representation [265]. Some implementations associate different current densities for each edge pair (compare solid red line in Fig. 7.6(b)) and do not assume it to be constant within the whole element. In this case there are independent rates for each box the element is contained in. A method using this techniques is the edge-pair method [266].

The third approach in this classification is depicted in Fig. 7.6(c). This approach was presented for triangles only [263], but extensions to tetrahedra are possible. Here, each triangle is split into three different regions associated with the three edges. In each of them, a generation term is estimated, and two of them contribute to the sum of one box. Therefore, the summation for the generation integral of the box requires to consider two contributions per element:

$$\int\limits_{V_i} g \ensuremath{ \mathrm{d}}V \approx \sum_{e \in \mathrm{E}_i} \left( g_{e,ij} V_{e,ij} + g_{e,ik} V_{e,ik} \right) .$$ (7.12)

This approach was used by Laux et al. [263] in the context of calculating the impact-ionization rate. More on the calculations and the applied vector discretization in this approach will be given in Section 7.3.1. Note that this further splitting of the domain implicitly leads to a higher resolution and increased accuracy, which can be seen in the comparison in Section 7.3.3.

7.2.4 Limitations of the Box Discretization Method

The box method is used in most numerical device simulation environments as it is particle conservative and has proven to deliver good results, is numerically very stable, and is relatively simple to implement. Problems arise when the Delaunay criterion is violated. This leads to obtuse elements which degenerate the accuracy due to negative flux areas $\ensuremath {A_{ij}}$ [267,264]. Also, use of the one-dimensional Scharfetter-Gummel discretization to solve multiple dimensional problems leads to the crosswind diffusion effect resulting in artificial current components perpendicular to the actual current direction [268]. The accuracy of the discretization also degrades if triangles are aligned with the hypotenuse along the current flow. As depicted in Fig. 7.7, a vanishing boundary area $A_{ki}$ leads according to (7.4) to a vanishing contribution of the current along this edge. A zig-zag characteristic of the discretized current is the result. There have been many proposals for more accurate discretizations (e.g. Patil in [267]). Some focus on the extension of the one-dimensional to a two-dimensional Scharfetter-Gummel current discretization [269,270]. But none of these extensions is as universal to use as the box integration method which is dimension independent and can be used for structured and unstructured meshes alike.

Figure 7.7: A triangle aligned with its hypotenuse along the current direction. Due to the vanishing Voronoi surface $A_{ki}$ there is no contribution to the current from the hypotenuse.