Subsections

• 7.3.1 Discretization Approach by Laux
• 7.3.2 Box Discretized Vector Quantities
  • 7.3.2.1 Scheme A
  • 7.3.2.2 Scheme B
  • 7.3.2.3 Differentiability of the Discretized Vector Quantities
  
  
• 7.3.3 Comparison of the Discretization Schemes

7.3 Vectors in Discretized Systems

There are different methods of how directed quantities can be represented in box discretized systems. In the initial divergence equation (7.1) the vector quantities $\ensuremath{\mathbf{F}}$ is transformed into one-dimensional edge discretized values $F_{ij}$ which are subsequently used to assemble the box and flux equations. The single fluxes are estimated using quantities defined at the end nodes of the edges, like in relation (7.6) for the dielectric flux and (7.7) for the current density. No vector quantities are required to assemble the basic equations. However, vector attributes like the electric field and the current density are needed for the calculation of physical models. The electric field and the current density are required, for example, to model the high-field mobility or the impact-ionization rate in the drift-diffusion framework. For the representation of attributes as vector quantity, additional methods are required, some are represented in this section.

7.3.1 Discretization Approach by Laux

In the work of Laux et al. [263] a discretization scheme is introduced for the accurate evaluation of the impact-ionization rate. The model assumes an unstructured triangular mesh. Since the potential and the electric field are based on linear functions, an electric field vector can be evaluated in a straightforward manner, consistent with the Poisson discretization. In a triangle $ijk,$ each pairwise linear combination of $\ensuremath{\bm{\xi}}_{ij}E_{ij},$ $\ensuremath{\bm{\xi}}_{jk}E_{jk},$ and $\ensuremath{\bm{\xi}}_{ki}E_{ki}$ gives the same constant value for the electric field $\ensuremath{\mathbf{E}}_{ijk}.$

Due to the non-linear dependence of the current density $\ensuremath{\mathbf{J}}$ on the carrier concentrations in the Scharfetter-Gummel discretization, each linear combination of $\ensuremath{\bm{\xi}}_{ij}J_{ij},$ $\ensuremath{\bm{\xi}}_{jk}J_{jk},$ and $\ensuremath{\bm{\xi}}_{ki}J_{ki}$ gives a different current density vector $\ensuremath{\mathbf{J}}.$ The system in the triangle is overdetermined. Laux suggested to partition each triangle into three avalanche regions associated with each edge as shown in Fig. 7.8.

Figure 7.8: A triangular element showing the avalanche region defined by Laux et al. [263] (black hatched area) and the weighting surfaces $A_{ij}'.$ The part contributing to the right hand side of box $i,$ $V_{e,ij},$ is bordered by the red solid line.

In this work, a current vector is defined for each avalanche region, i.e. for each edge. Note, that in the notations used here, the indices for the carrier type specification have been omitted. The current $\ensuremath{\mathbf{J}}_{ij}$ associated to the edge between the vertices $i$ and $j,$ for example, is defined as

$$\ensuremath{\mathbf{J}}_{ij} = \frac{\ensuremath{\mathbf{J}}_j^* A_{jk}' + \ensuremath{\mathbf{J}}_i^* A_{ik}'}{A_{jk}' + A_{ik}'} ,$$ (7.13)

weighted with the surface elements $A_{ij}'.$ Note that $A_{ij}$ is the whole surface between two boxes $i$ and $j,$ whereas $A_{ij}'$ is only the part within the current element. $\ensuremath{\mathbf{J}}_i^*$ is the current defined by the linear combination of the two to the vertex adjacent edges, given as

