Subsections

• 4.1 Random Numbers: Direct and Indirect Method
• 4.2 Piecewise Constant Gamma Scheme
• 4.3 Self-Scattering Scheme
• 4.4 Selecting a k-vector after Scattering
  • Finding the Final Wave Vector within the Chosen Tetrahedron
  
  
• 4.5 Meshing of the Brillouin Zone
• 4.6 Calculating the DOS
  • Calculation of the area coefficients
  
  
• 4.7 Precalculated Values related to the Mesh
• 4.8 Selfconsistent Monte Carlo Scheme

4. Monte Carlo Technique

4.1 Random Numbers: Direct and Indirect Method

Monte Carlo Methods require a large amount of uniformly distributed random numbers. In computer simulations usually a sufficient long sequence of pseudo random numbers is used, which makes the results reproducible and debugging of program code generally easier. The probability density of a uniformly distributed random number is

The uniformly distributed pseudo random numbers are used to obtain random numbers according to a given distribution function . The distribution function reads [Papoulis84]

The distribution function satisfies

If the distribution function is simple enough so that its inverse can be found the direct method [Kalos86] [Jacoboni83] can be used to obtain a random number from a uniformly distributed random number by

If the probability distribution is discrete - with being the probability of the th event - the integer number can be calculated with a discrete formulation of the direct method

With this procedure the th event of possible events is chosen according to the discrete probabilities of the events. For those cases where the inverse distribution cannot be obtained the rejection method can be used.

This method consists of two steps. It involves the probability density $\tilde p(x)$ which is chosen so that it reproduces the original probability density as close as possible but with the constraint that a closed form of the inverse distribution function can be found. In a first step a random number is generated with a probability density $\tilde p(x)$ with the condition

where is a positive constant. The closed form of the inverse distribution function allows to generate using the direct method from an uniformly distributed random number . In a second step it is determined if the random number is accepted by the condition

Here, is another uniformly distributed random number. If (4.7) does not hold, is rejected and the whole procedure is repeated with two new random numbers and .

4.2 Piecewise Constant Gamma Scheme

During a Monte Carlo simulation the scattering rate depends in a complex way on the carrier state. Therefore, calculating the time of the next scattering event demands for high computing time. To speed up the calculation of the scattering rate can be rendered to a constant upper bound value, , by introducing an artificial scattering process, the so called self-scattering [Jacoboni83] process. In the following a scheme using piecewise constant -values is illustrated, which is particularly suitable for FBMC.

For the sake of a simplified notation the variables of a particle state are written as

With (3.3) and (3.4) is formulated as

$$\dot \xi = \left(n,\mathbf{v}^\mathrm{T}_n(\mathbf{r},\mathbf{k})... ...f{F}^\mathrm{T}_n(\mathbf{r},\mathbf{k})}\frac {1} {\hbar} \right)^\mathrm{T} .$$ (4.9)

Furthermore, the integral over the phase space volume and the summation over the band index is shortened as

With the above formulations and the definitions [Jungemann03]

$$\dot \xi^\mathrm{T}\nabla_{\xi} = \mathbf{F}_n^\mathrm{T}(\mathbf... ...athbf{k} + \mathbf{v}_n^\mathrm{T}(\mathbf{r},\mathbf{k})\nabla_{\mathbf{r}} ,$$ (4.11)

the BTE (3.11) for the probability density can be written as

The probability density in the phase-space is related to the distribution function by

Here, N denotes the total number of carriers in the carrier gas. With the introduction of the conditional probability density , which gives the probability that a carrier in the state at the time is found in the state at time , the formal integration of the BTE gives for the probability density

$$p(\xi,t\vert\xi_0,t_0) = p_0(\xi,t\vert\xi_0,t_0) + \int_{t_0}^t ... ...xi_1',t_1) S(\xi_1'\vert\xi_1) p(\xi_1,t_1\vert\xi_0,t_0)d\xi_1' d\xi_1 dt_1 .$$ (4.15)

In (4.15) is the conditional probability density that a particle with the initial state at the time is found in the state at time without being scattered

