Quantum Enhancement of a S/D Tunneling Model in a 2D MS-EMC Nanodevice Simulator: NEGF Comparison and Impact of Effective Mass Variation

As complementary metal-oxide-semiconductor (CMOS) transistors approach the nanometer scale, it has become mandatory to incorporate suitable quantum formalism into electron transport simulators. In this work, we present the quantum enhancement of a 2D Multi-Subband Ensemble Monte Carlo (MS-EMC) simulator, which includes a novel module for the direct Source-to-Drain tunneling (S/D tunneling), and its verification in the simulation of Double-Gate Silicon-On-Insulator (DGSOI) transistors and FinFETs. Compared to ballistic Non-Equilibrium Green’s Function (NEGF) simulations, our results show accurate ID vs. VGS and subthreshold characteristics for both devices. Besides, we investigate the impact of the effective masses extracted Density Functional Theory (DFT) simulations, showing that they are the key of not only the general thermionic emission behavior of simulated devices, but also the electron probability of experiencing tunneling phenomena.


Introduction
Conventional and novel transistor architectures have been scaled down in the last decades to achieve better performance and larger integration with both lower power consumption and cost. However, conventional bulk complementary metal-oxide-semiconductor (CMOS) technologies present different problems with scaling, such as short-channel effects, reduction of the mobility, leakage current, degradation of the ON and OFF currents (I ON /I OFF ), or variability issues. Novel CMOS transistor architectures, such as Full-Depleted Silicon-On-Insulator (FDSOI) and FinFET, were introduced to mitigate the unwanted effects. Furthermore, in the area of nanotransistor transport simulations, one needs to assess the importance of new phenomena that were not relevant in previous technological nodes [1] in order to explain the electrical behavior of aggressively scaled nanodevices. The simulation of these new phenomena is therefore mandatory to investigate and design the next technology generations and to extend the end of the scaling Roadmap.
At present, different approaches incorporating quantum confinement and tunneling into semi-classical models have become popular due to their modular implementation and reduced computational cost in comparison to purely quantum transport simulation techniques. In particular, the direct Source-to-Drain tunneling (S/D tunneling) starts to play an important role degrading the subthreshold behavior when the channel length is reduced to below 10 nm [2,3], being traditionally considered as a scaling limit in ballistic Non-Equilibrium Green's Function (NEGF) calculations [4], distorting the MOSFET operation at transistor channel lengths around 3nm [2]. This tunneling mechanism allows electrons to tunnel from the source to the drain through the narrow potential barrier existing between both regions, which is controlled by the gate. As a result, the current is increased, eroding the gate control and the subthreshold slope and increasing the leakage. In the simulation of the direct tunneling phenomena, the employed band structure model must accurately represent the experimental energy gaps and effective masses for the most relevant subbands.
The aim of this work is twofold. First, we will discuss the quantum upgrade of our semi-classical 2D Multi-Subband Ensemble Monte Carlo (MS-EMC) simulation tool [5] through the inclusion of a novel S/D tunneling model. Second, we will perform a comprehensive study of how the effective mass variation in confined channels impacts the transport properties and the S/D tunneling. In particular, we have calibrated our tunneling model against the 2D NEGF solver included in the new simulation environment Nano-Electronic Simulation Software (NESS) [6] bearing in mind that the S/D tunneling is naturally included in the quantum NEGF approach. To understand the impact of electron effective mass variation, the bulk effective masses (m bulk ) and calibrated effective masses (m e f f ) from Density Functional Theory (DFT) are used in the MS-EMC model, which is the preferred technique to calculate the electronic band structure of confined nanostructures.
The paper is organized as follows. Section 2 provides a general overview of the simulation framework, including the outline of the simulated devices (Section 2.1), a brief description of NEGF-NESS (Section 2.2) and MS-EMC tools (Section 2.3). It also reports the S/D tunneling model incorporated in the MS-EMC tool (Section 2.4), together with the effective mass calculation and the corresponding extracted values (Section 2.5). The main findings are presented in Section 3 considering ballistic transport for MS-EMC vs. NEGF comparison (Section 3.1) as well as diffusive simulations for the study of the effective mass variation impact (Section 3.2). Finally, conclusions are drawn in Section 4.