$$\begin{pmatrix}\ensuremath{\bm{\xi}}_{ij}^T \ensuremath{\bm{\x... ...} \ensuremath{\mathbf{J}}_i^* = \begin{pmatrix}J_{ij} J_{ik} \end{pmatrix} .$$ (7.14)

The electric field $\ensuremath{\mathbf{E}}$ and the current density $\ensuremath{\mathbf{J}}_{ij}$ can now be used to calculate a field dependent model, in this context the impact-ionization rate, for the avalanche region. This region contributes as $g_{e,ij}$ to the two boxes $i$ and $j$ which have to be assembled according to (7.12).

This approach for vector discretization gives good results as discussed in Section 7.3.3. The limitation to triangular grids (although an extension to tetrahedra is possible) and the higher complexity during equation assembly seem to be the main disadvantages of this method compared the box based vector discretization schemes presented in the next section.

7.3.2 Box Discretized Vector Quantities

Considering the calculation of the summand in the box approach (7.10), it seems to be most convenient to estimate the vectorial attribute for each box volume, i.e. for each vertex. The model evaluation within the box can then be performed straightforwardly, since all quantities, scalar and vectorial, are then available for the whole box. The results of model evaluation can be directly applied to the box integration equation. In the following, two schemes of vector discretization within boxes are presented. In addition to the simple coupling to the box discretization method, both approaches give accurate approximations for homogenous fields and are numerically stable.

7.3.2.1 Scheme A

The first scheme follows the derivation of the box discretization scheme itself [271]. Similar to the discussion regarding the right hand side of the box integral, the discretized vector quantities are assumed constant over the whole Voronoi box volume. The derivation of an according discretization scheme is shown for the electric field and can be generalized to auxiliary gradient fields. The electric field is defined as

$$\ensuremath{\mathbf{E}}= - \nabla \psi = - \nabla (\psi - \psi_\mathrm{ref}) ,$$ (7.15)

where $\psi_\mathrm{ref}$ is an arbitrary reference potential. Similar to the method used to derive the box discretization method, an integration over a volume $V$ is performed which leads to

$$\int_V \ensuremath{\mathbf{E}}\ensuremath{ \mathrm{d}}V = - \int_V \nabla (\psi - \psi_\mathrm{ref}) \ensuremath{ \mathrm{d}}V .$$ (7.16)

Applying Gauss's theorem the equation can be written as

$$\int_V \ensuremath{\mathbf{E}}\ensuremath{ \mathrm{d}}V = - \oin... ...V} (\psi - \psi_\mathrm{ref}) \ensuremath{ \mathrm{d}}\ensuremath{\mathbf{A}},$$ (7.17)

where $\ensuremath{\mathbf{A}}$ is the outwardly oriented surface vector of $\partial V.$ Originally, Gauss's theorem is defined for divergence operation on vector fields, but it can be shown, that it is also applicable for gradients of potential fields. For the box $i$ with $N_i$ neighboring boxes $j,$ the equation can be estimated to

$$V_i \ensuremath{\mathbf{E}}_i = - \sum_{j \in N_i} \ensuremath{A_{ij}}\bm{\xi}_{ij} (\psi_{ij} - \psi_i) .$$ (7.18)

In this derivation the electric field is constant over the box volume and the integral was split into the single box boundary surfaces. $\psi_i$ was used as reference potential $\psi_\mathrm{ref}.$ The potential $\psi_{ij}$ is the potential at the boundary between the boxes, which is the mid-point of the edge and can be estimated as

$$\psi_{ij} = \frac{\psi_j + \psi_i}{2} .$$ (7.19)

Introducing $E_{ij}$ as the projection on the electric field by following (7.6) one gets

$$E_{ij} = \bm{\xi}_{ij} \cdot \ensuremath{\mathbf{E}}\approx -\frac{\psi_j - \psi_i}{d_{ij}} .$$ (7.20)

Combining equations (7.18) to (7.20), the final relation for the electric field can be given as

$$\ensuremath{\mathbf{E}}_i = \frac{1}{2 V_i} \sum_j \ensuremath{A_{ij}}\ensuremath{d_{ij}}\bm{\xi}_{ij} E_{ij} .$$ (7.21)

