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.
The box discretization method reformulates the equations based on this Voronoi tessellation and can be equally used in two and threedimensional simulation domains. Since the twodimensional 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
(7.2) 
(7.5) 

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 driftdiffusion 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 introduced in (7.1) are the dielectric flux density for Poisson's equation and the current densities for the continuity equation, respectively. The generation terms in (7.1) are the charge density and the carrier generation/recombination rates / 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.e. the vertices, and the box volumes ( ). Second, a list of all edges connecting the vertices, together with the edge length ( ), i.e. the distance between the vertices, and the common surface between the neighboring boxes ( ). 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 setup of the simulator can be realized.
The calculation of the onedimensional projection of the flux on the edge depends on the type of the differential equation. For the dielectric flux density in Poisson's equation the flux from box to box along the connecting edge through the common boundary area is This density is commonly approximated using the directional derivative of the electrostatic potential:
For the current densities of the carrier type the discretization of is not as straightforward. Insertion of the driftdiffusion 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 ScharfetterGummel method [39] instead. Here, the driftdiffusion current equations (4.7) and (4.8) are used to solve the onedimensional carrier concentrations, and respectively, along the edge. The boundary conditions of the carrier concentrations are given using the values at the corresponding vertices and The values of and and and are considered constant along the edge. Solving this onedimensional differential equation results in
(7.9) 
Up to now, the generation term on the right hand side of (7.4), which is the charge density 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

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. MINIMOSNT [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 Using the naming conventions from Fig. 7.6(b), this reads
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:
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 [267,264]. Also, use of the onedimensional ScharfetterGummel 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 leads according to (7.4) to a vanishing contribution of the current along this edge. A zigzag 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 onedimensional to a twodimensional ScharfetterGummel 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.
