The DD formalism for carrier transport has been the main workhorse in the TCAD industry for many decades. It is indispensable for simulating bulk CMOS transistors and relatively larger devices where a more sophisticated approach is neither desired nor practical.

In NESS, we have implemented the DD module using a finite volume discretization scheme for the current continuity equation following the Scharfetter-Gummel approach ^[26] using Bernoulli functions. The 3D current continuity equation is self-consistently solved with the 3D Poisson equation until convergence. Different mobility models are included in the current continuity equation. Convergence for potential and charge is reached when the max norm of the difference between two successive Gummel iterations reaches the preset criteria. At present, we have included doping dependence of the mobility using the Masetti model ^[27]. The transverse and longitudinal electric field (,, respectively) dependence of the mobility has been included by means of the Yamaguchi model ^[28] and the Caughey-Thomas ^[29] model, respectively. As examples, simulation results for a conventional bulk MOSFET with channel length of 25nm for two oxide thicknesses are shown in Figure 6(a) for low and high drain bias conditions considering constant bulk mobility. In Figure 6(b), we have shown the cumulative impact of the mobility models on the transfer characteristics for a nanowire transistor with a circular cross-section of 5nm diameter and 15nm channel length.

A key issue with classical DD simulations is that they cannot capture the quantum confinement effects. A quantum-corrected DD simulator can ensure a correct charge profile in the device at a fraction of the computational cost of a full quantum simulator. We have developed and implemented Schrödinger equation-based quantum-corrected DD approach in NESS ^[17]. For this, we first self-consistently solve the 2D Schrödinger equation in planes perpendicular to transport and 3D Poisson equation in the whole device. The 3D quantum charge is calculated using a top of the barrier approach ^[30], summing over all subbands and valleys. At convergence, the quantum charge density(n_Q ) is used to calculate a quantum correction term where N_C is the conduction band density of states, T is temperature, k_B is Boltzmannconstant, and q is the electronic charge ^[31,32]. This term is then used to generate a corrected potential which (instead of the classical potential obtained from the Poisson equation) is used as a driving force in the continuity equation. This is repeated until the charge and the potential converge. The quantum correction can either be fixed for a bias point (for low drain voltage) or can be updated in each Gummel iteration. It is worthwhile to note that this approach does not use any fitting parameters in the quantum correction procedure unlike the density gradient or the effective potential method.

The quantum-corrected DD remedies the deficiency of the classical DD charge profile as can be seen in Figure 7(a) for a nanowire FET with an 5nm × 5nm square cross-section. Further, in contrast to the classical DD, current-voltage characteristics obtained using quantum-corrected DD display the shift in threshold voltage due to quantum confinement with is in an excellent match with the result obtained from ballistic NEGF as shown in Figure 7(b).

The KG solver implemented in NESS provides accurate electron mobility at low-field near-equilibrium conditions ^[33,34]. It combines the quantum effects based on the 1D multi-subband scattering rates of the most relevant scattering mechanisms in confined channels ^[19] and the semi-classical BTE by applying the KG formula within the relaxation time approximation ^[11]. In the first step, the NEGF module of NESS is used to extract the electron densities, subband levels (E_l ), and the corresponding wavefunctions (ξ_l ) at the cross-section area of a gated NW in the presence of a low electric field in the transport direction (the long-channel device approximation).

In the second step, the 1D rates for the dominant scattering mechanisms in silicon are calculated using the parameters from the first step. The scattering rates are directly derived from the Fermi Golden Rule, using the time-dependent perturbation theory and assuming that the transitions between two states occur instantaneously. In this paper, we present two of the implemented scattering mechanisms:

Acoustic (Ac) phonon scattering is considered to be elastic and within the short-wave vector limit. Its equivalent equation from an initial subband l and a final subband l’ is:

(1)

where is the acoustic deformation potential, is the material density, is the reduced Planck’s constant, is the speed of sound, m_l is the electron effective mass in the transport direction, are vectors normal to the transport direction, θ represents the step function, is the kinetic energy for a wavevector with magnitude k, is the energy separation between subbands l and l’, and .

Optical (Op) phonon scattering takes into account g-type and f-type transitions (intra- and inter-valley transitions, respectively) and the energies of the different branches of the optical deformation potential are considered constant (as used in most of the standard approaches). Accordingly, the optical phonon scattering rate for the phonon mode j can be written as:

(2)

where

(3)

with

(4)

Here, n_j is the phonon number, ω_j is the phonon energy, D_Op,j is the optical deformation potential, and m_v (m_v’ ) is the transport effective mass of the initial(final) valleys, respectively.

In the third step, the mobility () for the scattering mechanism i and subband l is calculated considering the semi-classical simulation of the transport properties of a 1D electron gas using the BTE within the relaxation time approximation ^[11] as a function of the relaxation time (), the 1D density of states (g_l ), the Fermi-Dirac function (f₀ ), and the 1D electron concentration (N_l ):

(5)

