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

Citation: Padilla, J.L; Medina-Bailon, C.; Palomares, A.; Donetti, L.; Navarro, C.; Sampedro, C.; Gamiz, F. Analysis of the Reformulated Source to Drain Tunneling Probability for Improving the Accuracy of a Multisubband Ensemble Monte Carlo Simulator. Micromachines 2022,13, 533. https://doi.org/10.3390/ mi13040533 Academic Editor: Niall Tait Received: 23 December 2021 Accepted: 25 March 2022 Published: 28 March 2022 Publisher’s Note: MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations. Copyright: © 2022 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https:// creativecommons.org/licenses/by/ 4.0/). micromachines Article Analysis of the Reformulated Source to Drain Tunneling Probability for Improving the Accuracy of a Multisubband Ensemble Monte Carlo Simulator Jose Luis Padilla 1,2,* , Cristina Medina-Bailon 1,2 , Antonio Palomares 3, Luca Donetti 1,2 , Carlos Navarro 1,2 , Carlos Sampedro 1,2 and Francisco Gamiz 1,2 1Nanoelectronics Research Group, Departamento de Electrónica y Tecnología de Computadores, Universidad de Granada, 18071 Granada, Spain; [email protected] (C.M.-B.); [email protected] (L.D.); [email protected] (C.N.); [email protected] (C.S.); [email protected] (F.G.) 2Centro de Investigación en Tecnologías de la Información y las Comunicaciones, Universidad de Granada, 18071 Granada, Spain 3Departamento de Matemática Aplicada, Universidad de Granada, 18071 Granada, Spain; [email protected] *Correspondence: [email protected] 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-KramersBrillouin (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. Keywords: direct source-to-drain tunneling; tunneling probability; multi-subband ensemble Monte Carlo 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 Micromachines 2022,13, 533. https://doi.org/10.3390/mi13040533 https://www.mdpi.com/journal/micromachines Micromachines 2022,13, 533 2 of 8 for 2D simulations. Labeling x as the transport direction and z the perpendicular direction affected by confinement, the modified tunneling probability reads as [8] TDT(Ex) = ∆y 2√π"¯hZb a dx p2mx(Ei(x)−Ex)#−1/2 ·TWKB(Ex), (1) where ∆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; mx is the tunneling effective mass of the electron; Ex 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 Ei(x) stands for the energy profile of the i -th subband. TWKB represents the conventional tunneling probability given by the WKB approximation and reads as TWKB(Ex) = exp−2 ¯hZb aq2mx(Ei(x)−Ex)dx. (2) In light of the expression for TDT , 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 TDT 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 TDT 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 TDT 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 Ex , it is obvious that Ei(x) equals Ex 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 −1 2 in TDT (which is something not to be concerned about if one simply handles with tunneling probabilities described by TWKB). In what follows, and for the ease of treatment, all auxiliary constants will be set to one and TDT will be written as TDT(Ex) = hG(Ex)i−1 2exp−H(Ex), (3) with G(Ex) = Zb(Ex) a(Ex)g(Ex,x)dx,H(Ex) = Zb(Ex) a(Ex)h(Ex,x)dx, (4) and g(Ex,x) = 1 h(Ex,x)=1 pEi(x)−Ex . (5) Micromachines 2022,13, 533 3 of 8 Let us now consider the elementary triangular barrier depicted in Figure 1given by Ei(x) =    x a ≤x≤c −(x−b)c<x≤b . (6) We will suppose that the barrier starts at x= 0 (i.e., a= 0) and that the carrier hits the barrier with Ex=0. Figure 1. Triangular barrier profile where the carrier impinges on the barrier with Ex= 0. If the carrier suceeds in traversing the barrier, it exits at the point x=bwith the same energy. In that case, G(Ex)would read as G(Ex)Ex=0 =Zb 0 dx pEi(x)=2Zc 0 dx √x=4√xc 0=4√c. (7) where the divergence disappears thanks to the integral. However, for this simple example, the resulting behavior of TDT is incompatible with that corresponding to a well behaved tunneling probability, namely: lim Ex→max(Ei)G(Ex) = 0⇒lim Ex→max(Ei)TDT(Ex) = ∞(8) Therefore, careful assessment needs to be carried out as to determine for which type of potential barriers the utilization of TDT 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(Ex , x) , and in order to keep its integral finite when evaluating it at x=afor a given Ex<max(Ei), we need that lim x→a(Ex)g(Ex,x)<lim x→a(Ex) 1 x−a(Ex). (9) Or, in other words, that lim x→a(Ex)x−a(Ex)g(Ex,x) = lim x→a(Ex) x−a(Ex) pEi(x)−Ex =0. (10) Micromachines 2022,13, 533 4 of 8 If we compute this limit for an energy barrier Ei(x) , we see that both numerator and denominator tend to zero as x approaches a(Ex) . For that reason, the limit can be computed by l’Hôpital as lim x→a(Ex) x−a(Ex) pEi(x)−Ex =lim x→a(Ex) 1 1 2√Ei(x)−Ex E0 i(x)=lim x→a(Ex) 2pEi(x)−Ex E0 i(x), (11) that equals zero provided that the barrier has a certain slope at the impinging point a(Ex) , which is inherent to the nature of the barrier itself. The reasoning is exactly analogous for the exit point b(Ex) . Observe that this result is independent of the shape of Ei(x) , the only requirement is that E0 i6=0 at a(Ex)and b(Ex). Up to this point, the only thing that is guaranteed is that G(Ex) will remain finite for Ex<max(Ei) provided that E0 i6= 0 at a(Ex) and b(Ex) . However, we still need G(Ex) not to go to zero or to infinity as Ex tends to max(Ei) . As we are supposing Ei(x) to be smooth, it can be expanded in Taylor series around the point x0 where Ei(x) attains its maximum. This leads to Ei(x) = c0+c1(x−x0) + c2(x−x0)2+c3(x−x0)3+c4(x−x0)4+. . . , (12) with c0=max(Ei) , c1= 0 and c2≤ 0. Now suppose that, for simplicity, the potential is centered at the origin so that x0=0. The above expression reduces to Ei(x) = max(Ei) + c2x2+c3x3+c4x4+. . . (13) Since we are interested in the behavior when Ex tends to max(Ei) , 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 Ei(x) approaches to its maximum. Let us suppose that c26= 0 (which leads to c2< 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(Ex), we obtain G(Ex) = Z√max(Ei)−Ex √−c2 −√max(Ei)−Ex √−c2 dx pmax(Ei) + c2x2−Ex =1 √−c2"arcsin √−c2 pmax(Ei)−Ex x# √max(Ei)−Ex √−c2 −√max(Ei)−Ex √−c2 (14) =2 √−c2 arcsin(1) = π √−c2 . This expression is independent of Ex and remains constant when Ex tends to max(Ei) , thus providing a finite non zero value of G(Ex)even at the top of the barrier. On the other hand, if c2= 0 (which means that the barrier will feature a certain plateau at its top), the next non vanishing (and negative) coeficient will necessarily be cj with j even, and again all orders above it can be neglected if we want to analyze the behavior when Extends to max(Ei). The corresponding expression for G(Ex)is then G(Ex) = Zb(Ex) a(Ex) dx qmax(Ei) + cjxj−Ex (cj<0, j=4, 6, 8, . . .), (15) Micromachines 2022,13, 533 5 of 8 with a(Ex) = −max(Ei)−Ex −cj1 j ,b(Ex) = max(Ei)−Ex −cj1 j . (16) The computation of G(Ex)by means of elliptic integrals leads to G(Ex) = 2√πΓj+1 j (−cj)1 jhmax(Ei)−Exi j 2−1 jΓj 2+1 j (cj<0, j=4, 6, 8, . . .), (17) where Γ stands for the gamma function. The presence of the term (max(Ei)−Ex) in the denominator makes G(Ex) diverge as Ex approaches max(Ei) . 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 TDT 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(Ex) remains bounded and non zero for any of the impinging energies of interest is not enough for guaranteeing an appropriate behavior of TDT . We need to demand monotonicity too. In other words, the tunneling probability needs to increase as Exapproaches max(Ei). Note that, in principle, and from the form of TDT in Equation (3), one cannot affirm that this will always be the case since TDT results from the product of a function of G(Ex) and a function of H(Ex) , and G(Ex) 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 TDT once the shape of the potential is known. -1.0 -0.8 -0.6 -0.4 -0.2 2.0 2.5 3.0 3.5 Figure 2. Behavior of the factor G(Ex) for different potential barriers Ei(x) with max(Ei) = 0. Notice how, as Ex→max(Ei),G(Ex)tends to the finite expected value of π √−c2=π. For that aim, we shall demand that T0 DT(Ex)≥ 0. In our case, and considering that TDT always evaluates to a positive number, we can use that T0 DT(Ex)≥0⇔ln TDT(Ex)0≥0, (18) Micromachines 2022,13, 533 6 of 8 which proves to be a more convenient condition for studying the behavior of TDT recalling its expression as a function of G(Ex) and H(Ex) . Therefore, using Equation (3) we obtain that ln TDT(Ex)0≥0⇔ −1 2 G0(Ex) G(Ex)−H0(Ex)≥0⇔1 2 G0(Ex) G(Ex)+H0(Ex)≤0. (19) For calculating the derivatives of G(Ex) and H(Ex) we need to apply the Leibniz rule, which leads to G0(Ex) = Zb(Ex) a(Ex) ∂g(Ex,x) ∂Ex dx +g(Ex,b(Ex))b0(Ex)−g(Ex,a(Ex))a0(Ex), (20) H0(Ex) = Zb(Ex) a(Ex) ∂h(Ex,x) ∂Ex dx +h(Ex,b(Ex))b0(Ex)−h(Ex,a(Ex))a0(Ex). (21) Taking the relationship between g(Ex , x) and h(Ex , x) of Equation (5), we observe that ∂h(Ex,x) ∂Ex =−1 2pEi(x)−Ex =−1 2g(Ex,x), (22) ∂g(Ex,x) ∂Ex =∂ ∂Ex1 h(Ex,x)=−1 h2(Ex,x) ∂h(Ex,x) ∂Ex =1 2g3(Ex,x). (23) Which provides G0(Ex) = 1 2Zb(Ex) a(Ex)g3(Ex,x)dx +g(Ex,b(Ex))b0(Ex)−g(Ex,a(Ex))a0(Ex), (24) H0(Ex) = −1 2G(Ex,x) + h(Ex,b(Ex))b0(Ex)−h(Ex,a(Ex))a0(Ex). (25) Now, since h(Ex , x) evaluates to zero when x=a(Ex) and x=b(Ex) , the last equation is simply H0(Ex) = −1 2G(Ex). (26) If we apply this result to Equation (19) and use that G(Ex) is always positive, we obtain that 1 2 G0(Ex) G(Ex)+H0(Ex)≤0⇔G0(Ex)≤G2(Ex)⇔1 G(Ex)0≥ −1. (27) This last inequality constitutes a simple way of formulating the condition that must be satisfied once the shape of the potential barrier Ei(x) is known in order to guarantee the monotonicity of TDT . In Figure 3, we observe that three of the example barriers of Figure 2 fulfill this condition, so that the probability TDT would be well behaved for them. On the other hand, for the dashed dotted curve, the probability TDT would not be monotonic. We conclude that boundedness does not necessarily imply monotonicity for TDT . For a given potential, both requirements would need to be verified separately. Micromachines 2022,13, 533 7 of 8 -1.0 -0.8 -0.6 -0.4 -0.2 -1.0 -0.5 0.5 Figure 3. Behavior of the derivative of 1 /G(Ex) for the potential barriers considered in Figure 2. Three of them fulfill the condition that ensures monotonicity for TDT, whereas one does not. 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 nonlocal 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. Author Contributions: Conceptualization, J.L.P. and C.M.-B.; methodology, J.L.P. and A.P.; validation, C.M.-B. and L.D.; formal analysis, J.L.P. and A.P.; investigation, J.L.P. and C.M.-B.; writing—original draft preparation, J.L.P.; writing—review and editing, L.D., A.P. and C.S.; visualization, C.M.-B.; supervision, F.G.; project administration, J.L.P.; funding acquisition, C.N. and L.D. All authors have read and agreed to the published version of the manuscript. Funding: This research was funded by Junta de Andalucia, grant number P18-RT-4826. Institutional Review Board Statement: Not applicable. Informed Consent Statement: Not applicable. Data Availability Statement: Not applicable. Conflicts of Interest: The authors declare no conflict of interest. The funders had no role in the design of the study; in the collection, analyses, or interpretation of data; in the writing of the manuscript, or in the decision to publish the results. Abbreviations The following abbreviations are used in this manuscript: S/D Source to drain MS-EMC Multisubband Ensemble Monte Carlo WKB Wentzel-Kramers-Brillouin References 1. Iwai, H. Future of nano CMOS technology. Solid-State Electron. 2015,112, 56–67. [CrossRef] 2. Grillet, C.; Logoteta, D.; Cresti, A.; Pala, M.G. Assessment of the Electrical Performance of Short Channel InAs and Strained Si Nanowire FETs. IEEE Trans. Electron. Devices 2017,64, 2425–2431. [CrossRef] Micromachines 2022,13, 533 8 of 8 3. Medina-Bailon, C.; Sampedro, C.; Padilla, J.L.; Godoy, A.; Donetti, L.; Gamiz, F.; Asenov, A. MS-EMC vs. NEGF: A comparative study accounting for transport quantum corrections. In Proceedings of the Joint International EUROSOI Workshop and International Conference on Ultimate Integration on Silicon (EUROSOI-ULIS 2018), Granada, Spain, 19–21 March 2018; pp. 1–4. 4. Sampedro, C.; Medina-Bailon, C.; Donetti, L.; Padilla, J.L.; Navarro, C.; Marquez, C.; Gamiz, F. Multi-Subband Ensemble Monte Carlo Simulator for Nanodevices in the End of the Roadmap. In Large-Scale Scientific Computing, Proceedings of the 12th International Conference, LSSC 2019, Sozopol, Bulgaria, 10–14 June 2019; Lecture Notes in Computer Science; Lirkov, I.; Margenov, S., Eds.; Springer: New York, NY, USA, 2019; Volume 11958, pp. 438–445. 5. Medina-Bailon, C.; Padilla, J.L.; Sadi, T.; Sampedro, C.; Godoy, A.; Donetti, L.; Georgiev, V.; Gamiz, F.; Asenov, A. Multisubband Ensemble Monte Carlo Analysis of Tunneling Leakage Mechanisms in Ultrascaled FDSOI, DGSOI and FinFET Devices. IEEE Trans. Electron. Devices 2019,66, 1145–1152. [CrossRef] 6. Medina-Bailon, C.; Carrillo-Nunez, H.; Lee, J.; Sampedro, C.; Padilla, J.L.; Donetti, L.; Georgiev, V.; Gamiz, F.; Asenov, A. Quantum Enhancement of a S/D Tunneling Model in a 2D MS-EMC Nanodevice Simulator: NEGF Comparison and Impact of Effective Mass Variation. Micromachines 2020,11, 204. [CrossRef] [PubMed] 7. Nunez, H.C.; Ziegler, A.; Luisier, M.; Schenk, A. Modeling direct band-to-band tunneling: From bulk to quantum-confined semiconductor devices. J. Appl. Phys. 2015,117, 243501. 8. Medina-Bailon, C.; Padilla, J.L.; Sampedro, C.; Donetti, L.; Georgiev, V.; Gamiz, F.; Asenov, A. Self-Consistent Enhanced S/D Tunneling Implementation in a 2D MS-EMC Nanodevice Simulator. Micromachines 2021,12, 601. [CrossRef] [PubMed]