« PreviousUpNext »Contents
Previous: 3.2 Spherical Harmonic Expansion    Top: 3 Transport Modeling Approaches    Next: 4 Backward Monte Carlo Method

Investigating Hot-Carrier Effects using the Backward Monte Carlo Method

3.3 Monte Carlo Method

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 \( \vec {r} \)- and \( \vec {k} \)-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 \( 1/\sqrt {N} \), where \( N \) is the number of random events. In other words, if one wants to reduce the error by a factor \( N \), the calculation cost will increase by \( N^2 \).

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 \( f(\vec {r},\vec {k},t) \). Because of these advantages the Monte Carlo method is used in many cases as a reference method for simpler transport models [60].

3.3.1 Monte Carlo Integration

As mentioned above, the Monte Carlo method is in general a stochastic method to solve integrals. A integral over a function \( \varphi (x) \),

(3.33) \{begin}{align} I=\int \limits _a^b \varphi (x) \dint x \{end}{align}

can be expressed as an expectation value of some random variable. For this purpose, the function \( \varphi (x) \) is factorized:

(3.34) \{begin}{align} \varphi (x) = f(x)\,p(x)\„\qquad I= \int \limits _a^b f(x)\, p(x) \dint x\„ \label {eq:MC-Int} \{end}{align}

where \( p(x) \) is a density function of the Monte Carlo samples satisfying

(3.35) \{begin}{align} p(x)\geq 0,& &\int \limits _a^b p(x) \dint x =1 .         \{end}{align}

The integral in (3.34) represents an expectation value \( I= E\{f\} \). The Monte Carlo method estimates the expectation value by a sample mean:

(3.36) \{begin}{align} I\approx \frac {1}{N}\,\sum \limits _{i=1}^N f(x_i), \{end}{align}

where N is the number of sampling points [61].

3.3.2 Generation of Random Numbers

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].

Inversion Method

The uniformly distributed pseudo-random number \( r \), used in this method, has the property [51]

(3.37) \{begin}{align} 0 \leq r < 1 .          \label {eq:rand_prop} \{end}{align}

For the generation of a random number with a given probability density \( p(x) \), its cumulative distribution function is needed [9, 46, 51, 52]:

(3.38) \{begin}{align} P(x)=\int \limits _{-\infty }^{x}p(x’)\dint x’ , \{end}{align}

with

(3.39) \{begin}{align} P(-\infty )=0, \hspace {2cm} 0\leq P(x) \leq 1, \hspace {2cm}& P(\infty )=1.   \{end}{align}

A \( p(x) \) distributed random number \( x \) can be calculated through the inverse of the cumulative distribution function

(3.40) \{begin}{align} x=P^{-1}(r), \{end}{align}

where \( r \) is a pseudo-random number of type (3.37). However, this method is only applicable if the analytic integral of the function \( p(x) \) is possible to evaluate. Otherwise, different approaches to generate \( p(x) \) distributed pseudo-random numbers are needed.

Rejection Method

For the generation of a random number with a given probability density \( p(x) \), where \( p(x) \) is bounded in the finite interval \( (x_{\text {min}},x_{\text {max}}) \), a uniform distributed random number \( r \) with the properties

(3.41) \{begin}{align} 0 \leq r < \text {max}[p(x)] \{end}{align}

is chosen for the sample point \( x \). If the random number \( r \) lies in the range of

(3.42) \{begin}{align} 0 \leq r < p(x) \{end}{align}

it gets accepted. Else the random number gets rejected [12, 46]. The accepted random number is \( p(x) \) distributed. This method is always applicable, with any bounded function \( p(x) \) in a finite interval. However, if the function \( p(x) \) is heavily peaked, many rejections will occur and the computational expenses will become high.

Combined Method

To keep the rejection rate small, the function which represents the maximum of the probability density, \( \text {max}[p(x)] \), could be replaced by an analytically integrate-able function \( g(x) \). In the interval \( (x_{\text {min}},x_{\text {max}}) \), \( g(x) \) must be always greater than \( p(x) \).

