Analysis of the Reformulated Source to Drain Tunneling Probability for Improving the Accuracy of a Multisubband Ensemble Monte Carlo Simulator

Abstract

As an attempt to improve the description of the tunneling current that arises in ultrascaled nanoelectronic devices when charge carriers succeed in traversing the potential barrier between source and drain, an alternative and more accurate non-local formulation of the tunneling probability was suggested. This improvement of the probability computation might result of particular interest in the context of Monte Carlo simulations where the utilization of the conventional Wentzel-Kramers-Brillouin (WKB) approximation tends to overestimate the number of particles experiencing this type of direct tunneling. However, in light of the reformulated expression for the tunneling probability, it becomes of paramount importance to assess the type of potentials for which it behaves adequately. We demonstrate that, for ensuring boundedness, the top of the potential barrier cannot feature a plateau, but rather has to behave quadratically as one approaches its maximum. Moreover, we show that monotonicity of the reformulated tunneling probability is not guaranteed by boundedness and requires an additional constraint regarding the derivative of the prefactor that modifies the traditional WKB tunneling probability.

1. Introduction

The physical phenomenon of quantum mechanical tunneling between source and drain (S/D tunneling) arises when the scaling process of electronic devices enters the regime of channel lengths around and below 10 nm [1,2]. Monte Carlo simulations aiming to account for the description of this type of structures (when this tunneling contribution requires attention in the subthreshold regime) were traditionally providing current levels above those obtained from simulators with full quantum transport treatment [3]. In spite of this drawback, the design and conception of a Multisubband Ensemble Monte Carlo (MS-EMC) tool still provides a very efficient computation technique based on a modular implementation of the specific quantum transport mechanisms involved in the assessment of each device of interest [4,5]. For that reason, the improvement of the MS-EMC tool in order to achieve a satisfactory description of the S/D tunneling phenomena is particularly valuable in the context of the present sizes of ultrascaled devices.

To do so, one of the proposed solutions turned out to be the reformulation of the tunneling probability across the potential barrier between source and drain [6]. This was implemented by means of a non-local approach following the formalism developed in [7], that introduces an additional term correcting the traditional WKB tunneling probability for 2D simulations. Labeling x as the transport direction and z the perpendicular direction affected by confinement, the modified tunneling probability reads as [8]

$$T_{\mathrm{DT}}\left(E_x\right)=\frac{{\mathrm{\Delta }}_y}{2\sqrt{\pi }}{\left[\hslash {\int }_a^b\frac{dx}{\sqrt{2m_x(E_i\left(x\right)-E_x)}}\right]}^{-1/2}·T_{\mathrm{WKB}}\left(E_x\right),$$

(1)

where ${\mathrm{\Delta }}_y$ is the mesh spacing in the periodic direction normal to both transport and confinement dimensions; a and b are the limits of the tunneling path; $m_x$ is the tunneling effective mass of the electron; $E_x$ is the total energy of the superparticle in the transport plane considering only the projection of the kinetic energy in the direction facing the barrier, and $E_i\left(x\right)$ stands for the energy profile of the i-th subband. $T_{\mathrm{WKB}}$ represents the conventional tunneling probability given by the WKB approximation and reads as

$$T_{\mathrm{WKB}}\left(E_x\right)=exp\left[-\frac{2}{\hslash }{\int }_a^b\sqrt{2m_x(E_i\left(x\right)-E_x)}dx\right].$$

(2)

In light of the expression for $T_{\mathrm{DT}}$, its behavior needs to be assessed (due to the term with the negative power) in order to determine for which types of potential barriers it can be used as an appropriate tunneling probability. Should this analysis not be carried out, one might risk to derive untrustworthy conclusions about the S/D importance in those devices where it plays a non-negligible role.

Notice that, in an alternative way to the stochastic conception implied bythe Monte Carlo method, the problem of the lack of validity of the WKB approximation at the turning points of the potential barrier could also be envisaged from a different theoretical perspective employing a wave function formalism. Should this approach be followed, the wave functions would need to be described in in the neighborhoods of the turning points by asymptotic formulas whose leading term would be expressed in terms of Bessel functions.