In the fourth step, we calculate in two strategies the total mobility for the l subband (μ^l ): (1) it is calculated as a function of the individual mobilities associated with each scattering mechanism () using the Matthiessen rule (); and (2) the scattering rates of all mechanisms are directly added to avoid the Matthiessen rule and thereby μ^l is computed using Equation (5). The former strategy is of special interest for devices with large cross sections because the error induced by the Matthiessen rule in narrower devices is comparable to MS-MC and NEGF approaches. Finally, the average mobility of a NW structure is calculated accounting for all the subbands: . The advantage of both semi-classical alternatives in comparison to purely quantum transport simulations is that the rates are individually computed and then combined, reducing dramatically the computational cost.

Figure 8 shows the scattering rates and mobility for cylindrical Si NWs, which main parameters are summarized in Table 1. The results from the KG module have been compared to the results of an external to NESS 1DMC simulator ^[10], where the mobility is extracted after applying a small constant electric field by fitting the average velocity versus field dependence. In general, the 1DMC and KG scattering rates for the lowest subband of the 3nm nanowire (Figure 8(a)) are in very good agreement especially at low energy levels, the most relevant region which determines the accuracy of the low-field mobility calculations. Moreover, the phonon-limited mobility computed with both approaches shows a very good agreement as a function of the line density (Figure 8(b)) for a 3nm and 6nm circular NW with ⟨100⟩, and as a function of the NW widths (Figure 8(c)) at a fixed line density for ⟨100⟩ and ⟨110⟩ orientations.

The so-called NEGF formalism, which is derived based on the Keldysh technique ^[35], is a widely applied framework for analyzing the electronic transport in non-equilibrium many-body systems. This method allows a quantum treatment of charge transport in order to capture quantum phenomena such as tunneling, coherence, and particle-particle interactions in mesoscopic and nanoscale devices. We obtain the charge density, potential profile, and the current flow in the system by performing a self-consistent solution of the Poisson equation and the NEGF transport equations in coupled-mode space (CMS). We can either consider diffusive transport by switching on the acoustic- and/or optical-phonon scattering ^[36,37] to enable electron-phonon (e–ph) interactions within the self-consistent Born approximation (SCBA) or neglect them to investigate merely the ballistic transport ^[13]. Moreover, we can simulate 2D planar structures such as DGSOI ^[38], and the NEGF solver implemented in NESS also allows calculation of the band-to-band tunneling by using the Flietner model to compute the current in heterostructures with a direct bandgap ^[39].

Adopting the notation of Reference [14], we will summarize the main concepts required to understand the NEGF formalism. Having the system in a steady state, the retarded, advanced, and lesser/greater Green’s functions in real space representation read:

(6)

(7)

where, h, and represent the one-particle Hamiltonian, and the retarded (lesser/greater) self-energies accounting for electrons interactions with their surroundings, respectively. The charge at position r and the current take the forms:

(8)

. (9)

Here indicates the matrix elements of the Hamiltonian (lesser Green’s function) between the basis states in layer l +1 (l) and l (l +1) ^[12,40].

Before considering the e–ph interactions, let us briefly discuss the CMS approximation. The single-particle Hamiltonian in the EM approximation can be expressed as:

. (10)

We can obtain the CMS representation by projecting each diagonal block of on a subspace spanned by some chosen eigenmodes of ^[41] . The transformation matrix is unitary in the limit where all the transverse modes are selected and, consequently, the CMS Hamiltonian is exactly equivalent to the real space Hamiltonian. On the other hand, the CMS Hamiltonian with few chosen modes is equal to the full rank EM Hamiltonian on the chosen subspace, as it reproduces by construction the exact selected EM sub-bands and their wavefunctions. Therefore, CMS offers the possibility to reproduce the effect of roughness or ionized impurities if enough modes are chosen. In this approximation, the matrix elements between the modes i and j read

. (11)

To study the diffusive transport, the interactions of the electrons with phonons is implemented within the NESS via introducing the corresponding self-energies in real space ^[42,43]:

(12)

(13)

Where v, q, and refer to the electronic valley index, optical phonon with energy , and the Bose-Einstein occupation number. The coupling constant of acoustic phonons , and the coupling strength of e–ph interaction are obtained from the deformation potential theory ^[44]. The retarded component of the self-energy stemming from e-ph interactions may be expressed as:

(14)

Its CMS counterpart has the same form, whereas the real space self-energies are replaced by the CMS ones. Following the same notation as in Equation (11), and assuming that the self-energies are local in both space and time, the self-energies due to e-ph interactions in CMS representation read ^[45]:

, (15)

(16)

We can define the total retarded (lesser) self-energy as , where refers to the impact of electron exchange with the contacts ^[14].

In Figure 9(a), we show the ON-state-current spectrum resulting from the simulations for a 3nm × 3nm square NW transistor including scattering processes at V_DS = 0.6V. The tunnelling current reaches high values up to 30 μA/eV. Overall current damping is observed due to acoustic phonons and energy relaxation of carriers as they approach the drain due to optical phonons emission. Figures 9(b) and (c) show the I_DS−V_GS characteristics for a n-type 3nm×3nm square Si NW assuming ballistic and dissipative NEGF transport simulations. Figure 9(b) shows the results with L_G=20nm using the classical DD and the NEGF modules (Ac, g-type optical Op, and a combination of both phonon scattering mechanisms) at low drain voltage, whereas Figure 9(c) compares the results with L_G=10nm and L_G=20nm at V_DS=0.6V. More results from the NEGF module of NESS are presented in ^{[1,14,15,16,18,20,22,39]}.