Here, a random number \( r \) can be directly obtained from \( g(x) \) with the inversion method. This random number \( r \) can than be applied to the rejection method, where the number of rejections can be significantly reduced by finding a suitable function \( g(x) \). In the case, that \( g(x) \) is chosen as a constant, this method will result in the rejection method, described above.

3.3.3 Duration of the Free Flight

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 \( \lambda (t) \). The probability for the scattering of a carrier in an interval \( \Delta t \) at a time \( t \) is \( \lambda (t)\,\Delta t \), where \( \lambda (t) \) represents the \( \vec {r} \) and \( \vec {k} \) dependent scattering rate at the time \( t \):

(3.43) \{begin}{align} \lambda (t)=\lambda (\vec {k}{(t)},\vec {r}{(t)}) \{end}{align}

With the assumption, that the carrier scatters at \( t=0 \), the cumulative distribution function \( P(t) \) can be calculated. In this case, no other scattering process takes place until the time \( t \). Thus, the carrier is in free-flight for the duration of \( t \) [60]:

(3.44) \{begin}{align} P(t)=\begin {cases} 1-\e ^{-\int \limits _{0}^{t}\lambda (\tau )\dint \tau } &\text {:   } t \geq 0\\ 0 &\text {:   }t < 0 \end {cases}.
\{end}{align}

Using the inversion method, the duration of the free-flight \( t_f \) can be calculated from the equation \( P(t_f)=r \) as

(3.45) \{begin}{align} \int \limits _{0}^{t_f}\lambda (\vec {k}{(t)},\vec {r}{(t)})\dint t =-\ln (1-r), \{end}{align}

where \( r \) represents a uniform distributed random number [46]. The integration takes place along the trajectory of the carrier \( (\vec {k}{(t)},\vec {r}{(t)}) \), which can be acquired by integrating the equations of motion (2.22) and (2.23).

The full-band simulations use constant total scattering rates \( \Gamma _\mathrm {max} \), as shown in Fig. 3.1 [110]. In case \( \Gamma \) is constant, the duration of the free-flight is [P3]:

(3.46) \{begin}{align} t_f=-\frac {1}{\Gamma }\ln (1-r).   \{end}{align}

Because \( \Gamma _\mathrm {max} \) must be larger than the sum of the scattering rates of all physical scattering processes \( \lambda (t) \), the self-scattering is more likely to occur in areas where \( \lambda (t) \) is significant smaller than \( \Gamma _\mathrm {max} \). For this reason \( \lambda (t) \) can be approximated by local, picewise constant values \( \Gamma _\mathrm {max}(t) \). Since these local values are smaller than the global \( \Gamma _\mathrm {max} \), self-scattering is reduced. The duration of a collision-less free-flight is given by [51]

(3.47) \{begin}{align} \Gamma _{\mathrm {max}}(t_f - t_j) = -\ln (1-r) - \sum \limits _{k=1}^j \Gamma _{\mathrm {max}}(t_k -t_{k-1}) \qquad \mathrm {for}\qquad t_j < t_f <
t_{j+1}\„ \{end}{align}

where \( t_j \) is the time at which the particle passes a change in \( \Gamma _\mathrm {max}(t) \) during the free-flight \( t_f \).

(-tikz- diagram)

Figure 3.1: The scattering rate \( \lambda (t) \), the picewise constant scattering rate \( \Gamma (t) \) and the global maximum \( \Gamma _\mathrm {max} \).

3.3.4 Selection of the Scattering Process

After the free-flight, another uniformly distributed random number,

(3.48) \{begin}{align} 0 \leq r < \Gamma _\mathrm {max} , \{end}{align}

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 \( \Gamma _\mathrm {max} \),

(3.49) \{begin}{align} \Gamma _\mathrm {max} \geq \sum \limits _{i=1}^{N}\lambda _i\left (t_f\right ), \{end}{align}