$$p_0(\xi,t\vert\xi_0,t_0) = \delta(\xi-\xi_\mathrm{Mot}(t\vert\xi_... ...\exp{\left(-\int_{t_0}^tS(\xi_\mathrm{Mot}(\tau\vert\xi_0,t_0))d\tau\right)} .$$ (4.16)

Here, is the solution of the equations of motion (3.3) and (3.4) for a carrier tracked from time to time with the initial position in phase-space. The evaluation of (4.16) is simplified by adding a self-scattering process to the scattering integral, which has the transition rate [Jacoboni83]

This self-scattering does not change the carrier state and therefore leaves the solution of the BTE unchanged, but renders the total scattering rate to a constant .

must be larger than the scattering rate for all possible states . The probability that a scattering event occurs at time in an interval $d \tau$ after being propagated since time without scattering is

With the direct method the time is obtained by integration of (4.19), where is a uniformly distributed random number

Solving (4.20) for yields

Evaluation of (4.21) allows a much more efficient calculation of the time of flight compared to solving the integral (4.16). As a tradeoff it introduces new artificial scattering events. Since must be larger than the physical scattering rate for all possible states , and usually varies over several orders of magnitude, a constant leads to a dominating self scattering rate for some regions. To overcome this drawback the phase-space is partitioned into small regions and for each region a local is determined

Within FBMC it is convenient to choose the region to be the volume of a mesh element of the phase space mesh. The duration of a collisionless flight is then [Jungemann03]

where is the time at which a particle passes from mesh element to . The sum on the right hand side of (4.23) is evaluated during the collisionless flight until is smaller than the time the carrier needs to reach the border of the next mesh element. The random number remains unchanged during this procedure. After the element where the collisionless flight ends is determined, the new state

is obtained by evaluating the equations of motion (3.3) and (3.4). Next one of the scattering processes (including self scattering) is selected by using the direct method

where is another uniformly distributed random number. If self scattering is selected the carrier state remains unchanged but particle statistics are updated and evaluation of equation (4.23) is started again with a new random number. If a physical scattering process is selected a position within the -space has to be determined for the scattered particle.

4.3 Self-Scattering Scheme for Ionized Impurity Scattering

Instead of evaluating the ionized impurity scattering rate with the complicated potential (3.48), a more efficient self-scattering scheme can be implemented in a Monte Carlo simulator. Within this scheme the Brooks-Herring rate multiplied by a factor is evaluated. Using (3.35) and (3.48) is chosen so that the following relation is fulfilled for any in the interval

A scattering rate of Brooks-Herring type is obtained which is always larger than the scattering rate resulting from the potential (3.48)

Here, is the scattering rate, which is evaluated with the potential (3.48). Therefore a well defined probability exists that a Coulomb scattering event is accepted [Kosina97]

Instead of solving these integrals numerically a solution by means of Monte Carlo integration can be found the following way. First a random number with the probability density is chosen. Then a second random number is chosen, which is evenly distributed between . If $p_\mathrm{r} < \vert V_\mathrm{t}(q_\mathrm{r})\vert^2$ , the scattering event is accepted, otherwise it is rejected and a self-scattering event is performed.

4.4 Selecting a k-vector after Scattering

In this section, the basic layout of a rejection algorithm to select the final state of a carrier after a scattering event is shown.

For scattering models with constant matrix elements‚ the scattering rates are proportional to the DOS at the carrier's final energy (see Section 3.3). Selection of a state after scattering then relies mainly on the calculation of the contribution to the DOS of the mesh elements including the particle's final energy. This is achieved by using tetrahedral mesh elements and linear interpolation of the energy within the mesh elements. The contribution to the DOS of the -th tetrahedron is proportional to the intersection of the equi-energy surface  [Lehmann72][Jungemann03]

The group velocity defined by equation (3.5) is constant within a tetrahedron if we consider a linear interpolation of the energy. The group velocity can be precalculated and stored in a table. The total DOS for a band is given by

$$g_n\left(\varepsilon _\mathrm{f}\right) =\frac{1}{(2\pi)^{3}}\int... ...\mathrm{f})}\frac{1}{\vert\mathbf{v}(\mathbf{k})\vert} d^{2}k=\sum_{i} g_i .$$ (4.30)