Description of the Simulated Devices
In this work, we have compared two double-gate device architectures, with the main difference related to the confinement direction: a planar Double-Gate Silicon-On-Insulator (DGSOI) transistor and a vertical FinFET. Their description including orientation and design parameters can be found in Figure 1. The corresponding bulk effective masses are summarized in Table 1. The confinement direction for these devices on standard wafers [100] changes from (100) for DGSOI to (011) for FinFET, whereas the transport direction <011> is the same for both. The difference in the confinement direction modifies the electron distribution in the subbands and, consequently, the electrostatic potential profile. In addition, the carrier transport effective mass is also modified [7] as it is shown in Section 2.5.  Table 1. Silicon bulk effective masses (m bulk ) for the different crystallographic directions considered in the DGSOI and FinFET devices. Herein, m x and m z are the transport and confinement masses, respectively; m y is the effective mass in the periodic transverse direction; m 0 is the free electron mass; and the subindex of ∆ represents the degeneracy factor associated with the conduction band valley. At this stage, it is important to highlight that, although the FinFET is a 3D structure and our simulation approach is 2D, it has been shown that FinFETs with fin heights much higher than the corresponding thicknesses show similar behavior in all transport regimes when using 2D and 3D simulations [8].

Device
The devices under consideration have been parametrized for gate lengths ranging from 5 to 15 nm. A channel thickness T Si = 3 nm has been chosen for the MS-EMC vs. NEGF comparison. As for the effective mass variation impact, these devices have been simulated for two different channel thickness T Si = 5 nm and T Si = 3 nm in order to capture the effect of the channel thickness reduction. The rest of the technological parameters remains constant: a SiO 2 gate oxide with EOT = 1 nm and a metal gate work function of 4.385 eV.

Description of the 2D NEGF Module Inside NESS
The effective-mass real-space Hamiltonian can be expressed as, assuming y-direction as the periodic transverse direction. The total energy E can also be written where E x is the electron energy in the transport direction. The Hamiltonian in Equation (1) is then transformed to the mode-space representation in order to reduce the computational cost of quantum transport simulations [9]. This was carried out by means of a recursive NEGF approach [10] as implemented in NESS [6] to extract the most relevant physical quantities such as the carrier charge and current. Further, we briefly summarize the main expressions of the NEGF formalism.
For 2D devices and by exploiting the effective-mass approximation, all the transverse modes k y can be treated as independent devices in parallel. Then, within the ballistic regime and under steady-state conditions, the retarded and lesser Green's function for the active device region are written, in matrix notation, as: where H M and Σ R/< C are the Hamiltonian and the retarded/lesser contact self-energies (C = S/D) in the mode-space representation, respectively. The retarded Green's function at the contacts G R = g R C is calculated by means of the Sancho-Lopez-Rubio recursive method [11], allowing straightforwardly the where the mode-space hopping parameters t M are computed as in Ref. [9]. The lesser contact self-energy Σ < C can be then computed from: with where f S/D is the Fermi-Dirac distribution and L y is the periodic length in y-direction. Finally, the 3D carrier concentration and current are calculated in the mode-space representation as follows: where φ n (z) is the confinement wave-function for the subband n, whereas, matrices t M (i) couple two successive layers. Finally, Equations (1) to (3), (5) and (7) are solved self-consistently with Poisson's equation.