where \( \lambda _i \) is the rate of the \( i^\mathrm {th} \) scattering process. The scattering process \( m \) is chosen when

(3.50) \{begin}{align} \sum \limits _{i=1}^{m-1}\lambda _i\left (t_f\right ) \leq r < \sum \limits _{i=1}^{m}\lambda _i\left (t_f\right )\,.   \{end}{align}

If the random number is in the range of

(3.51) \{begin}{align} \Gamma _\mathrm {max} \geq r > \sum \limits _{i=1}^{N}\lambda _i\left (t_f\right ), \{end}{align}

self scattering occurs [60].

3.3.5 Generation of an Equilibrium Distribution

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.

Box Method

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]

(3.52) \{begin}{align} f_v(\vec {k})=|v_\perp (\vec {k})| \, f(\vec {k})\„ \{end}{align}

where \( |v_\perp (\vec {k})| \) is the velocity component perpendicular to the boundary.

In the statistical average of an attribute \( A(\vec {k}) \) the weighting factor \( |v_\perp (\vec {k})|^{-1} \) is obtained by replacing \( f(\vec {k}) \) by \( f_v(\vec {k})\,|v_\perp (\vec {k})|^{-1} \) [P3]:

(3.53) \{begin}{align} \braket {A}_\mathrm {Box}=&\frac {\int f(\vec {k}) A(\vec {k})\dk }{\int f(\vec {k}) \, \dk } =\frac {\int {f_v(\vec {k})}\,|v_\perp (\vec {k})|^{-1}
A(\vec {k}) \, \dk }{\int {f_v(\vec {k})}\,|v_\perp (\vec {k})|^{-1}\dk } \approx \frac {\sum \limits _{i=1}^{N} {A(\vec {k}_i)}\,|v_\perp (\vec {k}_i)|^{-1}}{\sum \limits _{i=1}^{N} |v_\perp
(\vec {k}_i)|^{-1}}, \{end}{align}

where \( N \) represents the number of sampled states.

(-tikz- diagram)

Figure 3.2: Sampling using the box method.

Bulk Method

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 \( \vec {k} \)-values is a Maxwellian distribution \( f(\vec {k}) \) weighted with the scattering rate \( \Gamma (\vec {k}) \) [46, 60].

(3.54) \{begin}{align} f_\Gamma (\vec {k})=\Gamma (\vec {k}) \, f(\vec {k}), \{end}{align}

When the mean value of an attribute \( A(\vec {k}) \) is calculated, the weighting factor \( \Gamma ^{-1} \) has to be taken into account:

(3.55) \{begin}{align} \braket {A}_\mathrm {Bulk}=&\frac {\int f(\vec {k}) A(\vec {k}) \, \dk }{\int f(\vec {k})\dk } =\frac {\int {f_\Gamma (\vec {k})}\,{\Gamma (\vec {k})}^{-1}
A(\vec {k})\dk }{\int {f_\Gamma (\vec {k})}\,{\Gamma (\vec {k})}^{-1}\dk } \approx \frac {\sum \limits _{i=1}^{N} {A(\vec {k}_i)}\,{\Gamma (\vec {k}_i)}^{-1}}{\sum \limits _{i=1}^{N} {\Gamma
(\vec {k}_i)}^{-1}} . \label {eq:before-scattering-sampeling} \{end}{align}

(-tikz- diagram)

Figure 3.3: Sampling using the before-scattering method.

Constant Time Sampling Method

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.

(3.56) \{begin}{align} \braket {A}_\mathrm {Time}=&\frac {\int f(\vec {k}) A(\vec {k}) \, \dk }{\int f(\vec {k})\dk } \approx \frac {1}{N}\,\sum \limits _{i=1}^{N} A(\vec {k}_i).
\{end}{align}

(-tikz- diagram)

Figure 3.4: Constant time sampling method.

