« PreviousUpNext »Contents
Previous: 4 Backward Monte Carlo Method    Top: 4 Backward Monte Carlo Method    Next: 4.3 Multi-Band Semiconductors

Investigating Hot-Carrier Effects using the Backward Monte Carlo Method

4.2 Current Estimators

The main difference between the current definition of the FMC method (3.64) and the BMC method is, that the BMC method is calculating the current not through a contact but through an area \( y_\mathrm {m}\times W \) in the \( yz \)-plane located at \( x=x_0 \) [P6]:

(4.24) \{begin}{align} I =\frac {q}{4\pi ^3}\,W \int \limits _{0}^{y_\mathrm {m}} \int \limits _\mathrm {BZ} v_x(\vec {k}_0) \, f(\vec {k}_0,x_0,y_0) \dk _0 \dint y_0 \„ \label {eq:current-1} \{end}{align}

In principle, the integral in (4.24) could be evaluated by numerical integration, whereby the values of the distribution function at the discrete points \( (\kv _0,\rv _0) \) are estimated by (4.23). However, it is more convenient to employ Monte Carlo integration instead. For this purpose, the current has to be expressed as an expectation value. This is accomplished by introducing a PDF \( p_0(\kv _0,y_0) \) which can be chosen freely, and reformulating (4.24) as [P4]:

(4.25) \{begin}{align} I &= q W\!\int \limits _{0}^{y_\mathrm {m}}\!\int \limits _\mathrm {BZ}\! \mu (\kv _0,y_0; x_0) \,p_0(\kv _0,y_0) \dk _0 \dint y_0 \equiv q W E\{\mu \} \{end}{align}

with

(4.26) \{begin}{align} \mu (\kv _0,y_0; x_0) = \frac {v_x(\kv _0) \, f(\vec {k}_0,x_0,y_0)}{4\pi ^3\, p_0(\kv _0,y_0)}. \label {eq:BMC-current-est} \{end}{align}

In a Monte Carlo simulation, the expectation value is estimated by a sample mean [83].

(4.27) \{begin}{align} I = q W E\{\mu \} \approx q\,W\,\frac {1}{N}\sum \limits _{i=1}^N \mu (\vec {k}_{0,i},y_{0,i}; x_0) \label {eq:sample-mean-I} \{end}{align}

Here, \( N \) denotes the number of sampling points. With the estimated distribution function (4.23) the estimator (4.26) takes the form

(4.28) \begin{equation} \mu (\kv _0,y_0; x_0) = \frac {v_x(\kv _0) \, f_b(\vec {k}_b,\rv _b)}{4\pi ^3\, p_0(\kv _0,y_0)} \e ^{\,\beta _D (H(\kv _{b},\rv _{b}) - H(\kv _0,\rv _0))}. \label {eqn:nu-definition-2} \end{equation}

Consider, that with this definition only one trajectory per sampling point \( (\kv _{0,i}, y_{0,i}) \) is started, corresponding to \( M=1 \) in (4.23).

The estimator (4.28) is the fundamental equation from which various variants can be derived. The following sections discuss different choices of the injection PDF \( p_0 \) and the properties of the resulting current estimators [P4]. In (4.28), \( \kv _0 \) and \( y_0 \) are random variables, whereas \( x_0 \) is a given parameter. In the following, \( x_0 \) is omitted from the argument list for the sake of readability.

4.2.1 The Boundary Distribution

Distribution functions at ohmic contacts are close to thermal equilibrium. Therefore, a Maxwell-Boltzmann or Fermi-Dirac distribution is an appropriate choice for the the boundary distribution \( f_b \). Here, an equilibrium Boltzmann distribution is assumed:

(4.29) \{begin}{align} f_b(\kv ,\rv )= \mathcal {C}(\rv )\, \e ^{-\beta _\mathrm {D} \, \Epsilon (\kv )} \label {eq:f_b} \{end}{align}

Two normalization integrals will be needed in the following. The first one is the partition function \( Z(T) \) (3.58) and the second integral is defined as

(4.30) \begin{equation} V(T)= \int \limits _\mathrm {BZ}|v_x(\kv _0)|\,\e ^{-\beta (T)\,\Epsilon (\kv _0)} \dk \,. \label {eq:V} \end{equation}