General Overview of the 2D MS-EMC Tool
The 2D MS-EMC simulation framework [5] used in this work is based on a decoupled mode-space quantum transport [12] and a semi-classical approach. The simulator solves the Schrödinger equation in the discretization slices along the confinement direction and the Boltzmann Transport Equation (BTE) in the transport plane ( Figure 1). Both equations are coupled through the 2D Poisson equation in the whole 2D simulation domain to keep the self-consistency of the solution. This tool has been widely used in different scenarios including the study of different tunneling mechanisms in similar devices [13]. Due to the modular design of our MS-EMC tool, the inclusion of these tunneling phenomena can be successfully included via additional modules that treat them as separate transport mechanisms without increasing the computational time in comparison to purely quantum simulators. These modules can be switched on or off depending on the simulation scenario, offering the possibility of studying each tunneling mechanism independently.

Description of the S/D Tunneling Model Inside the 2D MC-EMC Tool
S/D tunneling has been included as a separated transport mechanism in the 2D MS-EMC tool described in Section 2.3. It has been implemented as a stochastic mechanism evaluated for each superparticle at the end of the Monte Carlo cycle [14]. When this tunneling mechanism is considered, an electron near the S/D potential barrier will be either reflected or transmitted through it.
The first step is to calculate the tunneling probability by the Wentzel-Kramers-Brillouin (WKB) approximation [15]. It mainly depends on the energy and position of the carrier in the device; the transport effective mass (namely m x in Tables 1 and 2); and the energy profile of the i-th subband determining the shape of the tunneling barrier (E i (x)), which is calculated as a solution of the 1D Schrödinger equation. The probability of tunneling through the barrier is equivalent to the transmission coefficient, and determines the fraction of electrons experiencing S/D tunneling at a given energy below the potential barrier. The tunneling probability of the electron for a given energy (T WKB (E x )) is: where a and b are the limits of the tunneling path, and E x is the total energy in the transport plane considering only the projection of the kinetic energy in the direction that faces the potential barrier. It has been reported for the short-gate length devices that this model overesitimated the number of superparticles experiencing S/D tunneling compared to NEGF approach [16]. This model was compared to NEGF simulations showing an overestimation of the number of superparticles experiencing S/D tunneling, especially for short-gate length devices. In order to reduce this discrepancy, the tunneling model in Equation (9) has been reformulated following a non-local WKB probability approach as stated in Appendix B of Ref. [17]. In the context of a 2D simulation domain, the new S/D tunneling probability for a given energy (T DT (E x )) is now defined as: where ∆ y is the mesh spacing in the direction normal to transport. As this direction is not taken into account in our 2D MS-EMC tool, the value of ∆ y has been calibrated to fulfill the following conditions: (i) the force T DT (E x ) to be in the range [0-1]; (ii) to be small enough to be consistent with the periodic boundary condition in the y direction; and (iii) to have similar degradation in the subthreshold region compared to NEGF calculations for the device with L G = 7.5 nm (shown in Section 3.1). In order to assess the S/D tunneling impact as a function of the gate length, ∆ y has been calculated according to the mesh spacing in the transport direction (∆ x ). A fixed number of mesh points is taken in our calculation in the transport direction regardless of the gate length of the considered device, so that ∆ x varies as L G does so. In this particular study, we have chosen ∆ y = 0.05∆ x , which corresponds to ∆ y = 0.01 nm for the device with L G = 7.5 nm. The second step is to determine whether the particle tunnels or not by using a rejection criterion. To do so, a uniformly distributed random number r DT is generated between 0 and 1 and compared to T DT (E x ). If r DT ≤ T DT (E x ), the superparticle will cross the barrier; otherwise, it will turn back suffering a back-scattering. Finally, if the superparticle undergoes S/D tunneling, its motion inside the potential barrier is described using Newton's mechanics considering an inverted potential profile and ballistic transport [18].