It should be noted, that with linear energy interpolation an equi-energy surface within the Brillouin zone is continuous.

The most time consuming task while computing a scattering event is the selection of a final tetrahedron containing the final energy . Several tables are calculated once at start time of the simulation to speed up this selection process.

• Table of Sorted Tetrahedrons: the tetrahedral mesh elements representing the considered part of the first Brillouin zone are sorted with respect to their lowest energy values. Additionally, the numbering of the vertices of each tetrahedron is sorted with respect to increasing energy values.
• Table of Upper Bounds for DOS: the range from the lowest to the highest occuring energy within a band is divided into small energy intervals . For each energy interval the highest DOS within all tetrahedrons is stored to build a table of upper bounds .
• Upper Bound for Energy Differences: for each energy interval the largest energy difference within the tetrahedrons which intersect a specific energy interval is stored to build a list of upper bounds .

After a scattering process is selected the final energy of the scattered particle is known and so the search for a final tetrahedron can be constricted to the tetrahedrons containing . Actually only tetrahedrons with a minimum energy in the range between and from the table of sorted tetrahedrons are considered in this search.

$$\begin{array}{l} \varepsilon _\mathrm{min} = (\varepsilon _\m... ...\varepsilon _\mathrm{max} = \varepsilon _\mathrm{f} \end{array}$$ (4.31)

Figure 4.1: Rejection technique for the selection of the final state in -space. and are uniformly distributed random numbers between zero and one.

The corresponding table indices and are calculated by a binary search. Then a tetrahedron is chosen randomly from this interval. There is a small number of tetrahedrons in the considered interval which do not contain . These are sorted out within a first rejection step.

During a second rejection step the DOS within the tetrahedron has to be evaluated using equation (4.29). Figure 4.1 shows a flow-chart of the selection procedure with its two step rejection technique. Since the whole procedure is repeated until the second rejection step is passed, has to be evaluated frequently during the simulation.

Finding the Final Wave Vector within the Chosen Tetrahedron

After selecting the final tetrahedron , a final wave vector has to be chosen within this tetrahedron. Since equi-energy areas are planes and the group velocity is constant within a tetrahedron, this problem is reduced to the determination of a uniformly distributed random point on the equi-energy plane defined by .

The equi-energy plane has the shape of either a triangle or a quadrangle as depicted in Figure 4.9 and Figure 4.10, respectively. In the case of a quadrangle, it is split into two triangles and one triangle is randomly chosen according to the surface ratio. The final -vector is then determined with

where , , and are the -vectors to the vertices of the triangle. The random numbers and are derived as [Jungemann03]

$$\begin{gather*}\begin{split}\lambda_1 &= 1 - \sqrt{1 - r_1} \lambda_2 &= r_2 (1 - \lambda_1) . \end{split}\end{gather*}$$ (4.33)

Here, and uniformly distributed random numbers.

4.5 Meshing of the Brillouin Zone

Within this work we mesh the first octant of the first Brillouin zone to represent the band structure of relaxed and biaxially strained Si (see Section 2.6). Because an octant is a larger volume than the irreducible wedge, this approach obviously increases memory consumption. However, it simplifies the manipulations needed when particles reach a boundary [Karlowatz07].

The energy bands are calculated for the irreducible wedge using the empirical pseudopotential method [Ungersboeck07a][Rieger93] and then transformed by coordinate permutation to completely fill the first octant. Two approaches of mesh generation employing structured and unstructured tetrahedral meshes are explored for several simulation setups.

The structured mesh is based on an octree approach. The domain to be meshed is devided into cubes and afterwards the cubes are divided into six tetrahedrons. To mesh the {111} surface of the Brillouin zone either one or five tetrahedrons have to be cut away from the six tetrahedrons forming one cube. The result is a structured tetrahedral mesh, whose surface is conform with the Brillouin zone boundary.

Figure 4.2(a) shows the - plane of the first octant of the Brillouin zone with the contour plot of the energy of the first conduction band. The main drawback of the octree approach is the difficulty to implement a sufficiently flexible refinement strategy to accomodate different mesh densities. In [Fischer99] an octree algorithm is proposed which can deal with different mesh densities, but the refinement zone is limited to a cubical region and is therefore not very flexible.

Figure 4.2: - plane of the first conduction band in the first Brillouin zone with (a) structured and (b) unstructured meshes. Shown is only one octant.

A more flexible way of generating unstructured meshes is to use a mesh generator which can handle arbitrary point clouds with different point densities. In this particular work DELINK [Fleischmann99] was used to generate an initial, very coarse, unstructured mesh. For different energy bands this initial mesh was refined by the so-called tetrahedral bisection method. The basic idea of this method is to insert a new vertex on a particular edge, the refinement edge, of a tetrahedron, and to cut the tetrahedron into two pieces.

In literature one can find different improvements and specifications for this algorithm (see, for example [Arnold00]). One common problem is the detection of the refinement edge. In a mesh, tetrahedrons are not isolated and inserting a new vertex influences the whole refinement edge batch of the tetrahedron if the conformity of the mesh after the refinement step should be kept, which is normally the case. To guarantee good shaped elements, a recursive refinement mechanism was chosen. This approach produces very regular, almost isotropic elements.

To guarantee a conforming mesh during the refinement procedure, all tetrahedrons sharing a common refinement edge have to be divided. A tetrahedron is said to be compatibly divisible if its refinement edge is either the refinement edge of all other tetrahedrons sharing that edge or the edge is part of the boundary of the domain. If a tetrahedron is compatibly divisible, we divide the tetrahedron and all other tetrahdedrons sharing the refinement edge simultaneously. If a tetrahedron is not compatibly divisible, we ignore it temporarily and divide a neighbor tetrahedron by the same process first. This leads to the atomic algorithm [Kossaczký94]. Figure 4.3 illustrates the recursive refinement process, where one new vertex is inserted.

Figure 4.3: The recursive refinement procedure involving four tetrahedrons with one common refinement edge.

As an input to the mesher, regions of high point densities are pre-defined by the known positions of the band-minima in the Brillouin zone of Si. The dimension of this region is chosen such that the shifted band-minima of strained Si are taken into account and so the same mesh is usable with recalculated energy values for different amounts of strain. As every band has its unique position of the band-minina, the meshes are calculated for each band seperately. For the valence bands only one one mesh is used, which is refined around the -point.

To demonstrate the benefits of mesh refinement four meshes are compared. A fine one and a coarse one of each the structured and the unstructured mesh type were generated. The number of points and tetrahedrons in the first conduction band of these meshes are shown in Table .

Figure 4.4 shows the mean electron energy as a function of the electric field for bulk Si for both structured and unstructured tetrahedral meshes. As the curves for 300K and also for 77K are grouped very close together above 10kV/cm, it can be concluded that for practical purposes the accuracy of the results in the high field regime is about the same for all meshes. Figure 4.5 shows a similar result for the velocity as a function of the electric field. These results demonstrate that the unstructured meshes perform very well in the high energy regimes, despite they contain less mesh elements than the structured meshes in that areas.

Figure 4.4: Electron mean energy versus field along [100] direction for Si at 77K and 300K.

Figure 4.5: Electron velocity versus field along [100] direction for Si at 77K and 300K.

Figure 4.6: Normalized mean energy of electrons in Si at thermal equilibrium versus temperature.

Figure 4.7: Low field mobility of electrons in Si versus temperature.

Figure 4.6 shows the normalized mean energy of electrons obtained from FBMC simulation at thermal equilibrium. The result for the unstructured meshes are in good agreement and converge for low temperatures to the theoretical equilibrium value of . While the fine structured mesh is sufficiently accurate at high temperatures, both structured meshes fail completely at low temperatures. Figure 4.7 shows the low field mobility of electrons. Again the coarse structured mesh fails, while fine structured mesh is in fair agreement with the unstructured meshes.

Table  gives an overview about computation times for different meshes. The computation times are separated into the times for mesh data structure build-up, which is required once at the beginning of the simulation, and two typical velocity calculations, one in the low field regime at kV/cm and a second one at kV/cm. For every velocity calculated the total amount of scattering events was set to . For the calculations a computer system with a Intel Pentium 4 CPU with GHz was used and the user process CPU time was measured.

One can clearly observe that the CPU time consumption is high for the structured meshes. This is mainly due the higher build-up times. With the much higher amount of mesh elements for the fine structured mesh, it takes a long time to compute the precalculated tables shown in the last section. The computation time using the unstructured fine mesh is approximately in the same range as for the coarse structured mesh, but one has to keep in mind, that the structured mesh fails completely for the average kinetic energy at temperatures less than room temperature, where the coarse unstructured mesh still gives reasonable results (see Figure 4.6).

Table 4.1: Parameter values of the meshes used for the first conduction band.

		data structure
group	granularity	points	tetrahedrons
	fine
structured	coarse
	fine
unstructured	coarse

In conclusion, the tetrahedral meshes offer a very good potential for refinement techniques. Simulation results in the high field regime show similar accuracy for properly refined meshes as for structured octree-based meshes with more than ten times the amount of tetrahedral elements. Simulations of electron mobility and mean electron energy in the low field regime show much better accuracy for refined meshes than octree-based meshes, particularly for simulations at low lattice temperatures.

Table 4.2: CPU time consumption for electron velocity simulations for a set of four different test meshes.

group	granularity	build-up time
	fine
structured	coarse
	fine
unstructured	coarse

4.6 Calculating the DOS

Figure 4.8: Tetrahedron with its defining location vectors

As shown in the last section the contribution of the tetrahedrons at a certain energy to the DOS has to be evaluated very frequently during a FBMC simulation. Therefore it is important to implement a CPU time efficient calculation algorithm. According to (4.30) the crucial step for the calculation of the DOS is to obtain the areas of the intersecting equi-energy planes of all contributing tetrahedrons .

In a first step the table of sorted tetrahedrons as defined in the last section is used to efficiently find tetrahedrons which intersect the considered energy surface. In the following both a conventional way to obtain this area and a more efficient approach based on the use of pre-calculated coefficients are presented. For the conventional calculation of the equi-energy plane the intersection points of the plane with the tetrahedron's edges are determined first.

With the vectors to of the tetrahedron's corners as depicted in Figure 4.8 the edges are derived as

$$\vect k_{02} = \vect k_2 - \vect k_0, \vect k_{03} = \vect k_3 - \vect k_0,$$ (4.34)

$$\vect k_{13} = \vect k_3 - \vect k_1, \vect k_{23} = \vect k_3 - \vect k_2.$$

The edge equation in parameter form is then

The energy is interpolated linearly

For a given energy the parameter evaluates to

The are the sought intersection points of the euqi-energy plane with the edge vectors. With (3.5) and (4.37) equation (4.35) evaluates to

$$\vect s_{ij}(\varepsilon ) = \vect k_i + \frac{\varepsilon - \var... ...ct k_i)}(\vect k_j - \vect k_i) , \quad i,j=0,1,2,3 , i \neq j .$$ (4.38)