The integrations over the Brillouin zone (BZ) are carried out numerically, see Appendix A.1. From these two quantities the injection velocity is obtained:

(4.31) \begin{equation} v_\mathrm {inj}(T) = \frac {V(T)}{Z(T)} \label {eq:injection-velocity} \end{equation}

The definition of the electron concentration

(4.32) \{begin}{align} n(\rv )=&\frac {1}{4\pi ^3}\int \limits _\mathrm {BZ} f_b(\kv ,\rv ) \dk \label {eq:n} \{end}{align}

is used to determine the normalization constant \( \mathcal {C} \) in (4.29):

(4.33) \begin{equation} \mathcal {C}(\rv ) = \frac {4\pi ^3 n(\rv )}{Z(T_\mathrm {D})} \label {eq:C} \end{equation}

4.2.2 Injection from an Equilibrium Maxwellian

The starting points of the backward trajectories are generated by the PDF \( p_0 \). We express \( p_0 \) as a product of two independent PDFs:

(4.34) \{begin}{align} p_0(\vec {k}_0,y_0) =& f_0(\vec {k}_0)\, p_y(y_0) \label {eq:p_0} \{end}{align}

The PDF of the injection coordianate \( y_0 \) is assumed to be proportional to the electron concentration at the injection coordinate \( x_0 \).

(4.35) \begin{equation} p_y(y_0) = \frac {n(x_0, y_0)}{\int _0^{y_\mathrm {m}} n(x_0, y)\,\d y} \label {eq:p_y} \end{equation}

For the injection distribution \( f_0 \), a normalized Boltzmann distribution at device temperature \( T_D \) is chosen.

(4.36) \{begin}{align} f_0(\vec {k}_0) =& \frac {1}{Z(T_D)}\, \e ^{-\beta _D\, \Epsilon (\vec {k}_0)} \label {eq:equil-maxwell} \{end}{align}

Inserting the boundary distribution (4.29) and the injection distribution (4.36) in (4.28) gives the following current estimator [P4]:

(4.37) \{begin}{align} \label {eq:est1} \mu = v_x(\kv _{0}) \frac { \e ^{ \beta _\mathrm {D}\,\left (E_\mathrm {C}(\vec {r}_{\mathrm {b}}) - E_\mathrm {C}(\vec {r}_{0})\right ) } } { p_y(y_{0}) } n(\vec {r}_{\mathrm {b}}) \{end}{align}

Note that both the Boltzmann factors \( \e ^{-\beta _D\,\Epsilon (\kv _0)} \) and \( \e ^{-\beta _D\,\Epsilon (\kv _b)} \) have canceled out of this expression. To generate wave vectors from the equilibrium distribution (4.36), the constant time sampling method described in Section 3.3.5 can be used.

4.2.3 Injection from a Velocity-weighted Maxwellian

Another choice for \( f_0 \) is a velocity-weighted Maxwellian at equilibrium temperature \( T_D \) [P4]:

(4.38) \{begin}{align} f_0(\kv _0) =& \frac {1}{V(T_D)}\, |v_x(\kv _0)|\,\e ^{-\beta _D \,\Epsilon (\kv _0)} \label {eq:vel-weighted-df} \{end}{align}

This choice is motivated by the fact that in the numerator of (4.28) a term \( v_x(k_0)\, \e ^{-\beta _D\,\Epsilon (k_0)} \) occurs. Division by (4.38) will essentially cancel out this term. This reduces the \( k_0 \)-dependence of the estimator which is expected to reduce its variance. Inserting the boundary distribution (4.29) and the injection distribution (4.38) in (4.28) yields

(4.39) \{begin}{align} \label {eq:Iest} \mu = \mathrm {sign}(v_x(\kv _0))\, v_\mathrm {inj}(T_D) \frac {\e ^{\,\beta _\mathrm {D}\,\left (E_\mathrm {C}(\rv _b) - E_\mathrm {C}(\rv _0)\right )}}{p_y(y_0)} n(\rv _b)\,. \{end}{align}

In this equation, \( \mathrm {sign}(v_x) \) denotes the sign of the velocity component \( v_x \), and \( v_\mathrm {inj} \) is the injection velocity defined by (4.31). To generate wave vectors from the equilibrium distribution (4.38), the box sampling method described in Section 3.3.5 can be used.