This scheme can be applied in Voronoi grids. The unstructured neighborhood information is sufficient to calculate the vectors. For the special case of homogenous fields $\ensuremath{\mathbf{E}}= \ensuremath{\mathbf{E}}_H$ with the electrostatic potential $\psi(\ensuremath{\mathbf{x}})=-\ensuremath{\mathbf{E}}_H \cdot \ensuremath{\mathbf{x}},$ the discretization step in (7.18), the estimation of the mid potential in (7.19), and the projected electric field in (7.20) are exact. This can be verified for orthogonal grids by entering $E_{ij} = \ensuremath{\bm{\xi}}_{ij}^T \ensuremath{\mathbf{E}}_H$ into (7.21) and using the relation $\ensuremath{\bm{\xi}}_{ij} \left( \ensuremath{\bm{\xi}}_{ij}^T \ensuremath{\ma... ...xi}}\ensuremath{\otimes}\ensuremath{\bm{\xi}} \right) \ensuremath{\mathbf{E}}_H$ which leads to

$$\ensuremath{\mathbf{E}}_i = \frac{1}{2V_i} \left[ \sum_j \ensurem... ...}}_{ij} \right) \right] \ensuremath{\mathbf{E}}_H = \ensuremath{\mathbf{E}}_H .$$ (7.22)

The discretization scheme is also analyzed in a one-dimensional formulation applying a linearly changing electric field $E(x) = -2 \alpha x$ and an corresponding quadratic electrostatic potential $\psi(x) = \alpha x^2.$ Using only the x-axis and the naming convention from Fig. 7.9, (7.21) can be reduced for the linearly changing field to

$$E(x_0) =-\alpha (x_1 + x_2) .$$ (7.23)

This fits exactly in the equidistant case, where $x_0 = \left.\left(x_1 + x_2\right)\middle/2\right.,$ but for non-equidistant cases, the result is not exact.

7.3.2.2 Scheme B

The second discretization scheme is an extension of the finite difference method and is based on a scheme suggested by Fischer [272]. The evaluation concentrates on the vector at the vertex itself. Considering the box 0 in a non-equidistant orthogonal mesh depicted in Fig. 7.9 and its neighboring box 1 (not shown explicitly), the electric field along the edge $d_{01}$ can be expressed as $E_{01}=-\mathrm{d}\psi(x)/\mathrm{d}x.$

Figure 7.9: Rectangular Voronoi box of mesh point 0 with its neighboring points $1$ -$4$ in a non-equidistant orthogonal mesh. The contributing field components $E_{ij}$ from the edges are depicted. (Idea of this picture taken from [272])

At the boundary between the two boxes, i.e. the midpoint between 0 and 1, the finite difference method gives

$$E_{01} = - \frac{\psi_1 - \psi_0}{d_{01}} ,$$ (7.24)

where $\psi_i=\psi(x_i).$ Note that this is the same result as in (7.20). The electric field $E^x_0$ in direction $\ensuremath{\bm{\xi}}_x$ at mesh point 0 is expressed as a linear interpolation between $E_{01}$ and $E_{02},$ the same is done for the component in direction $\ensuremath{\bm{\xi}}_y$ for $E^y_0,$ which leads to

$$E^x_0 = \frac{ \frac{E_{01}}{x_1-x_0}+\frac{E_{02}}{x_2-x_0}}{ \frac{1}{x_1-x_0}+\frac{1}{x_2-x_0}} ,$$ (7.25)

$$E^y_0 = \frac{ \frac{E_{03}}{y_3-y_0}+\frac{E_{04}}{y_4-y_0}}{ \frac{1}{y_3-y_0}+\frac{1}{y_4-y_0}} .$$ (7.26)

Following [272], equations (7.25) and (7.26) are extended to cover edges not aligned to the coordinate axis, reading

$$\frac{1}{x_j-x_i} \Rightarrow \frac{x_j-x_i}{(x_j-x_i)^2+(y_j-y_i)^2} = \frac{x_j-x_i}{\ensuremath{d_{ij}}^2} .$$ (7.27)

