« Previous Up Next »Contents
Previous: 3.3 Monte Carlo Method Top: Home Next: 4.2 Current Estimators

Investigating Hot-Carrier Effects using the Backward Monte Carlo Method

4 Backward Monte Carlo Method

The backward MC method (BMC) was introduced at the end of the 1980s [45, 73]. These early algorithms turned out to be numerically unstable, as the carrier energy tends to grow indefinitely on a trajectory that is followed back in time [P4]. A numerically stable algorithm was proposed in 2003 [62]. Since the backward transition rates are chosen to obey the principle of detailed balance, a runaway of the carrier energy along a backward trajectory is avoided. From a practical point of view, this means that the scattering rates of the forward method can be used in the backward method as well [62].

The principle of the BMC method for the solution of a boundary value problem is to choose a set of states in phase space and trace trajectories from these states back in time until a contact is reached. The value of the given distribution function (DF) at the contact determines the statistical weight of the backward trajectory and consequently its contribution to the estimator of interest [P4].

This method enables the calculation of a current, that is controlled by an energy barrier. The current through a device is typically determined by the states at the top of the barrier, see Fig. 4.1. If the barrier is high, a forward trajectory is very unlikely to reach the top of the barrier, whereas in the backward method only these unlikely states are considered, and no computation time is wasted with the vast majority of trajectories that do not overcome the barrier [P4].

It is also possible to combine the backward and the forward MC method. Once a backward trajectory with an initial state $(\vec {v}_0, \vec {r}_0)$ is calculated, and the statistical weight of that state is determined, a forward trajectory can be started from the very same state (Fig. 4.5). The mean values of interest are then calculated from a set of forward trajectories in the usual manner [P3, P4].

4.1 Theory of the Backward Monte Carlo Method

With the initial conditions

$(4.1) \{begin}{align} \vec {K}_0(t_0) = \vec {k}_0\qquad \text {and}\qquad \vec {R}_0(t_0) = \vec {r}_0, \{end}{align}$

a phase space trajectory can be obtained by formally integrating the equations of motion (2.22) and (2.23) [P4]:

$(4.2–4.3) \{begin}{align} \vec {K}_0(t) &= \vec {k}_0 + \int _{t_0}^{t} \vec {F}_e\left ( \vec {R}_0(\tau ),\tau \right ) \dint \tau \„\\ \vec {R}_0(t) &= \vec {r}_0 + \int _{t_0}^{t} \vec {v}\left ( \vec {K}_0(\tau )\right ) \dint \tau . \{end}{align}$

The BTE (2.26) can be integrated over the phase space trajectory in the manner of [59]:

$(4.4) \{begin}{align} f(\kv _0,\rv _0,t_0) = \int _0^{t_0} \!\d t_1\!\int \! \dk _1\, K(\kv _0,t_0,\kv _1, t_1)\,f(\vect {k}_1,\vect {R}(t_1),t_1) + f_0(\kv _0,\rv _0,t_0) \label {eq:int-form} \{end}{align}$

The resulting integral equation (4.4) represents the generalization of Chamber’s path integral [13, P4, 64]. The source term of this equation includes the initial distribution for a initial value problem [61], or the boundary distribution for a boundary value problem [59]. The kernel of the integral equation is of the form:

$(4.5) \{begin}{align} K(\kv _0,t_0,\kv _1,t_1) = S(\kv _1,\Kv _0(t_1)) \exp \biggl (-\!\!\int _{t_1}^{t_0}\!\lambda (\Kv _0(\tau ))\d \tau \biggr ) \label {eq:kernel} \{end}{align}$

The trajectory $\vec {K}_0(\tau )$ passes through $\vec {k}_0$ at the time $t_0$ . The kernel (4.5) in a physical sense describes a transition from $(\vec {k}_1,t_1)$ to $(\vec {k}_0,t_0)$ [P4].

4.1.1 Probability Density Functions

The components of the kernel (4.5) allow the construction of probability density functions (PDF). From the scattering rate $S$ a PDF of the after-scattering states $\vec {k}_a$ can be constructed [P4]:

$(4.6) \{begin}{align} p_k(\kv _a| \kv _b) = \frac {S(\kv _b, \kv _a)}{\lambda (\kv _b)} \label {eq:pk-forw} \{end}{align}$

The PDF of the before-scattering states $\vec {k}_b$ can be constructed in a reverse manner:

$(4.7) \{begin}{align} p_k^*(\kv _b| \kv _a) = \frac {S(\kv _b, \kv _a)}{\lambda ^*(\kv _a)}. \label {pk-back} \{end}{align}$