3.3.6 Calculation of Equilibrium Averages

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 \( A \) is calculated like [74]

(3.57) \{begin}{align} \label {eq:estimator} \braket {A}=\frac {\sum \limits _n \int \limits _\mathrm {BZ} A_n(\vec {k})\e ^{-\beta \Epsilon _n(\vec {k})}\dk }{\sum \limits _n \int
\limits _\mathrm {BZ} \e ^{-\beta \Epsilon _n(\vec {k})}\dk }, \{end}{align}

where BZ denotes the whole Brillouin zone, \( n \) is the band index, \( \beta = \left (k_\mathrm {B}T\right )^{-1} \) with the Boltzmann constant \( k_\mathrm {B} \) and the Temperature \( T \). In statistical mechanics a parameter \( Z \) called partition function is introduced [74] .

(3.58) \{begin}{align} \label {eq:partfunc} Z=\sum \limits _n \int \limits _\mathrm {BZ} \e ^{-\beta \Epsilon _n(\vec {k})}\dk \{end}{align}

With the partition function (3.58) the statistical average (3.57) becomes:

(3.59) \{begin}{align} \label {eq:est_short} \braket {A}=\frac {1}{Z}\sum \limits _n \int \limits _\mathrm {BZ} A_n(\vec {k})\e ^{-\beta \Epsilon _n(\vec {k})}\dk .   \{end}{align}

For a detailed description of the numerical integration see Appendix A.1 or [P3].

3.3.7 Estimation of Non-Equilibrium Averages

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.

Before scattering method

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 \( \Gamma (\vec {k}) \), as shown in (3.55) [51, 60]:

In a spatial discretization every grid point \( P_j \) in a device has a volume \( V_j \) assigned. Inside this Volume \( V_j \), the closest grid point is \( P_j \), where \( j \) is an index for every grid point in \( \vec {r} \) 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 \( V_j \) separately in the manner of [88]:

(3.60) \{begin}{align} \braket {A(\vec {k})}_j \approx \frac {\sum \limits _{\vec {r}_{i}\in V_j} {A(\vec {k}_i)}\,{\Gamma (\vec {k}_i)}^{-1}}{\sum \limits _{\vec {r}_{i}\in V_j}
{\Gamma (\vec {k}_i)}^{-1}} . \label {eq:spartial-bs-sampling} \{end}{align}

Here, the summation runs over the before scattering states \( (\vec {k}_i,\vec {r}_i) \), where the scattering takes place inside the Volume \( V_j \) [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).

Boundary method

The boundary method gathers statistical information when a particle crosses a boundary from one Voronoi volume to another. Therefore, the \( \vec {k} \)-values are weighted with the velocity component of the particle perpendicular to the boundary \( |v_\perp (\vec {k})| \) [65]. In the manner of the box method, the mean values of the local attributes can be calculated like:

(3.61) \{begin}{align} \braket {A(\vec {k})}_j \approx \frac {\sum \limits _{\vec {r}_{n}\in e_j} {A(\vec {k}_n)}\,|v_\perp (\vec {k}_n)|^{-1}}{\sum \limits _{\vec {r}_{n}\in e_j}
|v_\perp (\vec {k}_n)|^{-1}}, \{end}{align}

where \( \vec {r}_n \) is the point on the \( j \)-th edge, where one particle is crossing the \( n \)-th time the edge \( e_j \) 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.

(-tikz- diagram)

Figure 3.5: Boundary method sampling in a device. The dashed lines are the boundaries between the Voronoi volumes.

Combination of methods

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 \( n^+n^-n^+ \) diode with abrupt junctions was chosen. The doping levels are \( 10^{19} \)cm\( ^{-3} \) and \( 10^{15} \)cm\( ^{-3} \), respectively. Fig. 3.6 shows the conduction band edge and the electron density for an applied voltage of \( \SI {2}{V} \).