4.2.4 Injection from a Non-equilibrium Maxwellian

For some applications, it can be useful to generate the initial points \( \kv _0 \) from a Maxwellian at a temperature \( T_0 \) different from the device temperature \( T_D \) [P4]. When calculating quantities depending on the high energy tail of the distribution, an injection temperature \( T_0 > T_D \) will be beneficial as it enhances the number of initial points at higher energies. In this work a non-equilibrium Maxwellian of the form

(4.40) \begin{equation} f_0(\kv _0) = \frac {1}{Z(T_0)} \e ^{-\beta _0\, \Epsilon (\kv _0)} \label {eq:heated-maxwellian} \end{equation}

is considered. Injecting with this non-equilibrium distribution, equation (4.28) leads to an estimator \( \eta \):

(4.41) \{begin}{align} \eta (\kv _0,y_0) = \frac {Z(T_0)}{Z(T_\mathrm {D})} \,v_x(\kv _0)\, \frac {\e ^{\,\beta _D\,\left (E_\mathrm {C}(\rv _b) - E_\mathrm {C}(\rv _0)\right )}}{p_y(y_0)} \, \times \e ^{\left (\beta _0-\beta _\mathrm {D}\right )\,\Epsilon (\vec {k}_0)}\,n(\rv _b) \label {eq:nu-heated-Maxwellian} \{end}{align}

The sample mean of \( \eta \),

(4.42) \begin{equation} I = q W\, \frac {1}{N} \sum \limits _{i=1}^N \eta (\kv _{0,i}, y_{0,i}) \„ \end{equation}

can be reformulated as a weighted average of the form:

(4.43) \begin{equation} I = q W \,\frac {\sum \limits _{i=1}^N \mu (\kv _{0,i}, y_{0,i})\, w(\kv _{0,i})}{\sum \limits _{i=1}^N w(\kv _{0,i})} \label {eq:weighted-average} \end{equation}

Here, \( \mu \) is given by (4.37) and the weight factor \( w \) is defined as

(4.44) \begin{equation} w(k_0) = \e ^{(\beta _0-\beta _D)\, \Epsilon (\kv _0)}\,. \label {eq:definition-w} \end{equation}

The derivation of (4.43) is described in detail in [P4]. In the case of a velocity-weighted Maxwellian with \( T_0\ne T_D \),

\begin{equation*} f_0(\kv _0) = \frac {1}{V(T_0)} |v_x(\kv _0)| \,\e ^{-\beta _0\, \Epsilon (\kv _0)}\, \end{equation*}

a similar procedure can be applied. In the weighted average (4.43) the estimator (4.39) and the very same weight (4.44) have to be used [P4].

4.2.5 Injection from the Equilibrium Concentration

All previous estimators have a dependence on the injection coordinate \( y_0 \) through the term \( \e ^{\,\beta _D E_\mathrm {C}(\rv _0)}/p_y(y_0) \). However, this dependency is weak and can even be eliminated from the estimator by choosing the injection distribution as [P4]:

(4.45) \{begin}{align} p_0(\kv _0, y_0; x_0) = \frac {1}{A(T_D)} |v_x(\kv _0)| \e ^{-\beta _D H(\kv _0,x_0, y_0)} \label {eq:p_0-H} \{end}{align}

Inserting (4.20) in (4.45) again yields a product of two independent PDFs.

(4.46) \{begin}{align} p_0(\vec {k}_0,y_0) =& f_0(\vec {k}_0)\, \tilde p_y(y_0) \{end}{align}

Here, \( f_0 \) is given by (4.38), and \( \tilde p_y \) is defined as

(4.47) \{begin}{align} \tilde p_y(y_0) = \frac {\tilde n(x_0, y_0)}{B(T_D)} \label {eq:tilde-p_y} \{end}{align}

with \( \tilde n(x_0,y_0) = \e ^{-\beta _D E_C(x_0, y_0)} \). This quantity is up to a constant the equilibrium concentration determined by the band edge energy \( E_C \). On the other hand, \( n \) in (4.35) represents the actual carrier concentration as obtained from a device simulation. With the normalization integral in (4.47) defined as