The distance $\ensuremath{d_{ij}}$ between the points $i$ and $j$ is defined as $\ensuremath{d_{ij}}=\vert\ensuremath{\ensuremath{\mathbf{d}}_{ij}}\vert$ using

$$\ensuremath{\ensuremath{\mathbf{d}}_{ij}}= \ensuremath{\mathbf{x}... ...{x}}_i = \ensuremath{\begin{pmatrix}{x_j - x_i} {y_j - y_i} \end{pmatrix}} .$$ (7.28)

By using this generalized formulation (7.27) and assuming point 0 as a generalized point $i$ of a box $i$ with its neighboring boxes $j,$ equations (7.25) and (7.26) can be reformulated to

$$\left( \sum_j \frac{x_j-x_i}{\ensuremath{d_{ij}}^2} \right) E_i^x = \sum_j \frac{x_j-x_i}{\ensuremath{d_{ij}}^2} E_{ij} ,$$ (7.29)

$$\left( \sum_j \frac{y_j-y_i}{\ensuremath{d_{ij}}^2} \right) E_i^y = \sum_j \frac{y_j-y_i}{\ensuremath{d_{ij}}^2} E_{ij} .$$ (7.30)

Introducing the electric field vector $\ensuremath{\mathbf{E}}_i = ( E^x_i \; E^y_i )^{\textrm{T}},$ using $\ensuremath{\bm{\xi}}_x^T \ensuremath{\mathbf{E}}_i = E^x_i,$ and adding the flux area $\ensuremath {A_{ij}}$ as additional weight factor on both sides of the equations, the extended discretization scheme becomes

$$\sum_j \frac{x_j-x_i}{\ensuremath{d_{ij}}^2} \ensuremath{A_{ij}}\... ... \ensuremath{\ensuremath{\mathbf{d}}_{ij}}^\mathrm{T} \ensuremath{\mathbf{E}}_i = \sum_j \frac{x_j-x_i}{\ensuremath{d_{ij}}^2} \ensuremath{A_{ij}}E_{ij} ,$$ (7.31)

$$\sum_j \frac{y_j-y_i}{\ensuremath{d_{ij}}^2} \ensuremath{A_{ij}}\... ... \ensuremath{\ensuremath{\mathbf{d}}_{ij}}^\mathrm{T} \ensuremath{\mathbf{E}}_i = \sum_j \frac{y_j-y_i}{\ensuremath{d_{ij}}^2} \ensuremath{A_{ij}}E_{ij} ,$$ (7.32)

which can be combined to the single summation

$$\sum_j \frac{\ensuremath{\ensuremath{\mathbf{d}}_{ij}}}{\ensurema... ...emath{\mathbf{d}}_{ij}}}{\ensuremath{d_{ij}}^2} \ensuremath{A_{ij}}E_{ij} .$$ (7.33)