Fig. 3.7 compares the different averaging methods in the \( n^+n^-n^+ \) diode in the area of the \( n^+n^- \) 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 \( x > \SI {250}{\nano \meter } \) 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.

(image)

Figure 3.6: Conduction band edge (red) and the electron density (blue) in an \( n^+n^-n^+ \) diode with abrupt junctions.

(image)

(a) Constant grid-cell size. The number of entries on the edge is always higher than in the Voronoi volume (before-scattering method)

(image)

(b) Variable gird-cell size. For \( x < \SI {250}{nm} \), the number of entries on the edges is higher. For \( x > \SI {250}{nm} \), the number of entries in the Voronoi volumes is higher. This results from different grid-cell sizes.

Figure 3.7: Comparison of two different averaging methods in an silicon \( n^+n^-n^+ \) diode with constant and variable grid size. Both images show the region around the \( n^+n^- \) junction at 200 nm.

3.3.8 Current Estimation

(-tikz- diagram)

Figure 3.8: Sketch of a device contact in the \( x \)-\( z \) plane.

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:

(3.62) \{begin}{align} \vec {J}=\frac {2\,q}{\left (2\pi \right )^3} \, \int \limits _\mathrm {BZ} \vec {v}(\vec {k}) f(\vec {k},\vec {r}) \dk , \{end}{align}

where \( f(\vec {k},\vec {r}) \) is the unknown solution of the BTE and \( q \) is the particle charge. The current through a contact in a two-dimensional device can be calculated as

(3.63–3.64) \{begin}{align} I=&\int \limits _\mathrm {contact} \vec {J} \cdot \dA = \int \limits _{z=0}^W\int \limits _{x=0}^{L_\mathrm {C}} \vec {J}\cdot \vec {e}_y \dint x
\dint z = W\,\int \limits _{x=0}^{L_\mathrm {C}} J_y(x,0) \dint x \\ =&\frac {q}{4\pi ^3}\,W \int \limits _{x=0}^{L_\mathrm {C}} \int \limits _\mathrm {BZ} v_y(\vec {k}) f(\vec {k},x,0)
\dk \dint x , \label {eq:current-fmc} \{end}{align}

where \( \int \dA =\int \vec {e}_y\,W \dint x \) represents the integration over the contact area and \( W \) is the width of the contact, as shown in Figure 3.8. Introducing the velocity weighted distribution of the contact [65],

(3.65) \{begin}{align} f_v(\vec {k},\vec {r})=|v_y(\vec {k})| f(\vec {k},\vec {r}) , \{end}{align}

the current is reformulated as

(3.66) \{begin}{align} \label {eq:MC_Iest} I=&\frac {q}{4\pi ^3}\,W \int \limits _{x=0}^{L_\mathrm {C}} \int \limits _\mathrm {BZ} \mathrm {sign}\left (v_y(\vec {k}) \right )
f_v(\vec {k},x,0) \dk \dint x , \{end}{align}

where \( v_y \) is the velocity component normal to the contact. To estimate this integral with the Monte Carlo method (Section 3.3.1) the distribution function \( f_v(\vec {k},x,0) \) needs to be normalized [61]. We assume that the doping concentration under the contact is constant:

(3.67) \{begin}{align} N_\mathrm {D}(x,y=0) = N_\mathrm {surf} .    \{end}{align}

This results in a constant electron concentration under the contact (\( n_\mathrm {surf}=N_\mathrm {surf} \)) and, therefore, the boundary distribution is independent of the \( x \)-coordinate:

(3.68) \{begin}{align} f_v(\vec {k},x,0)=f_v^\mathrm {surf}(\vec {k}) \{end}{align}

The electron concentration in general is defined as

(3.69) \{begin}{align} n(\vec {r})=\frac {1}{4\pi ^3}\int \limits _\mathrm {BZ} f(\vec {k},\vec {r}) \dk =\frac {1}{4\pi ^3}\int \limits _\mathrm {BZ} \frac {1}{|v_y|} f_v(\vec
{k},\vec {r}) \dk . \label {eq:nsurf} \{end}{align}