This paper analyzes the mathematical behavior of the proposed $T_{\mathrm{DT}}$ tunneling probability in the context of the utilizaton of a MS-EMC simulator and delimits the type of potentials for which it proves to be suitable. In Section 2, we expose the behavior of $T_{\mathrm{DT}}$ for a very simple triangular potential profile, illustrating the nature of the problems to be faced. In Section 3, we focus on realistic smooth potentials like those provided by solving the Poisson equation and derive the conditions to be fulfilled by them in order to guarantee an adequate utilization of $T_{\mathrm{DT}}$ when simulating with the MS-EMC tool. Finally, the main conclusions are drawn in Section 4.

2. Reformulated Tunneling Probability Applied to a Simple Potential Barrier

When the carrier impinges on the energy barrier with a certain energy $E_x$, it is obvious that $E_i\left(x\right)$ equals $E_x$ at the point where the carrier enters the barrier (i.e., the point with $x=a$ in the integrals above). The problem arises as long as this leads to a divergence in the integrand of the term with power $-\frac{1}{2}$ in $T_{\mathrm{DT}}$ (which is something not to be concerned about if one simply handles with tunneling probabilities described by $T_{\mathrm{WKB}}$).

In what follows, and for the ease of treatment, all auxiliary constants will be set to one and $T_{\mathrm{DT}}$ will be written as

$$T_{\mathrm{DT}}\left(E_x\right)={\left[G\left(E_x\right)\right]}^{-\frac{1}{2}}{exp}^{-H\left(E_x\right)},$$

(3)

with

$$\begin{array}{ccc}\hfill G\left(E_x\right)={\int }_{a\left(E_x\right)}^{b\left(E_x\right)}g(E_x,x)dx,& \hfill & \hfill H\left(E_x\right)={\int }_{a\left(E_x\right)}^{b\left(E_x\right)}h(E_x,x)dx,\end{array}$$

(4)

and

$$g(E_x,x)=\frac{1}{h(E_x,x)}=\frac{1}{\sqrt{E_i\left(x\right)-E_x}}.$$

(5)

Let us now consider the elementary triangular barrier depicted in Figure 1 given by