Description of the Effective Mass Calculation
To adopt more reasonable conduction band structures in nanoscaled structures, we accurately calculate m e f f by using DFT implemented in QuantumATK tool of Synopsys [19]. Table 2 summarizes the values of the masses for both devices (DGSOI and FinFET) studied here. It is also important to highlight at this point that the lowest energy subband changes from ∆ 2 in the planar transistor to ∆ 4 in the vertical one. Figure 2 shows the difference of the longitudinal (m l ) and transverse (m t ) effective masses calculated as a function of the silicon thickness (T Si ) for the two different confinement orientations. It is clearly shown that the effective masses tend to m bulk for larger T Si . Although these masses (m l and m t ) are included in the 2D MS-EMC tool as input parameters, it is important to analyze the modification of m x (transport mass), m z (confinement mass) and m y (mass in the direction normal to transport). Their expressions are shown in Table 1, and their particular values are grouped in Table 2 for the two T Si values herein considered. In order to study the impact of T Si reduction, the deviations (in %) of m l , m t and their combinations included in Table 1 have been calculated (Figure 3) as 100 · |m bulk − m e f f |/m e f f . It is interesting to mention that the deviation in m t is more noticeable than that of (m t + m l )/2, which corresponds to m x in the S/D tunneling model for the fundamental subband of the planar and vertical devices, respectively. In particular, they drop from ∼35% (DGSOI devices) to ∼15% (FinFETs) for T Si = 3 nm and from ∼17.5% (DGSOI devices) to ∼2.5% (FinFETs) for T Si = 5 nm.  Table 2. Effective masses (m e f f ) considering the DGSOI and FinFET devices herein studied with silicon thickness T Si = 3-5 nm using DFT simulations included in QuantumATK of Synopsys [19]. Notice that m x and m z are the transport and confinement masses, respectively, m y is the mass in the direction normal to transport, m 0 is the free electron mass, and the subindex of ∆ represents the degeneracy factor associated with the conduction band valley.  Table 2 as a function of the silicon thickness (T Si ) for DGSOI ((100) confinement orientation) as well as FinFET ((011) confinement orientation) devices.

Comparison of MS-EMC with S/D Tunneling Models vs. NEGF
The I D vs. V GS characteristics obtained from ballistic simulations of the DGSOI and FinFET devices at V DS = 500 mV with gate length ranging from 5 nm to 15 nm are shown in Figure 4. Four types of simulations are displayed: (1) the NEGF approach in the NESS tool, (2) the MS-EMC tool without any type of tunneling, and the MS-EMC tool with the S/D tunneling module using (3) T WKB (E x ) and (4) T DT (E x ). In order to attain a good I ON /I OFF behavior, the work function for the devices with Lg = 5 nm was chosen to be 5 eV rather than 4.385 eV, as for the rest of cases. In general, the I D vs. V GS characteristics were shifted to have the same threshold current (I TH ), showing similar I ON in all the cases. The particular values of I TH as a function of the gate length are included in Figure 4 too. Regarding the OFF region, where S/D tunneling was more noticeable, we can reach the following conclusions when the MS-EMC results are compared against NEGF. First, the simulation without any tunneling reduced I OFF substantially due to the absence of particles below the barrier. Second, there was an overestimation of I OFF when the tunneling probability was calculated by T WKB (E x ), specially for L G ≤ 10 nm. Third, when the tunneling probability was chosen as T DT (E x ), the current was comparable to NEGF, showing a reduction of I OFF . In particular for L G = 7.5 nm, it was really similar to NEGF because, as anticipated earlier, the parameter ∆ y included in Equation (10) was calibrated against the NEGF results. Fourth, the S/D tunneling was important in ultra-scaled devices with L G ≤ 10 nm due to the dimensions of the potential barrier. Consequently, for L G = 15 nm, there was almost no difference in the OFF current among the different MS-EMC cases. Moreover, the inherent statistical nature of the MC method also manifested in Figure 4a by the fluctuations in the subthreshold regime for the simulation without any tunneling module. Figure 5 shows the average number of electrons affected by S/D tunneling for the simulations considered in Figure 4. In general, the drain current in Monte Carlo simulators was calculated by the spatial average of the electron current along the channel. Therefore, the number of electrons located inside the potential barrier due to the S/D tunneling model contributed to the increase of the total current. On the other hand, as depicted in Figure 5, the T DT (E x ) probability reduced the number of electrons crossing the potential barrier compared to the T WKB (E x ) case and thus there was a reduction of I OFF . It is also worth to mention that the number of electrons affected by S/D tunneling approached the same value at high V GS regime for both approaches (T WKB (E x ) and T DT (E x )) and for both devices regardless of the gate length. The subthreshold swing (SS) is one of the main parameters used to determine the behavior of electronic devices in the OFF region. In practice, the best MOSFET implementations cannot reduce SS < 60 mV/dec. For the SS calculation, we have considered the I D decade between 10 −3 mA/µm and 10 −2 mA/µm (or between 10 −2 mA/µm and 10 −1 mA/µm) in order to avoid the stochastic noise of the Monte Carlo simulations at very low V GS . Figure 6 shows the SS difference (∆SS) between MS-EMC and NEGF simulations. In general, and for all values of V DS , we have reached the following three conclusions. First, ∆SS was negative for the MC simulation without any tunneling due to its lower I OFF compared to NEGF case. Second, ∆SS tended to zero for the FinFET device with L G ≥ 7.5 nm, showing the excellent agreement between both approaches for that device. Third, for L G = 15 nm, ∆SS was again negative due to the lower I OFF of the different MS-EMC simulations.