The scattering rates

$(4.8–4.9) \{begin}{align} \lambda (\vect {k}_b) &= \int S(\kv _b,\kv _a)\,\dk _a \„\\ \lambda ^*(\vect {k}_a) &= \int S(\kv _b,\kv _a)\,\dk _b \„ \{end}{align}$

are serving as normalization factors. The path integral in (4.5) is leading to the PDF of the backward free-flight time $t_1$ [P4]:

$(4.10) \{begin}{align} p_t^*(t_1|t_0; \kv _0) = {\lambda (\Kv _0(t_1))} \exp \biggl (-\!\int _{t_1}^{t_0}\!\lambda (\Kv _0(\tau ))\,\d \tau \biggr )\„ \qquad t_1 < t_0 \label {eq:pb} \{end}{align}$

With the transformation $t^* = -t$ , $\kv ^*=-\kv$ , and $\vv ^* = -\vv$ , the equations of motion (2.22) and (2.23) can be shown to be form-invariant [P4]. Therefore, the equations of motion for the forward path can also be used for the backward path. The vector $\vec {r}$ and the force field $\vec {F}_e$ need not be inverted. Consequently, the substitution $\tau ^* = -\tau$ transforms the PDF of the backward free-flight time (4.10) to the PDF of the forward free-flight [P4]:

$(4.11) \{begin}{align} p_t(t_1^*|t_0^*; \kv _0^*) = {\lambda (\Kv _0^*(t_1^*))} \exp \biggl (-\!\int _{t_0^*}^{t_1^*}\!\lambda (\Kv _0^*(\tau ^*))\, \d \tau ^*\biggr ), \qquad t_1^* > t_0^* \label {eq:pf} \{end}{align}$

Both PDFs (4.10) and (4.11) are normalized:

$(4.12–4.13) \{begin}{align} \int _{-\infty }^{t_0} p_t^*(t_1|t_0)\,\d t_1 &= 1 \„\\ \int ^{\infty }_{t_0^*} p_t(t_1^*|t_0^*)\,\d t_1^* &= 1 \,. \{end}{align}$

4.1.2 The Backward MC Method

In the more familiar forward MC method, the scattering events occur in a ascending time sequence: $t_0 < t_1 <t_2 < \cdots$ . In the the backward MC method based on (4.4) the scattering events occur in a descending sequence : $t_0 > t_1 >t_2 >\cdots$ .

The distribution function $f$ in one point $(\vec {k}_0,\vec {r}_0)$ can be estimated by the BMC method by the following sample mean [P4]:

$(4.14) \{begin}{align} f(\kv _0,\rv _0,t_0) \simeq \frac {1}{N}\sum \limits _{s=1}^N \mu _s^{(n(s))}(\kv _0,\rv _0,t_0) \label {eq:sample-mean} \{end}{align}$

The number of trajectories is represented by $N$ , and $n(s)$ is denoting the order of the $s$ -th numerical trajectory, which is the number of scattering events occurring in the time interval $[0,t_0]$ . In literature, there are two different kinds of estimators $\mu ^{(n)}$ . One is based on mathematical considerations, the other on physical considerations. Both of them will be discussed below.

Transition Rate Derived from Mathematical Considerations

The first works regarding the BMC method, [45] and [73], interpreted $S(\kv _b, \kv _a)$ as the unnormalized distribution of the before-scattering states $\kv _b$ . Thus, the normalized PDF (4.10) is applied [P4]. With the transition density

$(4.15) \{begin}{align} P(\kv _1,t_1|\kv _0,t_0) = p_k^*(\kv _1|\Kv _0(t_1))\; p_t(t_1|t_0;\kv ) \{end}{align}$

the estimator in (4.14) becomes

$(4.16) \{begin}{align} \mu ^{(n)}(\kv _0,\rv _0,t_0) = \frac {\lambda ^*(\Kv _0(t_1))}{\lambda (\Kv _0(t_1))} \ldots \frac {\lambda ^*(\Kv _{n-1}(t_n))}{\lambda (\Kv _{n-1}(t_n))} f_\text {in}(\Kv _n(0),\Rv _n(0))\„ \label {estimator-mixi} \{end}{align}$

where $f_\text {in}$ denotes the initial distribution. A trajectory of second order ( $n=2$ ) is illustrated in Fig. 4.2.