$$\begin{array}{c}\hfill E_i\left(x\right)=\left\{\begin{array}{ccc}x& & a\le x\le c\hfill \\ \hfill & & \\ \hfill -(x-b)& & c<x\le b\hfill \end{array}\right..\end{array}$$

(6)

We will suppose that the barrier starts at $x=0$ (i.e., $a=0$) and that the carrier hits the barrier with $E_x=0$.

In that case, $G\left(E_x\right)$ would read as

$$G\left(E_x\right){|}_{E_x=0}={\int }_0^b\frac{dx}{\sqrt{E_i\left(x\right)}}=2{\int }_0^c\frac{dx}{\sqrt{x}}=4{\left[\sqrt{x}\right]}_0^c=4\sqrt{c}.$$

(7)

where the divergence disappears thanks to the integral. However, for this simple example, the resulting behavior of $T_{\mathrm{DT}}$ is incompatible with that corresponding to a well behaved tunneling probability, namely:

$$\begin{array}{c}\hfill \underset{E_x\to max\left(E_i\right)}{lim}G\left(E_x\right)=0\Rightarrow \underset{E_x\to max\left(E_i\right)}{lim}{T}_{\mathrm{DT}}\left(E_x\right)=\infty \end{array}$$

(8)

Therefore, careful assessment needs to be carried out as to determine for which type of potential barriers the utilization of $T_{\mathrm{DT}}$ might be suitable.

3. Reformulated Tunneling Probability Applied to Smooth Potential Barriers

Let us now consider the situation where the potential barrier can be described by a smooth differentiable function between the limits of integration. This situation differs from the simplified scenario analyzed in the previous section (as far as the triangular barrier was not smooth at its maximum), but proves to be more realistic taking into account the usual shape of the potential barriers to be found in scenarios where source to drain tunneling arises.

3.1. Boundedness

According to the form of $g(E_x,x)$, and in order to keep its integral finite when evaluating it at $x=a$ for a given $E_x<max\left(E_i\right)$, we need that

$$\underset{x\to a\left(E_x\right)}{lim}g(Ex,x)<\underset{x\to a\left(E_x\right)}{lim}\frac{1}{x-a\left(E_x\right)}.$$

(9)

Or, in other words, that

$$\underset{x\to a\left(E_x\right)}{lim}\left[x-a\left(E_x\right)\right]g(E_x,x)=\underset{x\to a\left(E_x\right)}{lim}\frac{x-a\left(E_x\right)}{\sqrt{E_i\left(x\right)-E_x}}=0.$$

(10)

If we compute this limit for an energy barrier $E_i\left(x\right)$, we see that both numerator and denominator tend to zero as x approaches $a\left(E_x\right)$. For that reason, the limit can be computed by l’Hôpital as

$$\underset{x\to a\left(E_x\right)}{lim}\frac{x-a\left(E_x\right)}{\sqrt{E_i\left(x\right)-E_x}}=\underset{x\to a\left(E_x\right)}{lim}\frac{1}{\frac{1}{2\sqrt{E_i\left(x\right)-E_x}}{E}_i^{\prime }\left(x\right)}=\underset{x\to a\left(E_x\right)}{lim}\frac{2\sqrt{E_i\left(x\right)-E_x}}{E_i^{\prime }\left(x\right)},$$

(11)

that equals zero provided that the barrier has a certain slope at the impinging point $a\left(E_x\right)$, which is inherent to the nature of the barrier itself. The reasoning is exactly analogous for the exit point $b\left(E_x\right)$. Observe that this result is independent of the shape of $E_i\left(x\right)$, the only requirement is that $E_i^{\prime }\ne 0$ at $a\left(E_x\right)$ and $b\left(E_x\right)$.

Up to this point, the only thing that is guaranteed is that $G\left(E_x\right)$ will remain finite for $E_x<max\left(E_i\right)$ provided that $E_i^{\prime }\ne 0$ at $a\left(E_x\right)$ and $b\left(E_x\right)$. However, we still need $G\left(E_x\right)$ not to go to zero or to infinity as $E_x$ tends to $max\left(E_i\right)$. As we are supposing $E_i\left(x\right)$ to be smooth, it can be expanded in Taylor series around the point $x_0$ where $E_i\left(x\right)$ attains its maximum. This leads to

$$E_i\left(x\right)=c_0+c_1(x-x_0)+c_2{(x-x_0)}^2+c_3{(x-x_0)}^3+c_4{(x-x_0)}^4+\dots ,$$

(12)

with $c_0=max\left(E_i\right)$, $c_1=0$ and $c_2\le 0$. Now suppose that, for simplicity, the potential is centered at the origin so that $x_0=0$. The above expression reduces to

$$E_i\left(x\right)=max\left(E_i\right)+c_2{x}^2+c_3{x}^3+c_4{x}^4+\dots$$

(13)

Since we are interested in the behavior when $E_x$ tends to $max\left(E_i\right)$, this is the same as saying that we are interested in the behavior of the Taylor series when x tends to zero or, in other words, that the key of the analysis relies on the order of the first non vanishing coeficient of the expansion, which is the one controlling how $E_i\left(x\right)$ approaches to its maximum.

Let us suppose that $c_2\ne 0$ (which leads to $c_2<0$ because the potential has a maximum). In that case, we can neglect all the terms of higher order and obtain a simple inverted parabolic potential. If we compute the value of $G\left(E_x\right)$, we obtain

$$\begin{array}{cc}\hfill G\left(E_x\right)& ={\int }_{-\frac{\sqrt{max\left(E_i\right)-E_x}}{\sqrt{-c_2}}}^{\frac{\sqrt{max\left(E_i\right)-E_x}}{\sqrt{-c_2}}}\frac{dx}{\sqrt{max\left(E_i\right)+c_2{x}^2-E_x}}\hfill \\ \hfill & =\frac{1}{\sqrt{-c_2}}{\left[arcsin\left(\frac{\sqrt{-c_2}}{\sqrt{max\left(E_i\right)-E_x}}x\right)\right]}_{-\frac{\sqrt{max\left(E_i\right)-E_x}}{\sqrt{-c_2}}}^{\frac{\sqrt{max\left(E_i\right)-E_x}}{\sqrt{-c_2}}}\hfill \\ \hfill & =\frac{2}{\sqrt{-c_2}}arcsin\left(1\right)=\frac{\pi }{\sqrt{-c_2}}.\hfill \end{array}$$

(14)

This expression is independent of $E_x$ and remains constant when $E_x$ tends to $max\left(E_i\right)$, thus providing a finite non zero value of $G\left(E_x\right)$ even at the top of the barrier.

On the other hand, if $c_2=0$ (which means that the barrier will feature a certain plateau at its top), the next non vanishing (and negative) coeficient will necessarily be $c_j$ with j even, and again all orders above it can be neglected if we want to analyze the behavior when $E_x$ tends to $max\left(E_i\right)$. The corresponding expression for $G\left(E_x\right)$ is then

$$\begin{array}{cc}\hfill G\left(E_x\right)& ={\int }_{a\left(E_x\right)}^{b\left(E_x\right)}\frac{dx}{\sqrt{max\left(E_i\right)+c_j{x}^j-E_x}}(c_j<0,j=4,6,8,\dots ),\hfill \end{array}$$

(15)

with

$$\begin{array}{ccc}\hfill a\left(E_x\right)=-{\left(\frac{max\left(E_i\right)-E_x}{-c_j}\right)}^{\frac{1}{j}},& \hfill & \hfill b\left(E_x\right)={\left(\frac{max\left(E_i\right)-E_x}{-c_j}\right)}^{\frac{1}{j}}.\end{array}$$

(16)

The computation of $G\left(E_x\right)$ by means of elliptic integrals leads to

$$G\left(E_x\right)=\frac{2\sqrt{\pi }\mathrm{\Gamma }\left(\frac{j+1}{j}\right)}{{(-c_j)}^{\frac{1}{j}}{\left[\left(max\left(E_i\right)-E_x\right)\right]}^{\frac{\frac{j}{2}-1}{j}}\mathrm{\Gamma }(\frac{\frac{j}{2}+1}{j})}(c_j<0,j=4,6,8,\dots ),$$

(17)

where $\mathrm{\Gamma }$ stands for the gamma function. The presence of the term $(max\left(E_i\right)-E_x)$ in the denominator makes $G\left(E_x\right)$ diverge as $E_x$ approaches $max\left(E_i\right)$. This result indicates that the potential barrier must necessarily possess a cuadratic behavior close to its maximum, and that the appearance of an upper plateau (for example, if the gate length is not too aggressively scaled) deprives $T_{\mathrm{DT}}$ from its physical meaning. Fortunately, given that S/D tunneling becomes relevant for narrow potentials, this requirement will very likely be satisfied in most of the cases of interest.

3.2. Monotonicity

The fact of assuring that $G\left(E_x\right)$ remains bounded and non zero for any of the impinging energies of interest is not enough for guaranteeing an appropriate behavior of $T_{\mathrm{DT}}$. We need to demand monotonicity too. In other words, the tunneling probability needs to increase as $E_x$ approaches $max\left(E_i\right)$.

Note that, in principle, and from the form of $T_{\mathrm{DT}}$ in Equation (3), one cannot affirm that this will always be the case since $T_{\mathrm{DT}}$ results from the product of a function of $G\left(E_x\right)$ and a function of $H\left(E_x\right)$, and $G\left(E_x\right)$ can very well not be monotonic as shown in Figure 2 for different barrier examples fulfilling boundedness. It is therefore desirable to derive a condition that ensures the monotonicity of $T_{\mathrm{DT}}$ once the shape of the potential is known.

For that aim, we shall demand that $T_{\mathrm{DT}}^{\prime }\left(E_x\right)\ge 0$. In our case, and considering that $T_{\mathrm{DT}}$ always evaluates to a positive number, we can use that

$$T_{\mathrm{DT}}^{\prime }\left(E_x\right)\ge 0\iff {\left[ln{T}_{\mathrm{DT}}\left(E_x\right)\right]}^{\prime }\ge 0,$$

(18)

which proves to be a more convenient condition for studying the behavior of $T_{\mathrm{DT}}$ recalling its expression as a function of $G\left(E_x\right)$ and $H\left(E_x\right)$. Therefore, using Equation (3) we obtain that

$${\left[ln{T}_{\mathrm{DT}}\left(E_x\right)\right]}^{\prime }\ge 0\iff -\frac{1}{2}\frac{G^{\prime }\left(E_x\right)}{G\left(E_x\right)}-H^{\prime }\left(E_x\right)\ge 0\iff \frac{1}{2}\frac{G^{\prime }\left(E_x\right)}{G\left(E_x\right)}+H^{\prime }\left(E_x\right)\le 0.$$

(19)

For calculating the derivatives of $G\left(E_x\right)$ and $H\left(E_x\right)$ we need to apply the Leibniz rule, which leads to

$$\begin{array}{c}\hfill G^{\prime }\left(E_x\right)={\int }_{a\left(E_x\right)}^{b\left(E_x\right)}\frac{\partial g(E_x,x)}{\partial E_x}dx+g(E_x,b\left(E_x\right))b^{\prime }\left(E_x\right)-g(E_x,a\left(E_x\right))a^{\prime }\left(E_x\right),\end{array}$$

(20)

$$\begin{array}{c}\hfill H^{\prime }\left(E_x\right)={\int }_{a\left(E_x\right)}^{b\left(E_x\right)}\frac{\partial h(E_x,x)}{\partial E_x}dx+h(E_x,b\left(E_x\right))b^{\prime }\left(E_x\right)-h(E_x,a\left(E_x\right))a^{\prime }\left(E_x\right).\end{array}$$

(21)

Taking the relationship between $g(E_x,x)$ and $h(E_x,x)$ of Equation (5), we observe that

$$\begin{array}{cc}\hfill \frac{\partial h(E_x,x)}{\partial E_x}& =-\frac{1}{2\sqrt{E_i\left(x\right)-E_x}}=-\frac{1}{2}g(E_x,x),\hfill \end{array}$$

(22)

$$\begin{array}{cc}\hfill \frac{\partial g(E_x,x)}{\partial E_x}& =\frac{\partial }{\partial E_x}\left(\frac{1}{h(E_x,x)}\right)=-\frac{1}{h^2(E_x,x)}\frac{\partial h(E_x,x)}{\partial E_x}=\frac{1}{2}{g}^3(E_x,x).\hfill \end{array}$$

(23)

Which provides

$$\begin{array}{cc}\hfill G^{\prime }\left(E_x\right)& =\frac{1}{2}{\int }_{a\left(E_x\right)}^{b\left(E_x\right)}{g}^3(E_x,x)dx+g(E_x,b\left(E_x\right))b^{\prime }\left(E_x\right)-g(E_x,a\left(E_x\right))a^{\prime }\left(E_x\right),\hfill \end{array}$$

(24)

$$\begin{array}{cc}\hfill H^{\prime }\left(E_x\right)& =-\frac{1}{2}G(E_x,x)+h(E_x,b\left(E_x\right))b^{\prime }\left(E_x\right)-h(E_x,a\left(E_x\right))a^{\prime }\left(E_x\right).\hfill \end{array}$$

(25)

Now, since $h(E_x,x)$ evaluates to zero when $x=a\left(E_x\right)$ and $x=b\left(E_x\right)$, the last equation is simply

$$H^{\prime }\left(E_x\right)=-\frac{1}{2}G\left(E_x\right).$$

(26)

If we apply this result to Equation (19) and use that $G\left(E_x\right)$ is always positive, we obtain that

$$\frac{1}{2}\frac{G^{\prime }\left(E_x\right)}{G\left(E_x\right)}+H^{\prime }\left(E_x\right)\le 0\iff G^{\prime }\left(E_x\right)\le G^2\left(E_x\right)\iff {\left(\frac{1}{G\left(E_x\right)}\right)}^{\prime }\ge -1.$$

(27)

This last inequality constitutes a simple way of formulating the condition that must be satisfied once the shape of the potential barrier $E_i\left(x\right)$ is known in order to guarantee the monotonicity of $T_{\mathrm{DT}}$. In Figure 3, we observe that three of the example barriers of Figure 2 fulfill this condition, so that the probability $T_{\mathrm{DT}}$ would be well behaved for them. On the other hand, for the dashed dotted curve, the probability $T_{\mathrm{DT}}$ would not be monotonic.

We conclude that boundedness does not necessarily imply monotonicity for $T_{\mathrm{DT}}$. For a given potential, both requirements would need to be verified separately.

4. Conclusions

One of the envisageable solutions to improve the accuracy of the S/D tunneling current computation in a MS-EMC simulator is to redefine the tunneling probability by which it is estimated. However, in light of the resulting expression obtained from a non-local reformulation, its applicability is not automatically well defined for any potential barrier that might arise between source and drain regions. We have demonstrated that boundedness of the tunneling probability is ensured for potentials not featuring a flat region around its maximum (i.e., those with non-vanishing quadratic terms when expanded in series). With respect to monotonicity, this property is not automatically ensued from the boundedness of the tunneling probability and requires, in turn, a certain constraint on the boundedness of the derivative of the prefactor modifying the WKB tunneling probability.