Expression (4.38) is symmetric, . The intersection plane has either the shape of a triangle or a rectangle. If the intersection area consists of three intersection points forming a triangle , its area is calculated by

In the case that the intersection is a quadrangle, it is split into two triangles the areas of which can be calculated as in (4.39).

Calculation of the area coefficients

Figure 4.9: Triangular shaped equi-energy slice of a tetrahedron

To save computation time the calculation of intersection planes can be optimized by introducing area coefficients [Jungemann03], which are calculated in a pre-processing step and stored for each tetrahedron. The tetrahedron's vertices , and carry the energies , , , and . The energy changes linearly along the edges and is proportional to the parameter (4.37). Therefore for the case of a triangular shaped slice as shown in Figure 4.9 we obtain for the area

For energies in the range between and the intersection is given by the triangle with the area

$$\frac{1}{2} \Bigl\vert\Bigl(\vect s_{02}(\varepsilon ) - \vect s_... ...s \Bigl(\vect s_{03}(\varepsilon ) - \vect s_{01}(\varepsilon )\Bigr)\Bigr\vert$$ (4.41)

where is the first area coefficient. In the case of energies between and the intersection is also a triangle specified by . Similar to the first case we can obtain the intersection area with another area coefficient .

$$\frac{1}{2} \Bigl\vert\Bigl(\vect s_{32}(\varepsilon ) - \vect s_... ...s \Bigl(\vect s_{31}(\varepsilon ) - \vect s_{30}(\varepsilon )\Bigr)\Bigr\vert$$ (4.42)

For energies between and the intersection has the shape of a quadrangle as depicted in 4.10. Its area is obtained by subtracting the area of the inner triangle from the area of the surrounding triangle $\{ \vect s_{10}, \vect s_{13}, \vect s_{14} \}$ . The inner triangle is parameterized by the introduction of another area coefficient

$$\alpha_1 (\varepsilon - \varepsilon _0)^2 - \frac {1}{2} \Bigl\ve... ...s \Bigl(\vect s_{13}(\varepsilon ) - \vect s_{14}(\varepsilon )\Bigr)\Bigr\vert$$ (4.43)

The area coefficients , and derive as

$$\alpha_1 = \frac { A(\varepsilon _\mathrm{u})(\varepsilon _\mathr... ...mathrm{v} - \varepsilon _0)^2(\varepsilon _\mathrm{u} - \varepsilon _1)^2} ,$$ (4.44)

$$\alpha_2 = \frac { A(\varepsilon _\mathrm{u})(\varepsilon _\mathr... ...mathrm{u} - \varepsilon _0)^2(\varepsilon _\mathrm{v} - \varepsilon _0)^2} ,$$ (4.45)

and

Here, , , and .

Figure 4.10: Quadrangular shaped equi-energy slice of a tetrahedron

For the special case , an area calculation with equation (4.44) fails because both triangle areas become infinite. In this case another parameterization with the parameters and can be used, which is derived in the following. First the quadrangle shaped intersection area is split into two triangles with areas and

$$A^\mathrm{T}_{1}=\frac{1}{2}\left\vert\left(\mathbf{s}_{11}-\math... ...}_{12}\right)\times\left(\mathbf{s}_{13}-\mathbf{s}_{12}\right)\right\vert .$$ (4.47)

The terms on the right hand sides of (4.47) are given by

$$\mathbf{s}_{11}-\mathbf{s}_{10}=\frac{\left(\mathbf{k}_{2}-\mathb... ...on -\varepsilon _{0}\right)}{\varepsilon _{2}-\varepsilon _{0}}-\mathbf{k}_{0}=$$ (4.48)

$$=\left(\mathbf{k}_{0}-\mathbf{k}_{1}\right)\frac{\varepsilon -\va... ...ght)-\mathbf{k}_{\mathrm{a}}\left(\varepsilon _{2}-\varepsilon _{0}\right) ,$$

$$\mathbf{s}_{13}-\mathbf{s}_{10}=\frac{\left(\mathbf{k}_{3}-\mathb... ...on -\varepsilon _{0}\right)}{\varepsilon _{2}-\varepsilon _{0}}-\mathbf{k}_{0}=$$ (4.49)

$$=\left(\frac{\mathbf{k}_{13}-\mathbf{k}_{10}}{\varepsilon _{3}-\v... ...\right)= \mathbf{k}_{\mathrm{b}}\left(\varepsilon -\varepsilon _{0}\right) ,$$

$$\mathbf{s}_{11}-\mathbf{s}_{12}=\frac{\left(\mathbf{k}_{2}-\mathb... ...on -\varepsilon _{0}\right)}{\varepsilon _{3}-\varepsilon _{0}}-\mathbf{k}_{1}=$$ (4.50)

$$=\left(\frac{\mathbf{k}_{2}-\mathbf{k}_{1}}{\varepsilon _{2}-\var... ...}\right)=\mathbf{k}_{\mathrm{c}}\left(\varepsilon -\varepsilon _{0}\right) ,$$

and

$$\mathbf{s}_{3}-\mathbf{s}_{2}=\frac{\left(\mathbf{k}_{3}-\mathbf{... ...on -\varepsilon _{0}\right)}{\varepsilon _{3}-\varepsilon _{0}}-\mathbf{k}_{1}=$$ (4.51)

$$=\left(\mathbf{k}_{1}-\mathbf{k}_{0}\right)\frac{\varepsilon -\va... ...ght)-\mathbf{k}_{\mathrm{d}}\left(\varepsilon _{3}-\varepsilon _{0}\right) .$$

Now, the area of the quadrangle can be written in a parameterized formulation

where and are pre-calculated parameters for the given tetrahedron:

$$\alpha_4=\left\vert\mathbf{k}_{\mathrm{a}}\times\mathbf{k}_{\math... ...\vert+\left\vert\mathbf{k}_{\mathrm{c}}\times\mathbf{k}_{\mathrm{d}}\right\vert$$ (4.53)

$$\alpha_5=\left(\varepsilon _{2}-\varepsilon _{0}\right)\left\vert... ...)\left\vert\mathbf{k}_{\mathrm{c}}\times\mathbf{k}_{\mathrm{d}}\right\vert .$$ (4.54)

Since the regular tetrahedron and the tetrahedron where equals are completely alternative appearances, the area coefficient is not needed in the latter case. Therefore the area calculation is in any case completely parameterized with four pre-calculated area coefficients per tetrahedron.

4.7 Precalculated Values related to the Mesh
Structure

As shown in Section 4.5 it is convenient to use tetrahedral meshes for discretization of the Brillouin zone. It has also been shown that frequently used values like the area coefficients can be precalculated to improve performance during simulation. These values are stored related to the tetrahedral mesh elements, where they are either connected to the vertices, the surfaces or the volume.

The numbering of the vertices of each tetrahedron is sorted after increasing energies as explained in Section 4.4. The energies und positions of the vertices are stored in a table.

There are also several precalculated values with are connected to the tetrahedron volumes: the energy gradient, which is used to calculate the group velocity, the area-coefficients to (see Section 4.6), and the upper bounds for the scattering rates for different types of scattering mechanisms.

For each of the four surfaces of a tetrahedron the type of the surface and a pointer to the neighbor tetrahedron are stored. Figure 4.11 shows the definition of the surface types. This value is used for efficient detection of the surface that would be crossed when a carrier has reached a border of the meshed part of the Brillouin zone. If the surface of the tetrahedron is inside of the meshed domain then the value of the surface type is defined as zero.

Figure 4.11: Surface types for one octant of the Brillouin zone in -space [Wagner04].

Figure 4.12: Symmetry properties of the surface 7. Solid lines are boundaries of irreducible wedges, with one wedge indicated in blue. Dashed-dotted lines indicate the additional symmetry needed for the mapping procedure.

The pointer to the neighbors fulfills a special purpose if a tetrahedron's surface is part of the boundary of the meshed domain. In this case the pointer is set to the tetrahedron to which the carrier's position is mapped back when leaving the domain. This mapping procedure requires tetrahedrons with a surface mesh structure where the initial surface elements are congruent with the final elements. Because the mesh is built by the permutation of irreducible wedges as explained in section Section 2.6.5, this feature is guaranteed for the surfaces one to six. For surface seven the initial mesh must provide an additional symmetry [Jungemann03] as shown in Figure 4.12.

4.8 Selfconsistent Monte Carlo Scheme

The Boltzmann equation can be solved either with a given electric field or with a selfconsistent field using an iteration scheme as shown in Figure 4.13. After each time step the carrier concentration is evaluated and the Poisson equation is solved

for the electrostatic potential which in turn gives a new field for the next iteration step. In (4.55) denotes the dielectric constant and is the space charge density

Here, denotes the electron concentration, is the hole concentration, and are the concentrations of the ionized donors and acceptors, respectively.

Figure 4.13: Flowchart of a selfconsistent Monte Carlo device simulation for a given total simulation time .