(4.48) \{begin}{align} B(T_D) = \int _0^{y_\mathrm {m}} \tilde n(x_0, y)\,\d y\„ \{end}{align}

the total normalization factor \( A \) in (4.45) becomes

(4.49) \{begin}{align} A(T_D) = V(T_D)\, B(T_D)\,. \{end}{align}

Using the boundary distribution (4.29) and the injection distribution (4.45) the current estimator (4.28) can be reformulated as:

(4.50) \{begin}{align} \mu = \mathrm {sign}(v_x(\kv _0))\, v_\mathrm {inj}(T_D)\, B(T_D) \,n(\rv _b)\,\e ^{\,\beta _D E_C(\rv _b)} \label {eq:nu-4} \{end}{align}

Other than the estimators discussed above, this estimator is independent of the injection coordinate \( y_0 \). A more transparent physical interpretation is achieved by expressing the equilibrium concentration \( n(\rv _b) \) as a function of the local quasi-Fermi level \( F_n \).

(4.51) \{begin}{align} n(\rv _b) = N_C(T_D)\,\e ^{\,\beta _D(F_n(\rv _b) - E_C(\rv _b)} \label {eq:equilibrium-conc} \{end}{align}

Here, the effective density of states \( N_C \) is related to the partition function \( Z \) by \( N_C = Z/(4\pi ^3) \). Also, the normalization factor \( B \) will be expressed through an energy \( \mean {E}_C \) defined as

(4.52) \{begin}{align} \mean {E}_C = -k_B T_D\mathrm {ln}\left (\frac {B}{y_\mathrm {m}}\right )\,. \{end}{align}

\( \mean {E}_C \) has the meaning of an average of the band edge energy over the injection coordinate \( y_0 \):

(4.53) \{begin}{align} \e ^{-\beta _D\mean {E}_C(x_0)} = \frac {1}{y_\mathrm {m}} \int \limits _0^{y_\mathrm {m}} \,\e ^{-\beta _D E_C(x_0, y_0)}\,\d y_0 \{end}{align}

Expressing the estimator (4.50) in terms of the parameters \( F_n \) and \( \mean {E}_C \) gives

(4.54) \{begin}{align} \mu = y_\mathrm {m} N_C(T_D)\, v_\mathrm {inj}(T_D)\,\mathrm {sign}(v_x(\kv _0)) \,\e ^{\,\beta _D(F_n(\rv _b) - \mean {E}_C)}\,. \label {eq:nu-5} \{end}{align}

This equation states that a backward trajectory represents an elementary particle flux \( N_C\, v_\mathrm {inj} \). This flux is multiplied by a statistical weight given by the e-function. The higher the energy of the starting point (\( \mean {E}_C \)) with respect to the Fermi level at the trajectory end point (\( F_n \)), the lower is the statistical weight. If a constant \( \Delta \Epsilon \) were added to \( \mean {E}_C \), the estimator \( \mu \) and subsequently also the current \( I \) would be scaled by the factor \( \e ^{-\beta _D \Delta \Epsilon } \). In other words, increasing the barrier height by some energy increment will result in an exponential decrease in current. This means that the exponential dependence of the thermionic current on the barrier height can be directly deduced from the current estimator (4.54) [P4].

4.2.6 Symmetric Sampling

In thermodynamic equilibrium the distribution function is symmetric. Because of this symmetry, the current will vanish. In a BMC simulation the current is not vanishing exactly because of the finite sample size. However, this type of statistical error can be eliminated by always generating positive and negative values of the estimator in pairs. When a backward trajectory is started from a state \( (\kv _0, y_0) \), also another one is started with opposite momentum from the state \( (-\kv _0, y_0) \). This procedure will give exactly \( I=0 \) in thermal equilibrium without statistical error and is reducing the statistical error in situations close to thermal equilibrium.

Every estimator described above can be used to define a new estimator by taking the algebraic mean value [P4]:

(4.55) \{begin}{align} \mu _\mathrm {symm}(\kv _0,y_0) =\frac {\mu (\kv _0,y_0) + \mu (-\kv _0,y_0)}{2} \{end}{align}

Using (4.54) the new estimator will be of the form:

(4.56) \{begin}{align} \mu _\mathrm {symm} = \frac {y_\mathrm {m} N_C(T_D)\, v_\mathrm {inj}(T_D)}{2} \left (\e ^{\,\beta _D\,F^+_n(\rv _b)} - \e ^{\,\beta _D\,F^-_n(\rv _b)} \right ) \e ^{-\beta _D \mean {E}_C} \label {eq:nu-symmetric} \{end}{align}

Here, \( F^+_n \) denotes the quasi-Fermi level of the contact where the trajectory injected from \( k_0 \) has terminated, whereas \( F^-_n \) is the quasi-Fermi level of the contact where the trajectory injected from \( -k_0 \) has terminated [P4].

4.2.7 Estimation of the Statistical Error

Because of the statistical independence of the backward trajectories, an expression for the statistical error of the simulation result is readily found. In the following, \( \mu \) is a function of the random variables \( (\kv _0, y_0) \). Several such functions \( \mu (\kv _0, y_0) \) have been discussed in the preceding sections.

In the case of the injection states \( \kv _0 \) being generated from an equilibrium distribution, the sample mean \( \mean \mu \) and the variance \( s_\mu ^2 \) of the samples \( \{\mu _1,\mu _2,\cdots ,\mu _N\} \) can be calculated straightforwardly. The sample mean

(4.57) \begin{equation} \mean \mu = \frac {1}{N}\sum \limits _{i=1}^N \mu _i \end{equation}

gives the current, \( I=qW\mean \mu \), whereas the sample variance \( s_\mu ^2 \) allows an estimate of the current’s statistical error.

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

The standard deviation of the current is estimated as

(4.59) \begin{equation} s_I = qW \frac {s_\mu }{\sqrt {N}}\quad \Rightarrow \quad \frac {s_I}{I} = \frac {1}{\sqrt {N}}\,\frac {s_\mu }{\mean \mu }\,. \end{equation}

The relative standard deviation \( s_I/I \) can be used as a measure for the statistical error.

In the case of the injection states \( \kv _0 \) being generated from a non-equilibrium distribution, the random variable \( w \) defined by (4.44) and the random variable \( \xi \) defined by

(4.60) \begin{equation} \xi = \mu \,w\„ \end{equation}

are needed. In the course of a Monte Carlo simulation, the sample means \( \mean \xi \) and \( \mean {w} \) have to be calculated in order to obtain the current [P4]:

(4.61) \begin{equation} \mean \xi = \frac {1}{N}\sum \limits _{i=1}^N \xi _i,\qquad \mean {w} = \frac {1}{N}\sum \limits _{i=1}^N w_i \quad \Rightarrow \quad I = qW \frac {\mean \xi }{\mean {w}} \end{equation}

In addition, the sample variances and the sample covariance have to be determined.

(4.62–4.64) \{begin}{align} s_\xi ^2 &= \frac {1}{N-1}\biggl (\sum \limits _{i=1}^N \xi _i^2 - N\,\mean \xi ^2\biggr ) \\ s_w^2 &= \frac {1}{N-1}\biggl (\sum \limits _{i=1}^N w_i^2 - N\,\mean {w}^2\biggr ) \\ s_{\xi w}^2 &= \frac {1}{N-1}\biggl (\sum \limits _{i=1}^N \xi _i\,w_i - N\,\mean \xi \,\mean {w}\biggr ) \{end}{align}

Using these parameters, the variance of the random variable \( \mu = \xi / w \) can be estimated as

(4.65) \begin{equation} s_\mu ^2 = s_\xi ^2 - 2 r s_{\xi w}^2 + r^2 s_w^2\„ \end{equation}

where \( r=\mean \xi /\mean w \) [83]. From \( s_\mu \), the standard deviation of the current can be computed.

(4.66) \begin{equation} s_I = qW \frac {s_\mu }{\sqrt {N}\,\mean {w}} \quad \Rightarrow \quad \frac {s_I}{I} = \frac {1}{\sqrt {N}}\,\frac {s_\mu }{\mean \xi } \end{equation}

« PreviousUpNext »Contents
Previous: 4 Backward Monte Carlo Method    Top: 4 Backward Monte Carlo Method    Next: 4.3 Multi-Band Semiconductors