With the definition of the normalized distribution function

(3.70) \{begin}{align} \label {eq:MC_psurf} p_v^\mathrm {surf}(\vec {k})=\frac {f_v^\mathrm {surf}(\vec {k})}{\mathcal {C}} , \{end}{align}

equations (3.66) and (3.69) become

(3.71–3.72) \{begin}{align} I=&\frac {q}{4\pi ^3}\,W\,L_\mathrm {C}\,\mathcal {C} \int \limits _\mathrm {BZ} \mathrm {sign}\left (\vec {v}_y(\vec {k}) \right ) p_v^\mathrm
{surf}(\vec {k}) \dk \„\\ n_\mathrm {surf}=&\frac {\mathcal {C}}{4\pi ^3}\int \limits _\mathrm {BZ} \frac {1}{|v_y|} p_v^\mathrm {surf}(\vec {k}) \dk \,. \{end}{align}

The normalization constant \( \mathcal {C} \) can be eliminated from this system of equations:

(3.73) \{begin}{align} I= q\,W\,L_\mathrm {C}\,n_\mathrm {surf} \frac {\int \limits _\mathrm {BZ} \mathrm {sign}(v_y)\,p_\mathrm {surf}(\vec {k}) \dk }{\int \limits _\mathrm {BZ}
\frac {1}{|v_y|}\,p_\mathrm {surf}(\vec {k}) \dk } . \{end}{align}

The integrals in the numerator and denominator represent expectation values, that can be estimated by sample means, see Section 3.3.1.

(3.74) \{begin}{align} I&= q\,W\,L_\mathrm {C}\,N_\mathrm {surf} \frac {\sum \limits _{i=1}^{N} \mathrm {sign}(v_y)}{\sum \limits _{i=1}^{N}\frac {1}{|v_y|}} = q\,W\,L_\mathrm
{C}\,n_\mathrm {surf} \frac {N_\mathrm {in}-N_\mathrm {out}}{\sum \limits _{i=1}^{N}\frac {1}{|v_y|}} \,. \{end}{align}

Here, \( N=N_\mathrm {in} + N_\mathrm {out} \), where \( N_\mathrm {in} \) is the number of injected particles (\( \mathrm {sign}(v_y)=1 \)), \( N_\mathrm {out} \) is the number of absorbed particles (\( \mathrm {sign}(v_y)=-1 \)), \( N_\mathrm {surf} \) is the constant doping concentration under the contact and \( W\,L_\mathrm {C} \) is the area of the contact as shown in Figure 3.8.
To estimate the error of the contact current, an estimator \( \nu _i \) is introduced:

(3.75) \{begin}{align} \nu _i=\begin {cases} -1 &\text {: particle gets absorbed but not injected at this contact}\\ 0 &\text {:   particle gets injected and absorbed at this
contact}\\ +1 &\text {: particle gets injected but not absorbed at this contact} \end {cases}, \{end}{align}

where the index \( i \) corresponds to one particular trajectory. The current is proportional to the sample mean of the \( \nu _i \), and the statistical error of the current is proportional to the sample variance:

(3.76) \begin{equation} s_\nu ^2 = \frac {1}{N-1}\biggl (\sum \limits _{i=1}^N \nu _i^2 - N\,\bar \nu ^2\biggr ) \end{equation}

The standard deviation of the current is estimated as

(3.77) \begin{equation} s_I = \frac {s_\nu }{\sqrt {N}}\quad \Rightarrow \quad s_\mathrm {rel}=\frac {s_I}{I}\,.   \end{equation}

The relative standard deviation \( s_\mathrm {rel} \) 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 \( \nu _i \).

« PreviousUpNext »Contents
Previous: 3.2 Spherical Harmonic Expansion    Top: 3 Transport Modeling Approaches    Next: 4 Backward Monte Carlo Method