Impact of the Effective Mass Choice
In general, the utilization of m e f f instead of m bulk results in a shift of the I D vs. V GS characteristics ( [20]). Accordingly, we have focused on the study of this impact on the threshold voltage shift (∆V TH ) calculated as the difference of V TH using m e f f and m bulk (Figure 7). V TH has been calculated in this work using the constant drain current method [21]. In this section, non-equilibrium simulations (including acoustic phonon, optical, phonon, surface roughness, and Coulomb scattering mechanisms) have been considered using the 2D MS-EMC tool with the three possible combinations: (1) without any tunneling module, (2) with S/D tunneling using T WKB (E x ), and (3) with S/D tunneling using T DT (E x ). Four observations can be made based on the ∆V TH curves displayed in Figure 7. First, ∆V TH was positive for the ultra-scaled devices (L G = 5 and 7.5 nm) because the use of m e f f increases the current, whereas the opposite trend (∆V TH < 0) was observed for devices with L G ≥ 10 nm. Second, ∆V TH was reduced for thicker devices (T Si = 5 nm) because m e f f tends to m bulk when T Si increases. Third, a similar behavior is shown in Figure 7 for the simulations without any tunneling and with T DT (E x ). However, when the tunneling probability was calculated using T WKB (E x ) the difference between using m e f f instead of m bulk is greater due to the overestimation of the superparticles experiencing S/D tunneling. This effect became more relevant when the device size is reduced. In fact, this influence was significant for L G = 5 nm at all V DS and it was extended to longer devices as V DS increased. Fourth, the impact of the effective mass choice was smaller in the FinFET compared to the DGSOI device due to the lower deviation of its effective transport mass (m x ) shown in Figure 3.

Conclusions
This work presents the quantum enhancement of a semi-classical 2D MS-EMC simulator and its application to DGSOI transistors and FinFETs. It has been demonstrated as a useful tool for the optimization of devices targeting sub-10 nm nodes thanks to its higher computational efficiency. Two different approaches to consider S/D tunneling in MC are described in this work and their results with FinFET and DGSOI are compared to those from NEGF formalism. One of these models needs to be calibrated against quantum transport simulations. Results obtained from the MS-EMC code show an excellent agreement with the NEGF simulations in the subthreshold region. The impact of realistic effective masses, calculated from first principles, on electron transport has also been studied by means of MS-EMC simulations. Our findings suggest that effective masses variation alters in a significant way the tunneling probability in the subthreshold regime, in agreement with reported results in the literature.