Though the estimator (4.16) is formally derived from the BTE, previous simulations revealed a stability problem. The energy of one particle becomes statistical very high when the trajectory is followed back in time, as sketched in Fig. 4.3. The initial distribution takes on very small values for high energies. Whereas the probability that the particle energy becomes low is very small, but the initial distribution for low energies is high. These rare events are contributing mainly to the estimator and causing a large variance. The simulations show that the variance is increasing quickly over time. However, for a finite time $t$ the variance of the estimator is finite [P4].

Transition Rate Derived from Physical Considerations, based on [P4]

The time evolution of the particle energy can be understood from a property of the scattering rate known as the principle of detailed balance. This property ensures that in any system particles scatter preferably to lower energies. If the backward transition rate (4.7) is employed for trajectory construction, in the simulation the principle of detailed balance is inverted, and scattering to higher energies is preferred.

The principle of detailed balance is reflected by the following symmetry property of the scattering rate [61]:

$(4.17) \begin{equation} S(\kv _i,\kv _j) = S(\kv _j,\kv _i)\, \e ^{\,\beta _D(\Epsilon (\kv _i) - \Epsilon (\kv _j))}\„ \label {detailed} \end{equation}$

where $\beta _D=(k_B T_D)^{-1}$ with $T_D$ being the device temperature, and $\Epsilon (\kv )$ denoting the carrier energy. The stability problem can be solved by using the forward scattering rate also for the construction of the backward trajectory and changing the estimator accordingly, as sketched in Fig. 4.4. In the transition density the forward PDF (4.6) is employed [61].

$(4.18) \begin{equation} P(\kv _1,t_1 | \kv _0, t_0) = p_k(\kv _1|\Kv _0(t_1))\;p_t(t_1|t_0;\kv _0) \end{equation}$

The estimator in (4.14) becomes

$(4.19) \begin{equation} \mu ^{(n)}(\kv ,\rv ,t) = \e ^{\,\beta _D\Delta \Epsilon _1}\ldots \e ^{\,\beta _D\Delta \Epsilon _n} f_\text {in}(\Kv _n(0),\Rv _n(0))\,. \label {estimator-stable} \end{equation}$

Here, the difference in carrier energy induced by the $l$ -th scattering event is denoted by $\Delta \Epsilon _l$ .

The result obtained for the initial value problem can be adapted straightforwardly for the stationary boundary value problem. The total electron energy is defined as

$(4.20) \begin{equation} H(\kv ,\rv ) = \Epsilon (\kv ) + E_C(\rv )\, \label {eq:total-energy} \end{equation}$

where $E_C(\rv )$ is the conduction band edge. In the following, the starting point of a backward trajectory is labeled as $(\kv _0$ , $\rv _0)$ , and the endpoint at a Dirichlet boundary as $(\kv _b$ , $\rv _b)$ , see Fig. 4.1. A Dirichlet boundary is imposed by an ohmic contact, where an equilibrium distribution can be assumed. The boundary distribution at a contact will be referred to as $f_b$ . Inelastic scattering events cause a difference in the total energy along a trajectory. Using the energy balance equation

$(4.21) \begin{equation} H(\kv _b, \rv _b) - H(\kv _0, \rv _0) = \sum \limits _l\Delta \Epsilon _l \end{equation}$

the estimator (4.19) becomes

$(4.22) \{begin}{align} \mu ^{(n)}(\kv ,\rv ,t) = \e ^{\,\beta _D(H(\kv _b, \rv _b) - H(\kv _0, \rv _0))} f_\text {b}(\kv _{b,i},\rv _{b,i})\„ \{end}{align}$

where the initial distribution $f_{in}$ has been replaced by the boundary distribution $f_b$ . Finally, the distribution function in a given point $(\kv _0,\rv _0)$ is estimated by the sample mean (4.14). The distribution function in a given point $(\kv _0, \rv _0)$ becomes

$(4.23) \begin{equation} f(\kv _0, \rv _0) = \frac {1}{M} \sum \limits _{i=1}^M f_b(\kv _{b,i}, \rv _{b.i})\, \e ^{\,\beta _D(H(\kv _{b,i},\rv _{b,i}) - H(\kv _0,\rv _0))}\,. \label {eq:bmcest} \end{equation}$

Here, $M$ is the number of backward trajectories started from the point $(\kv _0, \rv _0)$ . Note that the backward trajectory is constructed in the very same manner as a forward trajectory. Using the forward PDF (4.11) to generate the free-flight time means that we have inverted the time axis and are progressing along the negative time axis. The selection of the scattering mechanism and the calculation of the after-scattering state are also identical to the forward algorithm.

« Previous Up Next »Contents
Previous: 3.3 Monte Carlo Method Top: Home Next: 4.2 Current Estimators