The Monte Carlo method is a stochastic approach to integrate functions in general and to solve integral equations in particular [21, 34, 43, 51, 59, 61].
The Monte Carlo method for carrier transport in semiconductors calculates trajectories of carriers in - and -space under the influence of acceleration by the electric field and of scattering mechanisms. The duration of free-flight, the type of the scattering mechanism, and the final state after scattering are calculated using random numbers. With a sufficiently large number of trajectories, the averages of the attributes can be calculated [18, 60].
One drawback of stochastic simulations is that the statistical error of the estimator is declining with the factor , where is the number of random events. In other words, if one wants to reduce the error by a factor , the calculation cost will increase by .
An advantage of the Monte Carlo method is that individual trajectories of carriers are much simpler to calculate than solving the BTE with a deterministic approach. Furthermore, the Monte Carlo method allows one to use full-band structures to estimate the carrier distribution function . Because of these advantages the Monte Carlo method is used in many cases as a reference method for simpler transport models [60].
As mentioned above, the Monte Carlo method is in general a stochastic method to solve integrals. A integral over a function ,
can be expressed as an expectation value of some random variable. For this purpose, the function is factorized:
where is a density function of the Monte Carlo samples satisfying
The integral in (3.34) represents an expectation value . The Monte Carlo method estimates the expectation value by a sample mean:
where N is the number of sampling points [61].
Monte Carlo methods rely on uniformly distributed random numbers. Since the generation of random numbers on computers is difficult, pseudo-random numbers are used [9, 47, 75]. These pseudo-random numbers have one big advantage for testing: if the same seed is used every time, the sequence is reproducible [P3, 78].
The uniformly distributed pseudo-random number , used in this method, has the property [51]
For the generation of a random number with a given probability density , its cumulative distribution function is needed [9, 46, 51, 52]:
with
A distributed random number can be calculated through the inverse of the cumulative distribution function
where is a pseudo-random number of type (3.37). However, this method is only applicable if the analytic integral of the function is possible to evaluate. Otherwise, different approaches to generate distributed pseudo-random numbers are needed.
For the generation of a random number with a given probability density , where is bounded in the finite interval , a uniform distributed random number with the properties
is chosen for the sample point . If the random number lies in the range of
it gets accepted. Else the random number gets rejected [12, 46]. The accepted random number is distributed. This method is always applicable, with any bounded function in a finite interval. However, if the function is heavily peaked, many rejections will occur and the computational expenses will become high.
To keep the rejection rate small, the function which represents the maximum of the probability density, , could be replaced by an analytically integrate-able function . In the interval , must be always greater than .
Here, a random number can be directly obtained from with the inversion method. This random number can than be applied to the rejection method, where the number of rejections can be significantly reduced by finding a suitable function . In the case, that is chosen as a constant, this method will result in the rejection method, described above.
The solution of the BTE for carrier transport can be estimated with the Monte Carlo method. One random parameter in this method is the time of the free-flight of a carrier. This parameter depends on the total scattering rate . The probability for the scattering of a carrier in an interval at a time is , where represents the and dependent scattering rate at the time :
With the assumption, that the carrier scatters at , the cumulative distribution function can be calculated. In this case, no other scattering process takes place until the time . Thus, the carrier is in free-flight for the duration of [60]:
Using the inversion method, the duration of the free-flight can be calculated from the equation as
where represents a uniform distributed random number [46]. The integration takes place along the trajectory of the carrier , which can be acquired by integrating the equations of motion (2.22) and (2.23).
The full-band simulations use constant total scattering rates , as shown in Fig. 3.1 [110]. In case is constant, the duration of the free-flight is [P3]:
Because must be larger than the sum of the scattering rates of all physical scattering processes , the self-scattering is more likely to occur in areas where is significant smaller than . For this reason can be approximated by local, picewise constant values . Since these local values are smaller than the global , self-scattering is reduced. The duration of a collision-less free-flight is given by [51]
where is the time at which the particle passes a change in during the free-flight .
After the free-flight, another uniformly distributed random number,
is used to select the scattering process, as described in [51, 60]. The actual scattering rates from all the implemented processes must be smaller than or equal to the constant total scattering rate ,
where is the rate of the scattering process. The scattering process is chosen when
If the random number is in the range of
self scattering occurs [60].
Equilibrium distributions at given temperatures are needed for injecting particles in a Monte Carlo simulation. This section shows three different methods how an equilibrium distribution can be sampled.
A single-particle Monte Carlo simulation is performed in a box with constant length. To obtain an equilibrium distribution, the electric field is set to zero. Every time the carrier hits a boundary and gets reflected, its state is added to the sample, as illustrated in Fig. 3.2. The distribution of the generated sample represents a velocity weighted Maxwellian distribution [65]
where is the velocity component perpendicular to the boundary.
In the statistical average of an attribute the weighting factor is obtained by replacing by [P3]:
where represents the number of sampled states.
A single-particle Monte Carlo simulation is performed for a uniform semiconductor at zero field. The equilibrium distribution is obtained by sampling the particle state before a scattering event occurs (before scattering method). With this method the distribution of the obtained -values is a Maxwellian distribution weighted with the scattering rate [46, 60].
When the mean value of an attribute is calculated, the weighting factor has to be taken into account:
Again, a single particle Monte Carlo simulation is performed for a uniform semiconductor at zero field. Now, the trajectory is sampled at constant time intervals. In this case the equilibrium distribution is directly obtained [46] and the average can be calculated from the sampled states without any additional weighting factors.
With the numerical representation of the dispersion relation in the first Brillouin zone, averages such as mean kinetic energy or injection velocity at thermal equilibrium can be calculated by numerical integration over the equilibrium distribution [P3].
This method is faster than Monte Carlo integration. In contrast to the Monte Carlo methods, this method can be only applied in thermal equilibrium, where the distribution function is known to be a Maxwell-Boltzmann or a Fermi-Dirac distribution. Here we assume a Maxwell-Boltzmann (MB) distribution. The average of an attribute is calculated like [74]
where BZ denotes the whole Brillouin zone, is the band index, with the Boltzmann constant and the Temperature . In statistical mechanics a parameter called partition function is introduced [74] .
With the partition function (3.58) the statistical average (3.57) becomes:
For a detailed description of the numerical integration see Appendix A.1 or [P3].
With the Monte Carlo method, average values of attributes of interest can be estimated, such as electron density, carrier velocity, energy and many more. For the estimation of local attributes using the single particle method, there is commonly one method used: the before scattering method. This work presents a second method to estimate local attributes based on the box-sampling method. Further, the estimation of global attributes, for example, the current through a contact is shown.
The before scattering method is obtaining the values of the attributes of interest before a scattering event takes place. These values are weighted with the total scattering rate , as shown in (3.55) [51, 60]:
In a spatial discretization every grid point in a device has a volume assigned. Inside this Volume , the closest grid point is , where is an index for every grid point in space. This kind of discretization volume is also known as Voronoi volume. The averages of the local attributes can be built for every discretization volume separately in the manner of [88]:
Here, the summation runs over the before scattering states , where the scattering takes place inside the Volume [60].
One drawback of this method is, that if no scattering event happens while a particle traverses a discretization volume, no contributions are made to the sum in (3.60).
The boundary method gathers statistical information when a particle crosses a boundary from one Voronoi volume to another. Therefore, the -values are weighted with the velocity component of the particle perpendicular to the boundary [65]. In the manner of the box method, the mean values of the local attributes can be calculated like:
where is the point on the -th edge, where one particle is crossing the -th time the edge between two Voronoi volumes.
The benefit of this method is, that also in regimes where there is barely no scattering, sampling values can be obtained. Therefore, statistical averages of local attributes even in areas with ballistic transport can be estimated.
The before-scattering method and the boundary method can be combined. This combination of methods gives a smaller statistical error than the ones of each method individually.
To investigate this combined averaging method, a silicon diode with abrupt junctions was chosen. The doping levels are cm and cm, respectively. Fig. 3.6 shows the conduction band edge and the electron density for an applied voltage of .
Fig. 3.7 compares the different averaging methods in the diode in the area of the junction at 200 nm. The simulation with constant grid size shows a nearly constant factor between before-scattering and boundary method. The diode with variable grid size shows a clear dependence of the number of entries of the before-scattering method on the grid size. At the grid size is becoming so big, that the before-scattering method has an advantage over the boundary method. Nevertheless, the combination of the two methods always gathers more entries than one method alone.
The Monte Carlo method is capable of estimating global attributes as well as local ones. The current is chosen here exemplarily for a global attribute to estimate. The following calculations of the current in a two dimensional device can be found in [P3]. The definition of the current density reads:
where is the unknown solution of the BTE and is the particle charge. The current through a contact in a two-dimensional device can be calculated as
where represents the integration over the contact area and is the width of the contact, as shown in Figure 3.8. Introducing the velocity weighted distribution of the contact [65],
the current is reformulated as
where is the velocity component normal to the contact. To estimate this integral with the Monte Carlo method (Section 3.3.1) the distribution function needs to be normalized [61]. We assume that the doping concentration under the contact is constant:
This results in a constant electron concentration under the contact () and, therefore, the boundary distribution is independent of the -coordinate:
The electron concentration in general is defined as
With the definition of the normalized distribution function
equations (3.66) and (3.69) become
The normalization constant can be eliminated from this system of equations:
The integrals in the numerator and denominator represent expectation values, that can be estimated by sample means, see Section 3.3.1.
Here, , where is the number of injected particles (), is the number of absorbed particles (), is the constant doping concentration under the contact and is the area of the contact as shown in Figure 3.8.
To estimate the error of the contact current, an estimator is introduced:
where the index corresponds to one particular trajectory. The current is proportional to the sample mean of the , and the statistical error of the current is proportional to the sample variance:
The standard deviation of the current is estimated as
The relative standard deviation can be used as a measure for the statistical error. These equations are applicable because the starting points of all trajectories are statistically independent, and so are the random numbers .