Using $\ensuremath{\ensuremath{\bm{\xi}}_{ij}}= \ensuremath{\ensuremath{\mathbf{d}}_{ij}}/ \ensuremath{d_{ij}}$ and introducing the tensor product $\ensuremath{\ensuremath{\bm{\xi}}_{ij}}\ensuremath{\otimes} \ensuremath{\ensur... ...h{\ensuremath{\bm{\xi}}_{ij}}\ensuremath{\ensuremath{\bm{\xi}}_{ij}}^\mathrm{T}$ in the equation, the full formulation of the discretization scheme reads

$$\underbrace{\sum_j \frac{\ensuremath{A_{ij}}}{d_{ij}} \left( \ens... ...= \sum_j \frac{\ensuremath{A_{ij}}}{d_{ij}} \ensuremath{\bm{\xi}}_{ij} E_{ij} .$$ (7.34)

Note that (7.25) and (7.26) are still retained and can be extracted by using $\ensuremath{\bm{\xi}}_{ij}=( 1 \; 0 )^{\textrm{T}}$ and $\ensuremath{\bm{\xi}}_{ij}=( 0 \; 1 )^{\textrm{T}},$ respectively. $\ensuremath{\mathbf{E}}_i$ at the left side of (7.34) can be taken out of the sum and the remaining part of the sum results in a pure geometry dependent matrix, which is calculated once in the beginning of the simulation. This allows the convenient formulation of the final discretization rule for a vector $\ensuremath{\mathbf{E}}_i$ in point $i$ as

$$\ensuremath{\mathbf{E}}_i = \ensuremath{\bm{M} _i^{-1}}\sum_j g_{ij} \ensuremath{\bm{\xi}}_{ij} E_{ij} ,$$ (7.35)

$$\ensuremath{\bm{M} _i}= \sum_j g_{ij} \ensuremath{\bm{\xi}}_{ij} \ensuremath{\otimes}\ensuremath{\bm{\xi}}_{ij} ,$$ (7.36)

$$g_{ij} = \frac{\ensuremath{A_{ij}}}{\ensuremath{d_{ij}}} .$$ (7.37)

$\ensuremath{\bm{M} _i}$ is called the geometry matrix and $g_{ij}$ the geometry factor, both of which depend only on the geometry and the mesh and do not change during the simulation.

Similar to Scheme A, the validation for homogenous fields using $\ensuremath{\mathbf{E}}= \ensuremath{\mathbf{E}}_H$ and the electrostatic potential $\psi(\ensuremath{\mathbf{x}})=-\ensuremath{\mathbf{E}}_H \cdot \ensuremath{\mathbf{x}}$ is shown. Again, the relations $E_{ij} = \ensuremath{\bm{\xi}}_{ij}^T \ensuremath{\mathbf{E}}_H$ and $\ensuremath{\bm{\xi}}_{ij} \left( \ensuremath{\bm{\xi}}_{ij}^T \ensuremath{\ma... ...xi}}\ensuremath{\otimes}\ensuremath{\bm{\xi}} \right) \ensuremath{\mathbf{E}}_H$ are inserted into (7.35) which leads to

$$\ensuremath{\mathbf{E}}_i = \ensuremath{\bm{M} _i^{-1}}\sum_j g_{... ...}}\ensuremath{\bm{M} _i}\ensuremath{\mathbf{E}}_H = \ensuremath{\mathbf{E}}_H .$$ (7.38)

This box based discretization scheme therefore is always exact for homogenous fields. Also the one-dimensional verification of linearly changing fields $E(x) = -2 \alpha x$ and a corresponding quadratic electrostatic potential $\psi(x) = \alpha x^2$ has been done. Inserting these values in (7.35) and solving in one dimension using the labels as in Fig. 7.9 leads to the exact result

$$E(x_0) =-2 \alpha x_0 .$$ (7.39)

The geometry matrix $\ensuremath{\bm{M} _i}$ introduced in this discretization scheme needs to be inverted for evaluation. The matrix results from a sum of symmetric matrices $\ensuremath{\bm{\xi}}_{ij} \ensuremath{\otimes}\ensuremath{\bm{\xi}}_{ij}$ whose determinants equal to 0 and whose main diagonals are positive. The sum of matrices with these constraints and with non-negative determinants also result in a symmetric matrix with positive main diagonal and a non-negative determinant. If at least two of the participating matrices are linearly independent, the determinant of the geometry matrix is positive which is the case as long as the Delaunay criterion is fulfilled. The inverse of the geometry matrix can therefore always be calculated.

7.3.2.3 Differentiability of the Discretized Vector Quantities

For the numerical solution process using Newton's method it must be possible to calculate all partial derivatives on quantities in all neighboring boxes. The derivatives need to be added to the Jacobian matrix. This makes it necessary that the discretization schemes are also differentiable on quantities $\eta_k$ associated to a mesh point $k.$ This is indeed possible and one obtains for Scheme A

$$\ensuremath{\frac{\partial \ensuremath{\mathbf{E}}_i}{\partial \e... ...nsuremath{\mathbf{d}}_{ij} \ensuremath{\frac{\partial E_{ij}}{\partial \eta_k}}$$ (7.40)

and for Scheme B

$$\ensuremath{\frac{\partial \ensuremath{\mathbf{E}}_i}{\partial \e... ...nsuremath{\bm{\xi}}_{ij} \ensuremath{\frac{\partial E_{ij}}{\partial \eta_k}} .$$ (7.41)

In both discretization schemes, the existence of $\ensuremath{\frac{\partial E_{ij}}{\partial \eta_k}}$ is sufficient to calculate $\ensuremath{\frac{\partial \ensuremath{\mathbf{E}}_i}{\partial \eta_k}}.$

7.3.3 Comparison of the Discretization Schemes

In this comparison, the vector discretization scheme proposed by Laux and the two box discretized schemes where implemented in MINIMOS-NT and used to calculate the impact-ionization rate. The generation integral on the right hand side of the continuity equation is assembled using (7.10) and (7.12), respectively. The implementation of the scheme by Laux is only limited to two-dimensional domains using triangular meshes, whereas the two other schemes are dimension and mesh independent and the same program code can be used for two- and three-dimensional simulations.

Simulation results from two devices are presented. The first device is a diode which was selected to investigate effects in a simple one-dimensional device. The second device used for the comparison is a parasitic n$^+$ -p-n-n$^+$ structure of a smart power device which is significantly influenced by the two-dimensional extension. The diode uses an equidistant mesh and is investigated in reverse-biased operating condition for different mesh spacings. The parasitic bipolar smart power structure is simulated in snap-back, a state the device can be driven in during voltage peaks on the power line.

Simulations on the diode structure clearly show that the mesh dependency is higher for the two box based discretization schemes. Using a high mesh density results, as expected, in a comparable output for all three discretization methods. Using a coarser mesh spacing, the results using Laux's scheme change very little, whereas the two box based schemes show larger deviations. In Fig. 7.10 an example using the diode clearly shows that an increased mesh spacing leads to a shift of the breakdown voltage using the box based schemes, whereas only a very small shift is observed using the scheme by Laux.

Figure 7.10: Comparison of the simulation results (a) and the error (b) for a reverse-biased diode near breakdown. Only Scheme B has been used in the comparison.

(a) Diode Current

(b) Snap-Back

The results for the snap-back simulation in the smart power device are depicted in Fig. 7.11 and show only small deviations between the three schemes. The two box based schemes again show a voltage shift in comparison to the scheme by Laux. The latter one fits well to the reference solution which was generated using a high mesh density (not shown in the figure). The reason for the stronger mesh dependence of the box based schemes can be found in the implicitly finer discretization used in Laux's scheme.

Figure 7.11: Snap-back simulations in a smart power structure using different vector discretization schemes (a). The influence of the discretization scheme on the Newton iteration process (b) during the transition marked in (a) is very small.

(a) Snap-Back

(b) Convergence Behavior

The influence of the vector discretization scheme on the convergence behavior is also investigated on the reduced smart power structure as shown in Fig. 7.11(b). Only very little differences were noted and no trend favoring one or another scheme was observed. The convergence process is tracked using the norm of the right hand side. The investigated simulation step is a numerically critical current level step at the triggering phase of the snap-back. It can be clearly seen that the choice of discretization method has only very little influence, despite the critical simulation step.

Although the Laux scheme gives results with a higher accuracy, the influence on the convergence behavior seems to be negligible. On the other hand, the box based schemes can be coupled straightforwardly to the box discretization scheme. Additionally, it is possible to reuse the same implemented code for arbitrary meshes in all dimensions. Throughout this work, the box based Scheme B has been used for all drift-diffusion simulations.