A review of selected topics in physics based modeling for tunnel field-effect transistors

Kane proposed his BTBT theory in his seminal paper in 1959 [13]. The paper made use of perturbation theory to write the Schrödinger equation as a sum of interband and intraband coupling elements. Here we provide a different, simpler proof of his results based on the Landauer conduction formula and WKB approximation. We start with a so-called local model, that is a model that assumes a constant electric field along the tunneling path and an uniform structure.

For the derivation we start by writing the general expression of the current density using the Landauer's conduction formula [20]:

$$\begin{eqnarray}J & = & \displaystyle \frac{e}{\pi {\hslash }A}{\displaystyle \int }_{-\infty }^{\infty }\displaystyle \sum _{{{\bf{k}}}_{\perp }}T(E,{{\bf{k}}}_{\perp })[{f}_{v}-{f}_{c}]{\rm{d}}{E}\\ & = & \displaystyle \frac{e}{4{\pi }^{3}{\hslash }}{\displaystyle \int }_{-\infty }^{\infty }{\displaystyle \int }_{A}T(E,{{\bf{k}}}_{\perp })[{f}_{v}-{f}_{c}]{\rm{d}}{{\bf{k}}}_{\perp }{\rm{d}}{E},\end{eqnarray} \tag{ 2 }$$

where e is the positive electron charge, ℏ is the reduced Plank constant, A is the area of the device in the direction normal to tunneling, ${{\bf{k}}}_{\perp }$ is the transverse wave-vector, $T(E,{{\bf{k}}}_{\perp })$ is the transmission probability and f_v, f_c are the local occupation functions. For simplicity, hereafter we set f_v = 1 and f_c = 0, even if this implies that equation (2) gives non-zero current when zero bias is applied at the tunneling junction [21]; this problem is solved for example in [14]. The f_c and f_v will appear explicitly in the models in the section 2.2, but here we remain consistent with [13] and set f_v = 1 and f_c = 0. We now convert the tunneling current per unit energy in a generation rate according to equation (1), hence following the physical picture illustrated in figure 2. In particular, we evaluate the integral in equation (2) by moving to polar coordinates; we further assume that the transmission depends on the energy E and only on the magnitude ${k}_{\perp }$ of the wave-vector ${{\bf{k}}}_{\perp }$ and thus obtain

$$\begin{eqnarray}&&{G}_{{\rm{T}}}(E)=\displaystyle \frac{e| F| }{4{\pi }^{3}{\hslash }}2\pi {\int }_{0}^{+\infty }{\rm{T}}(E,{k}_{\perp }){k}_{\perp }{{\rm{d}}{k}}_{\perp },\end{eqnarray} \tag{ 3 }$$

where the tunneling probability in equation (3) can be expressed by using the WKB approximation [22]:

$$\begin{eqnarray}&&T(E,{k}_{\perp })=\displaystyle \frac{{\pi }^{2}}{9}\exp \left(-2{\int }_{{x}_{i}}^{{x}_{f}}\mathrm{Im}({k}_{x}){\rm{d}}{x}\right)\end{eqnarray} \tag{ 4 }$$

and the ${\pi }^{2}/9$ factor is discussed in [23]. In order to evaluate the integral in equation (4) we need an E versus ${\bf{k}}$ relation valid in the energy gap, that is connecting the valence band to the CB. Moreover, since direct BTBT is an elastic process, total transverse momentum and energy need to be conserved. We employ here the Kane's two-band dispersion relation:

$$\begin{eqnarray}&&{E}^{{\rm{valence}}}(k)=\displaystyle \frac{{E}_{G}}{2}+\displaystyle \frac{{{\hslash }}^{2}{k}^{2}}{2{m}_{0}}-\displaystyle \frac{1}{2}\sqrt{{E}_{G}^{2}+\displaystyle \frac{{E}_{G}{{\hslash }}^{2}{k}^{2}}{{m}_{r}}},\end{eqnarray} \tag{ 5a }$$

$$\begin{eqnarray}&&{E}^{{\rm{conduction}}}(k)=\displaystyle \frac{{E}_{G}}{2}+\displaystyle \frac{{{\hslash }}^{2}{k}^{2}}{2{m}_{0}}+\displaystyle \frac{1}{2}\sqrt{{E}_{G}^{2}+\displaystyle \frac{{E}_{G}{{\hslash }}^{2}{k}^{2}}{{m}_{r}}},\end{eqnarray} \tag{ 5b }$$

$$\begin{eqnarray}&&{E}^{{\rm{gap}}}(\kappa )=\displaystyle \frac{{E}_{G}}{2}-\displaystyle \frac{{{\hslash }}^{2}{\kappa }^{2}}{2{m}_{0}}\pm \displaystyle \frac{1}{2}\sqrt{{E}_{G}^{2}-\displaystyle \frac{{E}_{G}{{\hslash }}^{2}{\kappa }^{2}}{{m}_{r}}},\end{eqnarray} \tag{ 5c }$$

where m₀ is the rest electron mass, ${m}_{r}^{-1}={({m}_{c}^{-1}+{m}_{v}^{-1})}^{-1}$, $k=| ({k}_{x},{{\bf{k}}}_{\perp })|$ is the magnitude of the wave vector. Inside the gap, since k_x is imaginary (and thus ${k}_{x}^{2}=-\mathrm{Im}{({k}_{x})}^{2}$), it is easier to express the energy as a function of ${\kappa }^{2}=-{k}^{2}=[\mathrm{Im}{({k}_{x})}^{2}-{k}_{\perp }^{2}]$. The ± sign in equation (5) indicates that inside the gap we have two possible energy values for given κ.

We now set x  =  0 at E  =  0, that is at the point where the top of the valence band crosses the constant energy line that we use as reference. Assuming the electric field F is uniform and along x direction, the energy is given by ${E}^{{\rm{gap}}}=e| F| x$. Neglecting the ${{\hslash }}^{2}{\kappa }^{2}/(2{m}_{0})$ term in equation (5c), the imaginary part of k_x inside the bandgap can be written as:

$$\begin{eqnarray}&&\mathrm{Im}({k}_{x})=\sqrt{\displaystyle \frac{{m}_{r}}{{E}_{G}{{\hslash }}^{2}}}\sqrt{{E}_{G}^{2}+{E}_{G}\displaystyle \frac{{{\hslash }}^{2}{k}_{\perp }^{2}}{{m}_{r}}-4{\left(e| F| x-\displaystyle \frac{{E}_{G}}{2}\right)}^{2}}.\end{eqnarray} \tag{ 6 }$$

After determining the classical turning points of the valence and CBs, which deviate from x_i = 0 and ${x}_{f}={{eE}}_{G}/| F|$ due to non-null parallel momentum ${{\bf{k}}}_{\perp }$, the transmission coefficient can be obtained by evaluating the integral over x in equation (4), that leads to

$$\begin{eqnarray}&&T({k}_{\perp })=\displaystyle \frac{{\pi }^{2}}{9}\exp \left(-\pi \displaystyle \frac{\sqrt{{m}_{r}}{E}_{G}^{\tfrac{3}{2}}}{2{\hslash }e| F| }\right)\exp \left(-\displaystyle \frac{\pi {\hslash }{k}_{\perp }^{2}}{2e| F| }\sqrt{\displaystyle \frac{{E}_{G}}{{m}_{r}}}\right).\end{eqnarray} \tag{ 7 }$$

By substituting equation (7) in equation (3) we find that the integral over ${k}_{\perp }$ can be evaluated analytically and finally reach the well-known formula [13]:

$$\begin{eqnarray}&&{G}_{{\rm{T}}}=\displaystyle \frac{{e}^{2}| F{| }^{2}\sqrt{{m}_{r}}}{18\pi {{\hslash }}^{2}\sqrt{{E}_{G}}}\exp \left(-\pi \displaystyle \frac{\sqrt{{m}_{r}}{E}_{G}^{\tfrac{3}{2}}}{2{\hslash }e| F| }\right).\end{eqnarray} \tag{ 8 }$$

The local models, such as the one presented above, tend to grossly overestimate the current when the electric field rapidly changes over small distances. In order to overcome this issue within the semiclassical modeling framework, the so-called non-local models were introduced [24, 25]. One main goal of these models was to account for the creation of carriers at either side of the tunneling barrier (see also figure 2), in contrast to the local models which unphysically assign electron and hole generation at the same position [26].

As an illustrative example, a non-local direct BTBT model based on Kane's two band dispersion can be obtained using an analytical expression for the imaginary wave-vector κ inside the forbidden gap derived from equation (5) as:

$$\begin{eqnarray}\kappa (x) & = & \displaystyle \frac{1}{{\hslash }}\sqrt{{m}_{r}{E}_{G}(1-{\alpha }^{2}(x))}\\ \alpha (x) & = & -\displaystyle \frac{{m}_{0}}{2{m}_{r}}+2\sqrt{\displaystyle \frac{{m}_{0}}{2{m}_{r}}\left(\displaystyle \frac{E-{E}_{V}(x)}{{E}_{G}}-\displaystyle \frac{1}{2}\right)+\displaystyle \frac{{m}_{0}^{2}}{16{m}_{r}^{2}}+\displaystyle \frac{1}{4}},\end{eqnarray} \tag{ 9 }$$

where the energy E is inside the gap, namely we have ${E}_{V}(x)\leqslant E\leqslant {E}_{V}(x)+{E}_{G}$.

Making use of equation (4) and assuming that only small ${k}_{\perp }$ values result in significant tunneling probability, we expand $\mathrm{Im}({k}_{x})$ for small $[{k}_{\perp }/\kappa ]$ as $\mathrm{Im}({k}_{x})=\sqrt{{\kappa }^{2}+{k}_{\perp }^{2}}\,\simeq \kappa +\tfrac{{k}_{\perp }^{2}}{2\kappa }$ and find:

$$\begin{eqnarray}T & = & \displaystyle \frac{{\pi }^{2}}{9}\exp \left(-2{\displaystyle \int }_{{x}_{i}}^{{x}_{f}}\mathrm{Im}({k}_{x}){\rm{d}}{x}\right)\\ & = & \displaystyle \frac{{\pi }^{2}}{9}\exp \left(-2{\displaystyle \int }_{{x}_{i}}^{{x}_{f}}\kappa {\rm{d}}{x}\right)\exp \left(-{k}_{\perp }^{2}{\displaystyle \int }_{{x}_{i}}^{{x}_{f}}\displaystyle \frac{{\rm{d}}{x}}{\kappa }\right),\end{eqnarray} \tag{ 10 }$$

where we implicitly assume that the physical system is uniform along the y and z direction. In principle the integration extremes x_i and x_f should depend on ${k}_{\perp }$. A further simplification is thus introduced by evaluating x_i and x_f as the classical turning points corresponding to ${k}_{\perp }=0$. Inserting the expression for T into equation (3) and evaluating the integral over ${k}_{\perp }$, we have:

$$\begin{eqnarray}&&{G}_{{\rm{T}}}(E)={\left|\displaystyle \frac{{{\rm{d}}{E}}_{V}}{{\rm{d}}{x}}\right|}_{{x}_{i}}\displaystyle \frac{e}{36{\hslash }}{\left({\displaystyle \int }_{{x}_{i}}^{{x}_{f}}\displaystyle \frac{{\rm{d}}{x}}{\kappa }\right)}^{-1}\\ &&\,\times \left[1-\exp \left(-{k}_{m}^{2}{\displaystyle \int }_{{x}_{i}}^{{x}_{f}}\displaystyle \frac{{\rm{d}}{x}}{\kappa }\right)\right]\times \exp \left(-2{\displaystyle \int }_{{x}_{i}}^{{x}_{f}}\kappa {\rm{d}}{x}\right)({f}_{c}-{f}_{v}),\end{eqnarray} \tag{ 11 }$$

where we have now reintroduced the local occupation functions at the beginning and at the end of the tunneling path. The G_T from equation (11) corresponds to a generation of holes at position x_i and a generation of electrons at position x_f. The electric field F used in equation (1) is here replaced by the derivative of the valence band energy at the beginning of the tunneling path. By considering different energies E in the range $[{E}_{\min },{E}_{\mathrm{MAX}}]$ between minimum value of E_C(x) and the maximum value of E_V(x), equation (11) provides the generation rate of electrons and holes along the whole structure. In this respect, the parameter k_m represents the maximum value of ${k}_{\perp }$ and it is given by:

$$\begin{eqnarray}&&{k}_{m}^{2}=\min \left(\displaystyle \frac{2{m}_{v}({E}_{\mathrm{MAX}}-E)}{{{\hslash }}^{2}},\displaystyle \frac{2{m}_{c}(E-{E}_{\min })}{{{\hslash }}^{2}}\right).\end{eqnarray} \tag{ 12 }$$

Equation (11) is precisely the formula used in the dynamic-path non-local model implemented in some commercial simulators [25].

A similar expression has been used in [27], where the integral of the wave-vector along the tunneling path is performed by solving the equation of motion of the imaginary energy dispersion in the gap using the Monte Carlo method.

We will return to the discussion of non-local tunneling models in section 2.4.

In the presence of quantization, different modeling approaches can be followed depending on the different alignment between tunneling and quantization directions. In fact, in a planar device, tunneling can be in-line with the quantization direction, as in the EHBTFET [28], or it can be transverse to the quantization direction. In this paragraph we give more details for the case with tunneling aligned with quantization because the picture changes more dramatically compared to the transverse quantization case, for which expressions similar to equation (11) have been proposed, where only the prefactor is different due to the different dimensionality of ${{\bf{k}}}_{\perp }$ [29, 30]. We will return to this point in section 2.4.

The models we are reviewing here consider BTBT as a post-processing calculation after the electrostatics has been determined by the self-consistent solution of the closed boundary Schrödinger and Poisson equations (quantum mechanical approach), or the semi-classical continuity equations. This is unlike the NEGF approach or similar quantum transport formalisms (see section 2.5), where transport and electrostatics are inherently coupled and the Schrödinger equation is solved with open boundary conditions along the transport direction.

For tunneling along the quantization direction, the integral over the total energy in equation (2) is replaced by a discrete sum since the energy spectrum of the carrier gas consists of a set of discrete subband levels. Using the Fermi's golden rule approach, one can write the following general expression:

$$\begin{eqnarray}J & = & \displaystyle \frac{4\pi e}{{\hslash }}\displaystyle \sum _{n}\displaystyle \sum _{\alpha ^{\prime} =\{{\rm{HH}},{\rm{LH}}\},n^{\prime} }| {M}_{\alpha ^{\prime} ,n^{\prime} }^{n}{| }^{2}{\rm{JDOS}}({E}_{{\rm{T}}})\\ & & \times {\rm{\Theta }}({E}_{n^{\prime} \alpha ^{\prime} }-{E}_{n{\rm{\Gamma }}})[{f}_{c}({E}_{{\rm{T}}})-{f}_{v}({E}_{{\rm{T}}})],\end{eqnarray} \tag{ 13 }$$

where the summation runs over all possible electron–hole subband pairs (n = electrons, ${n}^{\prime}$ = hole). Note that simultaneous conservation of total energy and transverse momentum (since direct BTBT is an elastic process) results in a single tunneling energy E_T for each pair of subbands (see figure 3). The joint density of states (JDOS) for the case of tunneling between two 2D carrier gases is given by

$$\begin{eqnarray}{\mathrm{JDOS}}_{2{\rm{D}}}({E}_{{\rm{T}}}) & = & \displaystyle \frac{\overline{m}}{4\pi {{\hslash }}^{2}}{\rm{\Theta }}({E}_{n^{\prime} \alpha ^{\prime} }-{E}_{n{\rm{\Gamma }}})\hspace{0.5cm}\mathrm{with}\\ {E}_{{\rm{T}}} & = & {E}_{n{\rm{\Gamma }}}+\displaystyle \frac{{{\hslash }}^{2}{k}_{\perp }^{2}}{2{m}_{e}}={E}_{n^{\prime} \alpha ^{\prime} }-\displaystyle \frac{{{\hslash }}^{2}{k}_{\perp }^{2}}{2{m}_{{\rm{h}}}}.\end{eqnarray} \tag{ 14 }$$

Here $\overline{m}$ is defined as $\overline{m}=2{m}_{{\rm{e}}}{m}_{{\rm{h}}}/({m}_{{\rm{e}}}+{m}_{{\rm{h}}})$ and we see how E_T is determined by energy and parallel momentum conservation.

Zoom In Zoom Out Reset image size

There are two main approaches to calculate the interband coupling coefficient ${M}_{\alpha ^{\prime} ,n^{\prime} }^{n}$ in the presence of quantization.

One is based on Bardeen's transfer Hamiltonian approach [31], where the matrix element for BTBT between two 2D carrier gases can be evaluated as [32]:

$$\begin{eqnarray}&&| {M}_{\alpha ^{\prime} ,n^{\prime} }^{n}{| }^{2}={{\hslash }}^{2}\displaystyle \frac{{E}_{G}^{{\rm{\Gamma }}}}{2\overline{m}}{\delta }_{{k}_{\perp },{{k}^{\prime }}_{\perp }}| {\psi }_{n}{\xi }_{\alpha ^{\prime} ,n^{\prime} }{| }^{2}{C}_{0}(\theta ),\end{eqnarray} \tag{ 15 }$$

where ${E}_{G}^{{\rm{\Gamma }}}$ is the energy gap at the Γ point, ${\delta }_{{k}_{\perp },{{k}^{\prime }}_{\perp }}$ is the Kronecker symbol, ${\psi }_{n}$ and ${\xi }_{\alpha ^{\prime} ,n^{\prime} }$ are the envelope wave-functions for the two different carrier gases and the product $| {\psi }_{n}\,{\xi }_{\alpha ^{\prime} ,n^{\prime} }|$ is somewhat arbitrarily taken at the point where it is maximum. An important difference compared to the bulk case is the inclusion of the term ${C}_{0}(\theta )$, where the angle θ represents the direction of the hole state. In [32], for example, the angle θ is estimated as

$$\begin{eqnarray}&&{\cos }^{2}\theta =\displaystyle \frac{{E}_{{\rm{h}}}^{z}}{{E}_{{\rm{h}}}^{\perp }+{E}_{{\rm{h}}}^{z}},\end{eqnarray} \tag{ 16 }$$

where E_h^z and ${E}_{{\rm{h}}}^{\perp }$ are effective kinetic energies along the quantization and transverse directions, respectively. Expressions for ${C}_{0}(\theta )$ have been given in [33, 34]. This term needs to be introduced in the presence of quantization since it accounts for the dependence of the coupling element on the direction of the electric field. In fact, ${\psi }_{n}$ and ${\xi }_{\alpha ^{\prime} ,n^{\prime} }$ are the envelope wave-functions, whereas the coupling between the electron and hole states is ruled by the underlying Bloch functions. Due to the symmetry properties of the atomic-like orbital functions, this dependence on the angle θ does not appear in the bulk case since for each envelope function we perform an average over all directions of the underlying Bloch functions. In the 2D gas, instead, the average takes place only over the plane normal to quantization, which introduces a dependence on θ.

A second approach to evaluate ${M}_{\alpha ^{\prime} ,n^{\prime} }^{n}$ in equation (13) is similar to the perturbation method employed by Kane, where the presence of the electric field couples states in the conduction and valence band. This methodology has been presented in [35], and then recently employed to 2D–2D tunneling in [36]. The matrix element is given by:

$$\begin{eqnarray}&&{M}_{\alpha ^{\prime} ,n^{\prime} }^{n}=e\sqrt{\displaystyle \frac{{{\hslash }}^{2}}{4\overline{m}{E}_{G}^{{\rm{\Gamma }}}}}{\delta }_{{k}_{\perp },{k}_{\perp }^{\prime} }\int {\psi }_{n}^{* }(z)\sqrt{{C}_{0}(\theta )}| F(z)| {\xi }_{\alpha ^{\prime} ,n^{\prime} }(z){\rm{d}}{z},\end{eqnarray} \tag{ 17 }$$

where F(z) is the electric field, that is allowed to vary only along the quantization direction z.

Figure 4 compares the currents of a Ge quantum-well diode as obtained with either equation (17) or equation (15) for the interband coupling. Results are similar, however equation (17) can be extended quite easily to describe physical systems where the electric field is non-uniform in more than one dimension (i.e. planar TFETs [36]), or quantization is present in more than one dimension (i.e. nanowire TFETs [37]). Reference [37] proposes a model for direct tunneling in bulk, 2D and 1D carrier gasses. It also presents a derivation of the non-local model for direct BTBT in bulk structures similar to our derivation in section 2.1.

Zoom In Zoom Out Reset image size

An important point to note for tunneling in quantized gases is the large asymmetry between the effective masses of the real and the imaginary dispersion for quantized holes. While quantization typically favors heavy hole (HH) subbands to contribute to the current, the calculation of the imaginary energy dispersion in the gap shows that the effective mass for the imaginary band is actually much closer to the light hole (LH) mass in the bulk crystal [38]. This implies that the interband tunneling matrix element in the presence of quantization is close to the one in the corresponding bulk material even if LHs do not come into play [39].

Band structure effects, such as the anti-crossing described above and the influence of quantization on the energy dispersion in the gap, are naturally accounted for in the ${\bf{k}}\,\cdot \,{\bf{p}}$ approach. We will discuss the simulation approach based on the NEGF with a ${\bf{k}}\,\cdot \,{\bf{p}}$ Hamiltonian in section 2.5. It is however worth mentioning that the ${\bf{k}}\,\cdot \,{\bf{p}}$ can be employed also in non-self-consistent quantum models. A relevant example is given in [40, 41], where the ${\bf{k}}\,\cdot \,{\bf{p}}$ method and the corresponding envelope wave-functions have been used with a quantum-transmitting-boundary formalism to study the transmission probability in hetero-junction III–V TFETs. The methodology in [40, 41] can be labeled as non-self-consistent because the potential profile is taken from TCAD simulations and then the BTBT current is calculated, but the possible influence on the potential profile of the charge produced by BTBT is not accounted for.

As an example of results that can be obtained with the models described in this section, we show in figure 5 the I–V characteristic of an EHBTFET with In${}_{0.53}$Ga${}_{0.47}$As channel. In the EHBTFET we have two independent gates that create an electron and a hole inversion layer at the top and bottom of a thin film semiconductor [28], so that BTBT occurs between such two inversion layers. We assume that the potential profile is essentially uniform along the channel direction x, so that we can develop our calculations in a single vertical slice along z. The simulation methodology employed for the results in figure 5 starts by solving Schrödinger equation for electrons and holes for closed boundary conditions self-consistently with the Poisson equation. When computing the charge, electrons and holes are assumed to be at equilibrium with respectively the drain and source contacts: the hole Fermi-level E_Fp is set by the source contact, while electron Fermi-level E_Fn is set by the drain contact and is equal to ${E}_{{\rm{Fp}}}-{{eV}}_{{\rm{DS}}}$. The BTBT current is then computed as a post-processing using equations (13), (17), where the hole envelope functions are modified to account for the anti-crossing between LH and HH bands described above. The bottom of figure 5 shows the subband energy versus gate bias. We see that the current bumps in the top graph correspond to specific subband alignments which occur when the gate voltage increases and moves the electron subbands to lower energy; the hole subbands are also slightly affected by the gate voltage due to coupling between the top and the back interface. In particular, we observe a first bump in the current when e1 (lowest electron subband) crosses the hh1 (highest HH subband). The second bump takes place when e2 aligns with hh1. On the other hand, the alignments between e1 and hh2 and between e1 and lh1 do not produce any current bump because they occur at energy values outside the window between E_Fn and E_Fp, so that the term $[{f}_{c}-{f}_{v}]$ is practically null.

Zoom In Zoom Out Reset image size

Phonon-assisted tunneling occurs in indirect gap semiconductors whose CB minima do not lie at the Γ point: Si and Ge are prominent examples. Keldysh first proposed an expression describing the phonon-assisted BTBT by making use of perturbation theory [42]. Later on, one of the first non-local BTBT models was proposed by Tanaka for phonon-assisted and direct tunneling [24]. In this model, the interband coupling due to electron–phonon interaction in the deformable ion model is used for interband coupling elements of the Wannier formula. The wavefunctions are obtained patching the plane wave solutions (in classically allowed regions) with the exponentially decaying components (in classically forbidden regions) by using the WKB approximation. Compared to the case of direct tunneling, in PAT the energy dispersion in the gap can be obtained from the imaginary branch of the energy relation of the two individual bands, i.e. the phonon connects two branches that can be unequivocally identified as belonging to the conduction or to the valence band.

A recent model that is applicable to non-uniform potential profiles and to a carrier gas with different dimensionality has been proposed by Vandenberghe et al [43]. The approach makes use of the diagonal elements of the spectral functions that, for a 2D gas described by the parabolic effective mass approximation (EMA), can be written as

$$\begin{eqnarray}&&{A}_{c\alpha }(z,E)=\displaystyle \frac{{g}_{\alpha }{m}_{{xy},\alpha }}{{{\hslash }}^{2}}\displaystyle \sum _{{{\bf{k}}}_{\perp }}{\psi }_{{{\bf{k}}}_{\perp }\alpha }^{* }(z){\psi }_{{{\bf{k}}}_{\perp }\alpha }(z){\rm{\Theta }}(E-{E}_{{{\bf{k}}}_{\perp }\alpha }),\end{eqnarray} \tag{ 18a }$$

$$\begin{eqnarray}&&{A}_{v\alpha }(z,E)=\displaystyle \frac{{g}_{\alpha }{m}_{{xy},\alpha }}{{{\hslash }}^{2}}\displaystyle \sum _{{{\bf{k}}}_{\perp }}{\psi }_{{{\bf{k}}}_{\perp }\alpha }^{* }(z){\psi }_{{{\bf{k}}}_{\perp }\alpha }(z){\rm{\Theta }}({E}_{{{\bf{k}}}_{\perp }\alpha }-E),\end{eqnarray} \tag{ 18b }$$

where z is the quantization direction and the structure is uniform along x and y directions, ${\rm{\Theta }}(E)$ is the step function (corresponding to a step-like 2D DOS), ${g}_{\alpha }$, ${m}_{{xy},\alpha }$ are respectively the valley degeneracy and transverse effective mass for the subband α. Since wave-functions ψ can be assumed to be independent of ${{\bf{k}}}_{\perp }$, the spectral functions essentially reduce to summations over the 2D DOS functions of the different subbands weighted by the probability distribution of carriers.

The total phonon-assisted tunneling current is evaluated by first calculating the energy dependent transmission probability:

$$\begin{eqnarray}&&{T}_{v}^{{\rm{abs}},{\rm{em}}}(E)=\displaystyle \frac{{D}_{{\rm{ph}}}^{2}{{\hslash }}^{2}}{2\rho {E}_{{\rm{ph}}}}\displaystyle \sum _{\alpha ,\alpha ^{\prime} }\int {A}_{v\alpha ^{\prime} }(z,E){A}_{c\alpha }(z,E\pm {E}_{{\rm{ph}}}){\rm{d}}{z},\end{eqnarray} \tag{ 19 }$$

where D_ph and E_ph are the deformation potential and the energy of the phonon, respectively, and ρ is the mass density of the material. The phonon-assisted BTBT current is finally calculated as the summation of phonon emission and absorbtion terms:

$$\begin{eqnarray}{J}_{{\rm{ph}}-} & = & -\displaystyle \frac{e}{\pi {\hslash }}\displaystyle \int (\{{f}_{v}(E)[1-{f}_{c}(E-{E}_{{\rm{ph}}})][{n}_{{\rm{B}}}({E}_{{\rm{ph}}})+1]\\ & & -{f}_{c}(E-{E}_{{\rm{ph}}})[1-{f}_{v}(E)]{n}_{{\rm{B}}}({E}_{{\rm{ph}}})\}{T}_{v}^{{\rm{em}}}(E)\\ & & +\{{f}_{v}(E)[1-{f}_{c}(E+{E}_{{\rm{ph}}})]{n}_{{\rm{B}}}({E}_{{\rm{ph}}})\\ & & -{f}_{c}(E+{E}_{{\rm{ph}}})[1-{f}_{v}(E)][{n}_{{\rm{B}}}({E}_{{\rm{ph}}})+1]\}\\ & & \times {T}_{v}^{{\rm{abs}}}(E)){\rm{d}}{E},\end{eqnarray} \tag{ 20 }$$

where f_c and f_v are the occupation functions of the conduction and valence bands, respectively, and ${n}_{{\rm{B}}}({E}_{{\rm{ph}}})$ is the Bose–Einstein distribution for the phonons. For device structures such as EHBTFETs, f_c and f_v can be taken as Fermi–Dirac functions with respectively electron and hole Fermi levels (whose difference is set to ${{eV}}_{{\rm{DS}}}$), which is the assumption used also to obtain the results of figure 5 discussed at the end of section 2.2.

In [43] it is shown that, if the modeling framework described above is applied to the bulk crystal (using the corresponding spectral functions) and for a uniform electric field, then one obtains the same expression for the tunneling generation rate as in [42]. For a non uniform field, instead, it can be shown that, by introducing suitabe approximations for the envelope wave-functions, the formalism in equation (20) provides an expression for non-local tunneling in bulk systems similar to the dynamic-path, non-local-phonon-assisted model used in TCAD [25], and similar also to the model derived in [24].

Models for BTBT are present in many commercial TCAD tools. Besides the local models based on the expressions in [13, 42] (usually reformulated to assure a zero tunneling generation rate at equilibrium [14]), commercial TCAD tools have also recently implemented non-local models similar to equation (11) for direct tunneling and phonon-assisted tunneling [25], based on the approach developed in [24].

As discussed in section 2.1, these models have been derived for bulk structures, where the carriers are not confined. The effect of carrier quantization is manifold. First of all it increases the effective energy gap and changes the density of states (DOS) according to the dimensionality of the carrier gas. In addition, quantization changes the dispersion relationship inside the gap. Many efforts are being devoted at present to account for these effects in TCAD models.

As for the increase of the effective energy gap, in [44] the energy profile of the CB close to the channel/dielectric interface is modified to mimic the presence of the electron subband splitting. More precisely the part of the profile ${E}_{C}(z)$ that is lower than the lowest electron subband E_e0 is set to ${E}_{C}(z)={E}_{{\rm{e}}0}$. As an alternative approach, an effective gap obtained from the solution of the 1D Schrödinger equations for electrons and holes has been also used in [45, 46].

A similar methodology has been used in semi-classical Monte Carlo simulations [30] that implements non-local models for direct and phonon assisted BTBT as in [25]. First the 1D Schrödinger equation for electrons and holes is solved in each slice along the channel, then the conduction (E_C) and valence (E_V) band profiles are modified so as to suppress E_C values below the energy of the lowest electron subband E_e0 and the E_V values above the highest hole subband E_h0. The non-local model with the effective gap correction was found to match quite well the full-quantum non-self-consistent model described in section 2.3 [43]. More recently, the authors of [47] proposed a methodology based on the rejection of those tunneling paths that involve states in the CB below E_e0 or states in the valence band above E_v0. Thanks to the use of physical models interfaces, PMI, [25], this approach can be plugged directly in the TCAD, with no need for an external solver of the Schrödinger equation.

The effect of quantization on the energy dispersion inside the gap of III–V materials has been also analyzed using the ${\bf{k}}\cdot {\bf{p}}$ method [29]. A simple approach to modify the parameters of Kane's formula has been proposed, that requires only the knowledge of the effective gap (distance between E_e0 and E_h0), because the effective mass m_r has an almost linear dependence on the effective gap. As a result the prefactor in Kane's formula scales as ${({E}_{{\rm{e}}0}-{E}_{{\rm{h}}0})}^{2}$, whereas the exponential term remains unchanged. The expression of the prefactor is further modified based on the dimensionality of the carrier gases [29]. Similar expressions for low dimensionality gases have been reported also in [30].

The non-local models implemented in TCAD tools require the identification of a suitable tunneling path, defined as the line in real space over which performing, for example, the integral of equation (11) for the case of direct tunneling. In fact, when deriving equation (11) we asssumed a purely 1D electric field profile, but in a realistic device simulation the electric field may follow quite complex patterns, so that the integration direction must be carefully defined and selected. The effect of different choices for the tunneling path has been analyzed in [26], going from simple horizontal paths to a more complicated dynamic path [25], where different directions of the path are identified in the different points of the discretization mesh based on the gradient of the valence band energy. This algorithm leads to more precise estimates of the BTBT current in complex device architectures, but has the disadvantage of not being available in AC simulations [25]. Another algorithm making use of Newton's law in the forbidden gap to estimate the tunneling path has been proposed in [48], essentially similar to the Monte-Carlo procedure in [27].

Overall, the models implemented in commercial TCADs allow for an investigation of the design space for TFETs with a numerical efficiency that is not attainable with full-quantum tools (not even with the non-self-consistent models described in the previous section). With TCAD tools it is thus possible to analyze devices with a relatively large size, as most of the actual TFETs fabricated so far. Furthermore, models for trap-assisted-tunneling (TAT) are also available [25]. TAT has a strong effect on the subthreshold characteristic of many fabricated TFETs. Non-local models have been employed to analyze the subthreshold behavior of fabricated silicon [44, 49], as well as III–V TFETs [50, 51].

In some circumstances, however, the numerical efficiency of the TCAD tools comes at the cost of a limited accuracy and predictive capability. In fact, an inherently quantum effect is approximated by point-like generation of electrons and holes and the generation rate comes from approximated WKB integrals. Besides the approximations in the modeling of BTBT, the transport in TCAD tools is often described by using a drift-diffusion approach, whose results in terms of internal physical quantities are questionable in nanoscale transistors. For example, in nanoscale MOSFETs described with drift-diffusion, the carriers essentially move at the saturation velocity, whereas quasi-ballistic transport is expected to take place in nanoscale TFETs. Transport of the generated carriers affect the electrostatic of the device, that may influence the electric field in the BTBT region and thus the overall tunneling generation rate. In [26] a non-local BTBT model, similar to the models available in commercial TCAD simulators, has been implemented in a Monte-Carlo simulator. Comparison between different self-consistent schemes allowed authors to analyze the impact on the transistor electrostatics and current of the transport of generated carriers. In silicon TFETs this effect was found to be very limited, due to the low BTBT rates and low concentrations of generated carriers. However, in III–V hetero-junction TFETs where the BTBT rate is much higher (see section 4) this may not remain valid.

As a final remark, we notice that the equations for non-local BTBT (such as equation (11)) apply to the transitions from one specific branch of the valence band to one minimum of the CB. Consequently, in principle one should evaluate the BTBT rates considering all the possible combinations of valence and CBs, and each pair of bands should have specific values of the tunneling parameters entering the BTBT rate equations. We will return to this point in section 3, where we will analyze BTBT in strained silicon TFETs.

The NEGF formalism is considered in many aspects the most general and rigorous method to simulate quantum transport in electronic devices and it is thus particularly suited to study TFETs. The NEGF method solves the Schrödinger equation with open boundary conditions, that is electron transport through a quantum system connected to external electronic reservoirs and it can also consider the impact of different sources of scattering. Several reviews and textbooks describe the details of the theory [52–54], hence here we will only recall basic equations and discuss their use in the simulation of TFETs.

2.5.1. Coherent transport regime

When dealing with quantum transport problems, the most natural choice is to express the Green's functions in the real-space representation [55]. In particular, either by using a local atomic basis (as in tight-biding method) [56–59], or by discretizing with finite difference (or other discretization methods) an EMA or a ${\bf{k}}\,\cdot \,{\bf{p}}$ Hamiltonian [60–62], all operators take the form of matrices, and the retarded Green's function matrix G is obtained by solving the problem

$$\begin{eqnarray}&&({EI}-H-{{\rm{\Sigma }}}_{{\rm{S}}}-{{\rm{\Sigma }}}_{{\rm{D}}})G=I,\end{eqnarray} \tag{ 21 }$$

where I is the identity matrix, H is the Hamiltonian matrix describing the device and ${{\rm{\Sigma }}}_{{\rm{S}}},{{\rm{\Sigma }}}_{{\rm{D}}}$ are the so-called retarded self-energy matrices accounting for the injection/absorption of carriers from external leads, as illustrated in figure 6. The determination of the self-energies requires the knowledge of the Hamiltonian terms describing the coupling between the device and the contacts (or leads), as well as the Green's function of the leads themselves, which are described as semi-infinite periodic systems. The Green's function of the leads can be numerically evaluated with various approaches [63], among which we would like to mention the recursive Sancho–Rubio [64], and the eigenvalue method [65, 66]. From the knowledge of the Green's function and the contact self-energies, macroscopic quantities like DOS, carrier concentration and source-to-drain current I_DS can be readily evaluated. In the absence of inelastic interactions, the current is expressed by means of the Landauer–Büttiker formula [67, 68]

$$\begin{eqnarray}&&{I}_{{\rm{DS}}}=\displaystyle \frac{e}{\pi {\hslash }}\int {\rm{d}}{E}\ T(E)[{f}_{{\rm{S}}}(E)-{f}_{{\rm{D}}}(E)],\end{eqnarray} \tag{ 22 }$$

where ${f}_{{\rm{S}},{\rm{D}}}$ is the Fermi–Dirac function at source (S) and drain (D) and the transmission probability T(E) is evaluated from the retarded Green's function as

$$\begin{eqnarray}&&T(E)=\mathrm{Tr}[{{\rm{\Gamma }}}_{{\rm{S}}}(E)\,G(E)\,{{\rm{\Gamma }}}_{{\rm{D}}}(E){G}^{\dagger }(E)],\end{eqnarray} \tag{ 23 }$$

where $\mathrm{Tr}[\cdots ]$ is the trace operation, and ${{\rm{\Gamma }}}_{{\rm{S}},{\rm{D}}}\,={\rm{i}}[{{\rm{\Sigma }}}_{{\rm{S}},{\rm{D}}}-{{\rm{\Sigma }}}_{{\rm{S}},{\rm{D}}}^{\dagger }]$ is the broadening function of the source and drain leads.

Zoom In Zoom Out Reset image size

Equations (22) and (23) assume a coherent propagation of the electron wave function inside the device, hence they do not describe the phase-breaking interaction with phonons, that may be particularly important for indirect bandgap semiconductors or for PAT processes mediated by defects or interface states. Such a coherent quantum transport formalism has been used for instance to simulate the electric performance of carbon nanotube (CNT) [69, 70], graphene based TFETs [71] and InAs nanowire TFETs [57, 62]. Moreover, this approach can account for some sources of elastic scattering, such as surface roughness [72, 73], because these interactions do not break quantum-phase coherence and their effect can be included in the real-space electrostatic potential. We notice, however, that, in case a variability analysis is performed, this approach requires the generation of several different samples in order to accumulate a significant device statistics.

2.5.2. Possible forms for the matrix H

2.5.3. Reduction of the computational burden

2.5.4. Alternative approaches for coherent transport

The advantage of the NEGF method with respect to the QTB method is primarily in the possibility to include the effect of inelastic scattering mechanisms, such as the electron–phonon interaction [84]. This is obtained by introducing the concept of lesser(greater)-than Green's functions ${G}^{\lt (\gt )}$, which are associated to the electron(hole) statistics, whereas the retarded Green's fucntion G provides the information about the available electronic states and the carrier dynamics [54]. More precisely, in the steady-state regime the main diagonal of the lesser-than GF matrix (that is the ${G}^{\lt }({{\bf{r}}}_{i},{{\bf{r}}}_{i})$ in a real-space representation), is used to compute the electron concentration $n({{\bf{r}}}_{i})$, the main diagonal of the greater-than GF matrix, ${G}^{\gt }({{\bf{r}}}_{i},{{\bf{r}}}_{i})$, gives the hole concentration $p({{\bf{r}}}_{i})$ and the main diagonal of the retarded GF matrix, $G({{\bf{r}}}_{i},{{\bf{r}}}_{i})$, is linked to the local density of states (LDOS) $\rho ({{\bf{r}}}_{i})$.

The lesser(greater)-than self-energy describing the phonon scattering is evaluated as ${{\rm{\Sigma }}}_{\mathrm{PH}}^{\lt (\gt )}={D}^{\lt (\gt )}{G}^{\lt (\gt )}$ [84], where ${D}^{\lt (\gt )}$ is the phonon lesser(greater)-than Green's function. The problem is here most often simplified by assuming that the phonon bath is at equilibrium, so that ${D}^{\lt (\gt )}$ can be written in terms of the Bose–Einstein distribution function ${n}_{{\rm{B}}}(\omega )={[\exp ({\hslash }\omega /{k}_{{\rm{B}}}T)-1]}^{-1}$, where ${\hslash }\omega$ is the phonon energy. In this case the lesser-than and the greater-than self-energies for the jth phonon branch read [84]

$$\begin{eqnarray}{{\rm{\Sigma }}}_{\mathrm{PH},j}^{\lt }(E) & = & {| {M}_{j}| }^{2}\{{G}^{\lt }(E-{\hslash }{\omega }_{j}){n}_{{\rm{B}}}({\omega }_{j})\\ & & +{G}^{\lt }(E+{\hslash }{\omega }_{j})[{n}_{{\rm{B}}}({\omega }_{j})+1]\},\end{eqnarray} \tag{ 24a }$$

$$\begin{eqnarray}{{\rm{\Sigma }}}_{\mathrm{PH},j}^{\gt }(E) & = & {| {M}_{j}| }^{2}\{{G}^{\gt }(E+{\hslash }{\omega }_{j}){n}_{{\rm{B}}}({\omega }_{j})\\ & & +{G}^{\gt }(E-{\hslash }{\omega }_{j})[{n}_{{\rm{B}}}({\omega }_{j})+1]\},\end{eqnarray} \tag{ 24b }$$

where M_j is the matrix element expressing the microscopic details of the electron–phonon interaction [85]. Equation (24) describe a dissipative phenomenon and in fact the self-energy at energy E depends on the Green's function at a different energy $E\pm {\hslash }{\omega }_{j}$, where the first term in equations (24) represents a phonon-absorption and the second term a phonon-emission process. A very similar formalism can be used also to consider approximately elastic phonons, such as intra-valley acoustic phonons in the long-wavelength approximation.

2.6.1. Self-consistent Born approximation

Since the retarded phonon self-energy ${{\rm{\Sigma }}}_{\mathrm{PH}}$ has to be expressed in terms of the greater-than and lesser-than self-energies according to the relation ${{\rm{\Sigma }}}_{\mathrm{PH}}-{{\rm{\Sigma }}}_{\mathrm{PH}}^{\dagger }={{\rm{\Sigma }}}_{\mathrm{PH}}^{\gt }-{{\rm{\Sigma }}}_{\mathrm{PH}}^{\lt }$, the retarded and the lesser-than Green's function are non-linearly coupled and have to be calculated self-consistently. This is obtained by solving self-consistently the kinetic equations

$$\begin{eqnarray}&&[{EI}-H-{{\rm{\Sigma }}}_{{\rm{S}}}-{{\rm{\Sigma }}}_{{\rm{D}}}-{{\rm{\Sigma }}}_{\mathrm{PH}}]\,G=I,\end{eqnarray} \tag{ 25 }$$

$$\begin{eqnarray}&&[{EI}-H-{{\rm{\Sigma }}}_{{\rm{S}}}-{{\rm{\Sigma }}}_{{\rm{D}}}-{{\rm{\Sigma }}}_{\mathrm{PH}}]\,{G}^{\lt }=[{{\rm{\Sigma }}}_{{\rm{S}}}^{\lt }+{{\rm{\Sigma }}}_{{\rm{D}}}^{\lt }+{{\rm{\Sigma }}}_{\mathrm{PH}}^{\lt }]{G}^{\dagger }\end{eqnarray} \tag{ 26 }$$

until convergence is reached. This method to include phonon scattering within the NEGF formalism is called self-consistent Born approximation, and it has the remarkable merit that it guarantees current conservation along the device also in the presence of dissipative phenomena [54]. In case of dissipative transport the current can no longer be expressed in terms of the retarded Green's function alone, and the lesser-than function is also necessary [84]. For example, the steady-state current at the lead $L={\rm{S}},{\rm{D}}$ may be written as

$$\begin{eqnarray}&&{I}_{L}=\displaystyle \frac{e}{\pi {\hslash }}\int {\rm{d}}{E}\,\mathrm{Tr}[{{\rm{\Sigma }}}_{L}^{\lt }(G-{G}^{\dagger })+{\rm{i}}{{\rm{\Gamma }}}_{L}{G}^{\lt }]\,,\end{eqnarray} \tag{ 27 }$$

where the broadening function is given by ${{\rm{\Gamma }}}_{L}\,={\rm{i}}[{{\rm{\Sigma }}}_{L}-{({{\rm{\Sigma }}}_{L})}^{\dagger }]$.

The impact of phonon scattering on the transfer characteristics of TFETs within the NEGF formalism has been discussed with more details in [58, 86, 87].

2.6.2. Main difficulties and approximations

Due to the maturity of the silicon CMOS technology, at first the implementation of TFETs in Si has been investigated as a quite natural option. In fact most of the fabrication steps are shared with the ones used for conventional CMOS transistors and this allows for MOSFET-TFET co-integration [93].

A good electrostatic control is mandatory in order to enhance the electric field at the source-channel junction, increase the BTBT rate and obtain steep sub-threshold I_DS versus V_GS characteristic. For this reason TCAD simulations have investigated TFETs with thin SOI films and a double-gate biasing, as well as GAA NWs TFETs [94]. The use of thin silicon films on one side improves the electrostatic control (resulting in a shorter tunneling path), but on the other hand increases the effective gap (see section 2.4). To discuss this trade-off, we plot in figure 7 the trans-characteristic of double-gate Si TFETs as obtained with the multi-subband Monte Carlo described in [30]. The simulator implements a non-local model for BTBT similar to the one in TCAD [25], but with a quantum-corrected tunneling path (QCTP) where the effective gap is modified in each point of the real space according to the solution of the 1D Schrödinger equation (see discussion in section 2.4).

Zoom In Zoom Out Reset image size

We see in figure 7 that the effect of the QCTP is stronger at low V_GS and for thin films. However, the improved electrostatic control obtained with thin films is not offset by the increase of the effective gap, so that for T_Si = 3 nm we observe the best SS and on-current: the minimum point SS is around 20 mV/dec, however the on-current is only a few nA μm⁻¹. The low on-current of Si-TFETs found in simulations is consistent with available experimental data [9, 12, 93, 95]. In fact some experiments have reported currents approaching the μA μm⁻¹, but for very high biases. If extrapolated to supply voltages around 0.5 V, the experimental on-current is still in the few nA μm⁻¹ range.

It must be noted that device structures as the ones analyzed in figure 7 (and in several experimental studies [93, 95]) are dominated by point-tunneling, that is not very efficient since the tunneling direction is normal to the vertical electric field induced and modulated by the gate. Many structures based instead on line-tunneling have been proposed recently [49, 96, 97], where the tunneling direction is aligned with the electric field induced by the gate. One example is the EHBTFET that we have discussed in section 2.2, while a second example is the core–shell nanowire proposed in [98]. However, also these architectures are expected to result in poor on-current if implemented in silicon. For the Si EHBTFET this has been shown both by experiments [99] and by simulations [100]. The simulations for the core–shell nanowire indicate an on-current of about 30 nA μm⁻¹ if one sets I_off to 1 pA μm⁻¹ in figure 2 of [98], and then projects the on-current to a supply voltage below 0.5 V. Also the experimental data for line tunneling in [49] results in an on-current of a few tens nA μm⁻¹ if it is projected to a 0.5 V power supply. Core–shell TFETs using III–V materials, instead, look promising and relatively large on currents have been experimentally reported [101].

Strain modifies the band structure and suitable strain configurations may enhance the BTBT rate. Indeed, experimental data for strained Si and SiGe nanowires have shown promising on-currents [102–105]. The effect of strain on the current of TFETs has been analyzed experimentally in [106], while TCAD simulations of Si TFETs with non-uniform strain profiles have been reported in [107].

In this section we illustrate a worked out example of how to use a commercial TCAD environment [25] to investigate possible strain induced BTBT enhancements in silicon TFETs. The BTBT rate for phonon assisted transitions is given by [25]:

$$\begin{eqnarray}{G}_{{\rm{T}}} & = & | {\rm{\nabla }}{E}_{v}(0)| \cdot {C}_{p}{\rm{\exp }}\left(-2{\displaystyle \int }_{0}^{l}{\kappa }_{{\rm{im}}}{\rm{d}}{x}\right)[{f}_{v}-{f}_{c}]\\ {C}_{p} & = & \displaystyle \frac{{C}_{{\rm{path}}}}{{2}^{6}{\pi }^{2}{E}_{g}}\,\sqrt{\displaystyle \frac{{m}_{c}{m}_{v}}{h\sqrt{2{m}_{r}{E}_{g}}}}{\left({\displaystyle \int }_{0}^{l}\displaystyle \frac{{\rm{d}}{x}}{{\kappa }_{{\rm{im}}}}\right)}^{-1}\\ & & \times \left[1-{\rm{\exp }}\left(-{k}_{v,{\rm{\max }}}^{2}{\displaystyle \int }_{0}^{{l}_{c}}\displaystyle \frac{{\rm{d}}{x}}{{\kappa }_{v}}\right)\right]\\ & & \times \left[1-{\rm{\exp }}\left(-{k}_{c,{\rm{\max }}}^{2}{\displaystyle \int }_{{l}_{c}}^{l}\displaystyle \frac{{\rm{d}}{x}}{{\kappa }_{c}}\right)\right]\\ {\kappa }_{{\rm{im}}} & = & \min ({\kappa }_{c},{\kappa }_{v}),{\kappa }_{c}=\displaystyle \frac{1}{{\hslash }}\sqrt{2{m}_{c}({E}_{C}-E)},\\ {\kappa }_{v} & = & \displaystyle \frac{1}{{\hslash }}\sqrt{2{m}_{v}(E-{E}_{V})},\end{eqnarray} \tag{ 28 }$$

where 0 and l are the coordinate of the beginning and the end of the tunneling path, while l_c is the point where ${\kappa }_{c}={\kappa }_{v}$. The terms ${k}_{v,\max }$ and ${k}_{c,\max }$ have the same meaning as the term k_m in equation (11) and C_path contains the phonon parameters (see discussion below). The reduced mass is ${m}_{r}=2{m}_{c}{m}_{v}/({m}_{c}+{m}_{v})$. We see that the modeling ingredients for the BTBT rate are the effective masses for the valence and CBs, the energy gap and the phonon energy and deformation potential.

Strain changes both the energy gap, inducing splitting between the different valleys of the conduction and valence bands, and modifies the effective masses. With the tool sband included in [25] we computed the effective masses inside the band (real branch) and in the gap (imaginary branch). A good correspondence is found between (real) mass values in literature and sband for both unstrained and strained silicon. For the CB, the masses extracted from the imaginary branches of the energy relation in the gap are almost the same as the masses corresponding to the real energy branch. For the valence band, instead, the masses extracted in the energy gap are remarkably different compared to the values corresponding to real energy branches. More precisely, figure 8 shows that the effective mass of the imaginary branch corresponding to the HH valley has a small mass, similar to the one of the real branch of the LH valley, whereas the imaginary branch of the LH valley has a mass similar to the real branch of the HH valley. These imaginary valence bands can still be reproduced, in first approximation and for small energies, with a parabolic expression, as can be seen in figure 8. The effect of biaxial strain on the valence band parameters is shown in figure 9.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Equation (28) refers to a single type of tunneling transition, for example between LH or HH and one of the minima of the CB. To simulate more realistically the BTBT process in silicon, we used six different tunneling paths, from HH and LH to the three Δ valleys of the CB. All these paths produce different contributions to the generation rate, particularly when strain breaks the symmetry between the CB valleys. For each valley combination we used masses and energy gaps resulting from the sband analysis. In particular, for each path we used the effective masses (m_c and m_v) obtained from sband along the corresponding tunneling direction and corresponding to the imaginary energy relation in the gap. The phonon energy is assumed to be the same as in unstrained Si. The problem with equation (28) is that it does not take into account the fact that the CB minima are anisotropic. In fact, the parameter m_c in equation (28) is used to compute the transmission probability (and in this respect it should be the effective mass inside the gap), but also as a prefactor that accounts for the DOS involved in tunneling. To solve this issue, we set as m_c and m_v the masses inside the energy gap (that is for imaginary branches of the energy relation) from sband, and modified the C_path parameters by the ratio between the DOS masses and the imaginary masses:

$$\begin{eqnarray}&&{C}_{{\rm{path}}}={C}_{{\rm{path}}0}{\left(\displaystyle \frac{{m}_{{\rm{dos}},c}{m}_{{\rm{dos}},v}}{{m}_{c}{m}_{v}}\right)}^{3/2}.\end{eqnarray} \tag{ 29 }$$

Equation (29) implicitly assumes that the phonon deformation potential does not change with strain. The power of 3/2 stems from the fact that, under constant lateral electric field, equation (28) reduces to Keldysh formula with prefactor proportional to ${({m}_{c}{m}_{v})}^{3/2}{m}_{r}^{-5/4}$. While the ${m}_{r}^{-5/4}$ is related to the dispersion relationship in the gap, the term ${({m}_{c}{m}_{v})}^{3/2}$ instead results from the DOS. The correction in equation (29) is thus exact under uniform lateral field and approximated in general cases. The value of C_path0 has been calibrated to reproduce the experimental Si diodes in [9, 108], once the effective masses have been extrated from sband calculations.

To assess the impact of strain on BTBT, we have considered uniform structures with constant lateral electric field and different doping levels. Results are reported in figure 10, where we see that, after proper calibration of the C_path0, the approach consisting of six tunneling paths reproduces the experimental data with better accuracy compared to the default SDevice calibration. Biaxial strain of 2% is shown to enhance the BTBT rate by more than an order of magnitude.

Zoom In Zoom Out Reset image size

The contributions of the different tunneling paths are reported in table 1 for one particular value of F_MAX. For unstrained silicon the tunneling is dominated by the transitions from HH (low tunneling mass) to ${{\rm{\Delta }}}_{y}$ and ${{\rm{\Delta }}}_{z}$ (low tunneling mass along the electric field direction x). In the presence of biaxial tensile stress, the energy of the ${{\rm{\Delta }}}_{z}$ valley is reduced, and transitions between HH and ${{\rm{\Delta }}}_{z}$ are significantly enhanced. The change in the gap and effective masses also enhances the LH to ${{\rm{\Delta }}}_{z}$ transitions.

Table 1.  Contribution of different tunneling paths to the currents in figure 10.

	${J}_{{\rm{unstr}}}/V$	${J}_{{\rm{tens}},1 \% (001)}/V$
Paths	$@{F}_{\mathrm{MAX}}=1.45\,\mathrm{MV}\,{\mathrm{cm}}^{-1}$	$@{F}_{\mathrm{MAX}}=1.29\,\mathrm{MV}\,{\mathrm{cm}}^{-1}$
	$[m{\rm{A}}\,\mu {{\rm{m}}}^{-2}\,{{\rm{V}}}^{-1}]$	$[{\rm{m}}{\rm{A}}\,\mu {{\rm{m}}}^{-2}\,{{\rm{V}}}^{-1}]$
LH $\to \,{{\rm{\Delta }}}_{x}$	$5.5\times {10}^{-13}$	$5.4\times {10}^{-11}$
HH $\to \,{{\rm{\Delta }}}_{x}$	$5.9\times {10}^{-10}$	$5.4\times {10}^{-11}$
LH $\to \,{{\rm{\Delta }}}_{y}$	$9.3\times {10}^{-9}$	$9.0\times {10}^{-10}$
HH $\to \,{{\rm{\Delta }}}_{y}$	$4.5\times {10}^{-7}$	$3.9\times {10}^{-9}$
LH $\to \,{{\rm{\Delta }}}_{z}$	$9.3\times {10}^{-9}$	$1.1\times {10}^{-6}$
HH $\to \,{{\rm{\Delta }}}_{z}$	$4.5\times {10}^{-7}$	$2.9\times {10}^{-6}$

Germanium has been deeply investigated as channel material for TFETs due to the lower band-gap and smaller effective masses compared to Si. In Ge the direct gap is only 0.14 eV larger than the indirect gap between the top of the valence band and the CB minimum at the L point. For this reason direct and phonon-assisted BTBT may both contribute to the overall tunneling current. In this respect, figure 11 shows that the experimental data for the tunneling currents in the Ge diodes collected in [109] can be reproduced quite well by the model for direct BTBT, whereas the contribution of phonon-assisted transitions is one order of magnitude lower.

Zoom In Zoom Out Reset image size

The default calibration in TCAD tools is quite empirical and assumes that phonon-assisted tunneling is dominant in Ge [25]. This point is thoroughly discussed in [110], which proposes a more realistic calibration where Ge is instead dominated by direct tunneling. It can be verified that in uniform structures, the direct BTBT rate provided by the model in [110] is very close to the phonon-assisted BTBT rate obtained with the default calibration in TCAD, whereas the phonon assisted model from [110] provides a much lower BTBT rate.

To assess the effect of different calibrations on realistic devices, we simulated with TCAD the planar device with Ge source and Si channel fabricated in [111]. In figure 12 the experimental results of [111] are reported by black solid line and look promising, with an ${I}_{\mathrm{ON}}$ close to 1 μA μm⁻¹ and a minimum SS significantly below 60 mV/dec at low bias. Simulations using the default SDevice calibration (dashed lines) slightly underestimate the current but are, somewhat surprisingly, in quite good agreement with the experiments. The dotted curves with symbols in figure 12 have been obtained with the calibration proposed in [110]: as can be seen they are significantly lower than the experiments and predict a poor device performance. In particular, as discussed above, the calibration in [110] considers a much smaller phonon assisted BTBT (circles) than the one assumed as default in TCAD. Direct BTBT (triangles) is high, but its onset is at large gate voltage due to the Ge–Si hetero-junction that essentially inhibits direct tunneling from the Ge valence band to the Γ point in the CB of Si. In fact, direct BTBT becomes significant only when the gate voltage is large enough to allow for an essentially vertical tunneling within the Ge source. On the other hand, phonon assisted BTBT takes place between the Ge source and the Si channel.

Zoom In Zoom Out Reset image size

The default calibration in SDevice and the calibration in [110] for BTBT in Ge result in orders of magnitude difference in on-current and also large differences in SS. Notice that the phonon model in [110] has been recently revised in [112]. The new set of parameters (that we used in figure 11) features a significantly higher phonon deformation potential ($7.8\times {10}^{8}\,{\rm{eV}}\,{{\rm{cm}}}^{-1}$ instead of $0.8\times {10}^{8}\,\mathrm{eV}\,{\mathrm{cm}}^{-1}$ in [110]) and slightly different phonon energy (6 meV versus 8.6 meV). It is easy to show that these modifications result in a BTBT rate 200 times higher than using the set proposed in [110]. If we multiply the dotted curve with symbol in figure 12 (obtained with the set in [110]) by a factor of 200, the simulations approach the experiments (not shown), but are still more than one decade below them.

From the analysis above, we can conclude that there are still open issues in the calibration of TCAD BTBT models for Ge TFETs. In fact, the uniform structures used for model calibration are dominated by direct BTBT, whereas in realistic devices phonon-assisted tunneling dominates, and different sets of parameters have thus been proposed in the literature. As a further example of interplay between direct and phonon-assisted BTBT, we consider in figure 13 a Ge EHBTFET. We have employed the non self-consistent quantum model for direct and phonon-assisted BTBT described in sections 2.2 and 2.3 respectively. Due to quantization, the offset between the direct and the indirect gap becomes larger than in bulk Ge. In particular, we can see that the first bump in the current takes place when subband e1, originating from the L-valleys, gets aligned with the heavy-hole subband hh1. A second bump is due to alignment between e1 and hh2. We see direct tunneling only at much higher gate voltages, such that subband ${\rm{\Gamma }}1$ aligns with hh1, but in the same gate voltage range we also have another bump due to phonon-assisted transitions between e2 and hh1.

Zoom In Zoom Out Reset image size

The effect of strain in Ge TFETs has been analyzed in [113] by using quantum-corrected TCAD simulations to estimate the current and ${\bf{k}}\cdot {\bf{p}}$ calculations for the band parameters. It was found that tensile strain significantly enhances the BTBT rate, similarly to the results for strained Si discussed in section 3.2.

Due to their small energy gap and electron effective mass, InAs based TFETs have been extensively simulated with the NEGF approach by using either an atomistic, TB Hamiltonian or a ${\bf{k}}\,\cdot \,{\bf{p}}$ approach.

The TB method was used in [57] to analyze InAs p–i–n single-gate, dual-gate, ultrathin-body and gate-all-around nanowire (GAA NW) devices with a gate length of 20 nm, showing that a reduced subthreshold swing can only be achieved if the electrostatic potential is efficiently controlled by the gate. Similarly, in [58] it is shown that in direct-gap semiconductors, such as InAs, the BTBT is often dominated by a single branch of the energy relation in the gap, and therefore a fairly good estimate of the tunneling rate can be obtained with the WKB approximation.

In order to reduce the computational effort and simulate larger structures within the NEGF method, the ${\bf{k}}\,\cdot \,{\bf{p}}$ approximation has been successfully employed to simulate transport at the Γ point in homojunction InAs TFETs [40, 60, 61, 122, 123], as well as in hetero-junction TFETs such as GaSb/InAs based transistors [124]. In [124] it was shown that ideal InAs NW TFETs with gate length of 17 nm require a lateral width of 7 nm or less to achieve a SS smaller than 60 mV/dec due to the loss of electrostatic integrity at small gate lengths. In such narrow FETs the energy dispersion is strongly altered by quantum confinement compared to the energy dispersion in the corresponding bulk materials. An example is shown in figure 15 reporting the highest valence subband and the lowest conduction subband profiles along the channel, as well as the current spectra for a p–i–n hetero-junction nanowire FET having GaSb in the source and InAs in the channel and drain regions. As can be seen, for a lateral diameter of ${D}_{W}=7$ nm or smaller this material system is no longer broken-gap.

Zoom In Zoom Out Reset image size

The requirement for narrow multi-gate transistors necessary to preserve a strong control of the gate terminal on the channel region results in a significant carrier confinement that in turn degrades the on current delivered by TFETs. Consequently, several material and device design options have been investigated to improve the on current and preserve a steep subthreshold swing.

As a prominent example mechanical strain has been proposed, in analogy to its use as mobility booster in MOSFETs, as an efficient way to improve the on current of InAs TFETs [77]. In [61], the impact of different strain conditions on the bandstructure and ${I}_{{\rm{DS}}}$ versus ${{\bf{V}}}_{{\rm{GS}}}$ characteristics of an InAs nanowire TFET was investigated by using NEGF simulations based on an 8-band ${\bf{k}}\,\cdot \,{\bf{p}}$ Hamiltonian [125], and with a strain interaction matrix whose deformation potentials were taken from [114]. It was shown that both compressive uniaxial and tensile biaxial stress can be used to modify the imaginary branch of the InAs nanowire and reduce the effective gap between the highest valence and the lowest conduction subband. The tensile biaxial strain in the device cross-section was found to be particularly effective in order to reduce the energy gap and enhance BTBT. The introduction of strain, however, implies also a deterioration of the SS due to the change of the valence subband profile in the source region [124], so that specific design options at the source region may be needed to preserve a steep subthreshold characteristic in the presence of strain [126]. The impact of strain on the bandstructure and the turn-on characteristics of III–V based TFETs is still under investigation at the time of writing [127–129].

III–V compounds have also the remarkable merit to form hetero-junctions with interesting properties. In particular, the GaSb/InAs hetero-junction is approximately lattice matched and forms a broken bandgap material system, that is a system where the valence band edge of GaSb is higher than the CB edge of InAs [130]. Such a material system has attracted a lot of interest in the quest for large on currents in tunnel FETs [62, 131].

Heterojunction TFETs have been simulated by using the NEGF approach with both a TB [83, 115, 132, 133], and a ${\bf{k}}\,\cdot \,{\bf{p}}$ Hamiltonian [62, 124, 126, 131, 134]. The use of the GaSb/InAs hetero-junction usually provides a larger ${I}_{\mathrm{ON}}$ compared to InAs homojunction TFETs, but such a potential improvement is partly frustrated by quantum confinement effects that tend to turn the supposedly broken-bandgap GaSb/InAs hetero-junction into an actually staggered band alignment system, as illustrated in figure 15. These effects reduce the advantages of the GaSb/InAs hetero-junction compared to an InAs homo-junction TFET, and also compared to the ideal broken-bandgap, 1D system considered in [135]. According to our results, the GaSb/InAs hetero-junction TFET is still unable to reach the very demanding ITRS specifications for low power applications, which prescribe an I_ON/I_OFF ratio larger than 10⁷ and ${I}_{\mathrm{OFF}}=10$ pA μm⁻¹. Therefore several additional design options have been proposed to increase the on current, including the grading of the Al molar fraction in an AlGaSb source [131], or the insertion and engineering of quantum wells in the source region [78].

Recently an alterative device design based on the concept of polarization-engineered tunnel diodes has been proposed [136], which exploits III-nitrides hetero-junctions [82]. By means of TB NEGF simulations it was shown that interband tunneling can be significantly large in GaN/InN/GaN hetero-junctions, leading to an on current close to 100 μA μm⁻¹ and to a very competitive subthreshold swing.

Non-ideal properties of interfaces are considered the most relevant problem affecting the off-state behavior of TFETs and impeding the experimental realization of devices with a subthreshold swing well below 60 mV/dec. The impact of traps has been first included in NEGF simulations with a phenomenological description in [137], where interface traps were modeled as 0D electrically active states modifying simultaneously the electron transport and the device electrostatics. This was obtained by superimposing cubic potential wells to the CB with different volume size and potential height in order to tune the energy trap levels [138]. It was found that only a few traps can significantly deteriorate the SS of narrow InAs NW FETs with a 5 × 5 nm² cross-section, because traps can act as stepping stones for BTBT and thus enhance the corresponding current, with a tunneling mechanism that can be phonon assisted and thus inelastic. The phonon assisted tunneling also brings along a sensitivity to the temperature of the otherwise temperature independent ${I}_{{\rm{DS}}}$ versus ${{\bf{V}}}_{{\rm{GS}}}$ characteristics [138].

Several papers using a semi-classical approach to describe trap assisted tunneling by using a modified Shockley–Read–Hall formula found that traps can produce a sizable degradation of the SS, and investigated the minimum trap densities necessary for a steep slope behavior [51, 120, 139–141]. In very recent references, for instance, it has been argued that the trap density should be decreased by 100 times with respect to the state of the art values in order to obtain IIIV TFETs with an SS below 60 mV/dec at room temperature [142]. The key role played by traps and defects for the interpretation of the experimental IV characteristics of InAs nanowire TFETs has been clearly underlined in [51]. Analitycal formulas using a detailed balance analysis have been also used to include the effects of traps in the simulation of GaSb/InAs TFETs [143].

Interface traps can be also an important source of device variability in TFETs [144, 145], together with surface roughness [72, 73], dopant fluctuations in the source [146, 147] and work-function variation associated to the granularity of the metal gate [141, 148]. Traps typically influence the variability of ${I}_{\mathrm{OFF}}$ and subthreshold swing much more than the ${I}_{\mathrm{ON}}$ variability [144], because in the on state the direct tunneling is dominant with respect to trap-assisted tunneling.

From a modeling perspective it is interesting to notice that both resonant tunneling FETs and DOS switches based on van der Waals hetero-structures have been originally conceived and described using the Bardeen's transfer Hamiltonian approach [31, 160, 161]. Based on such formalism, the current per unit area can be written as [154]:

$$\begin{eqnarray}&&I={g}_{v}\displaystyle \frac{4\pi e}{{\hslash }}\displaystyle \sum _{{{\bf{k}}}_{{\rm{T}}},{{\bf{k}}}_{{\rm{B}}}}| M({{\bf{k}}}_{{\rm{T}}},{{\bf{k}}}_{{\rm{B}}}){| }^{2}\delta ({E}_{{\rm{B}}}({{\bf{k}}}_{{\rm{B}}})-{E}_{{\rm{T}}}({{\bf{k}}}_{{\rm{T}}}))({f}_{{\rm{B}}}-{f}_{{\rm{T}}}),\end{eqnarray} \tag{ 30 }$$

where g_v is the valley degeneracy, ${{\bf{k}}}_{{\rm{B}}}$, ${{\bf{k}}}_{{\rm{T}}}$ are the wave-vectors respectively in the bottom (source) and top (drain) 2D material, ${E}_{{\rm{B}}}({{\bf{k}}}_{{\rm{B}}})$ and ${E}_{{\rm{T}}}({{\bf{k}}}_{{\rm{T}}})$ denote the corresponding energies, while f_B and f_T are the Fermi occupation functions in the bottom and top layer that depend on the respective Fermi levels E_FB and E_FT. This expression for the current density assumes that the majority carriers of the two 2D materials are at thermodynamic equilibrium with their Fermi levels, where the difference between the Fermi levels, ${{E}}_{\mathrm{FB}}-{{E}}_{{\rm{FT}}}\,=$ eV_DS, is simply set by the drain to source voltage ${V}_{{\rm{DS}}}$.

The reader may notice the similarity between equations (30) and (13) introduced for 2D to 2D tunneling (e.g. EHBTFETs). However in equation (30) the wave-vectors ${\bf{k}}$_B, ${\bf{k}}$_T belong to two different materials and, furthermore, the equation does not enforce the conservation of the in plane wave-vector (as discussed in more details below), whereas the JDOS in equation (13) entails the conservation of both energy and in plane wave-vector in the BTBT process.

The matrix element $M({{\bf{k}}}_{{\rm{T}}},{{\bf{k}}}_{{\rm{B}}})$ expresses the coupling between the two layers and entirely governs the current, but the Dirac function $\delta ({E}_{{\rm{B}}}({{\bf{k}}}_{{\rm{B}}})-{E}_{{\rm{T}}}({{\bf{k}}}_{{\rm{T}}}))$ is also of utmost importance in equation (30), because it enforces energy conservation by prescribing that the sum over ${\bf{k}}$_B, ${\bf{k}}$_T in equation (30) be restricted to the states of the two layers that have the same energy. This implies that the current vanishes if the top of the valence band in the bottom layer is lower than the bottom of the CB in the top layer. Under these circumstances, in fact, no valence band electrons exist in the bottom layer that can transfer to the CB of the top layer via an elastic tunneling process.

As can be seen the current expression in the Bardeen's transfer Hamiltonian approach describes quite effectively the nature of DOS switch of a vdW-TFET. In fact, by modulating the band edge alignment in the top and bottom layer by means of an external gate bias, one can open or close the tunneling window and thus drive the device respectively in the on- or in the off-state.

It is usually assumed that in actual van der Waals hetero-structures several physical mechanisms can occur in the interlayer region that relax the conservation of the in plane wave-vector ${\bf{k}}$ in the tunneling process [157]. This implies that in equation (30) we have significant contributions to the current also for different ${\bf{k}}$_B and ${\bf{k}}$_T values, unlike in the physical picture described in section 2.2 and formally expressed by equation (13). Physical mechanisms that can assist the tunneling and relax momentum conservation may be charged impurities [162], short-range disorder [163], or Moiré patterns that have been observed at the graphene-hBN interface [164, 165]. When the tunneling is assisted by a perturbation potential ${U}_{{\rm{sc}}}(r,z)$ in the interlayer region, the matrix element in equation (30) can be expressed as [154]:

$$\begin{eqnarray}&&M({{\bf{k}}}_{{\rm{T}}},{{\bf{k}}}_{{\rm{B}}})={\displaystyle \int }_{A}{\rm{d}}{r}\displaystyle \int {\rm{d}}{z}\,{\psi }_{{\rm{T}},{{\bf{k}}}_{{\rm{T}}}}^{\dagger }(r,z)\,{U}_{{\rm{sc}}}(r,z)\,{\psi }_{{\rm{B}},{{\bf{k}}}_{{\rm{B}}}}(r,z),\end{eqnarray} \tag{ 31 }$$

where r and z are the in plane and vertical spatial coordinates, while ${\psi }_{{\rm{B}},{{\bf{k}}}_{{\rm{B}}}}$ and ${\psi }_{{\rm{T}},{{\bf{k}}}_{{\rm{T}}}}$ are the electron wave-function respectively in the bottom and top 2D layers.

The wave-functions ${\psi }_{{\rm{T}},{\bf{k}}}(r,z)$ and ${\psi }_{{\rm{B}},{\bf{k}}}(r,z)$ in equation (31) are the Bloch functions of the two 2D layers and they entail a decay in the out of plane direction that has not been made explicit in the expression of the matrix element. The calculation of such a decay is the most important and arguably most difficult problem in the application of Bardeen's transfer Hamiltonian approach to vdW-TFETs, and it is a necessary step if one wants to use equations (30), (31) as a quantitative and predictive model, as opposed to a conceptual tool for an insightful but essentially qualitative analysis. As a matter of fact, however, the decay of the wave-function in the interlayer or in the van der Waals gap has been frequently taken as an adjustable parameter to reproduce experimental data [156, 157].

The problem of the coupling between the wave-functions in the bottom and top layer remains a central and somewhat thorny problem even in different transport formalisms, such as the NEGF modeling based on the empirical TB method. In fact the parameters of a TB Hamiltonian are typically calibrated to reproduce the electronic band-structure, as determined via ab initio density functional theory (DFT) calculations or inferred from experiments. Such a calibration is effective for the in-plane coupling energies of a given 2D crystal, whereas it is far more difficult to determine the vertical coupling between the atomic sites of different materials stacked to form a van der Waals hetero-structure. An interesting alternative to the empirical TB approach is to transform the wave-functions obtained by DFT calculations, that are typically obtained using a plane-wave expansion [166], into maximally localized Wannier functions centered on the ions, which can be accomplished by using the wannier90 package [167]. This effectively produces a TB Hamiltonian matrix with no need of an empirical calibration of the parameters [168]. Such a methodology is becoming quite popular, and it is particularly effective for 2D crystals and vdW-TFETs[169, 170].

Also ab initio calculations for van der Waals hetero-structures, however, encounter some conceptual and practical difficulties, starting right from the possible lattice mismatch of the 2D crystals and the resulting problem in the definition of the super-cell to be used for the analysis of the physical system. As an example, if we consider a van der Waals hetero-structure consisting of ${\mathrm{WTe}}_{2}$ and ${\mathrm{MoS}}_{2}$ monolayers, then the lattice parameter for the two unstrained crystals is about 3.55 Å and 3.19 Å, respectively, so that the super-cell necessary to include an integer number of unstrained unit cells of the individual materials is very large. It is thus a common practice to assume a strain in the two layers so that a matching of the lattice parameters can be obtained. For the case of ${\mathrm{WTe}}_{2}$ and ${\mathrm{MoS}}_{2}$ monolayers, so as to conclude the example, authors have introduced a compressive strain on ${\mathrm{WTe}}_{2}$ and tensile strain on ${\mathrm{MoS}}_{2}$ layer to obtain a commensurate lattice parameter of 3.411 Å in the two materials [171–173]. As can be seen for the ${\mathrm{MoS}}_{2}$ this implies a strain larger than 6%, which alters significantly the bandstructure of an isolated ${\mathrm{MoS}}_{2}$ monolayer compared to the unstrained case. Thus, quite paradoxically, although it has been claimed that the weak van der Waals bonding in the vertical direction may ease the fabrication of high quality hetero-structures even in the presence of a significant lattice mismatch, in most DFT studies of van der Waals hetero-structures a significant strain has been introduced in order to achieve a lattice matching condition.

A cutting-edge methodology in the modeling of van der Waals hetero-structures and vdW-TFETs is the employment of ab initio methods to calculate directly the transmission across the physical structures [174], which is a capability that has become available even in some popular, open-source programs such as Quantum Espresso and Siesta [166, 175].

As of today, some aspects remain very hard to be adequately described in the numerical simulation of van der Waals hetero-structures and vdW-TFETs. As an example, it should be recalled that in most fabrication techniques it is very challenging to obtain an accurate rotational alignment of the two 2D materials. The rotational alignment of the two 2D crystals is crucial for the working principle of resonant tunneling vdW-TFETs [156, 157], whereas it is expected to be less critical for DOS switches targeting a small SS [154]. In any case it is difficult to account for rotational misalignment in transport simulations and even in bandstructure calculations, in fact such a misalignment poses new difficulties in the definition of the super-cell of the physical system, besides the already mentioned problems due to a possible lattice mismatch. Furthermore, similarly to the case of TFETs consisting of conventional 3D semiconductors, even for vdW-TFETs it is difficult to include in a self-consistent transport and device simulation the tails of the DOS in the energy gap, whose effect on the minimum SS of vdW-TFETs can be important.

While the original concept of a van der Waals TFET has been conceived as a device where the current should be fairly constant across the overlap area between the two 2D layers, numerical simulations have shown that the current injection from the bottom to the top layer can be strongly non uniform in the overlap area [170, 172, 173], and very sensitive to the length ${{L}}_{{\rm{ext}}}$ of the extension between the gate and the top 2D layer (see figure 16(b)). In this latter regard, figure 17 reports the ${I}_{{\rm{DS}}}$ versus top gate voltage, ${V}_{{\rm{TG}}}$, characteristic for a van der Waals TFET consisting of a WTe₂ bottom layer (being the source region), and an MoS₂ top layer (being the drain); the simulation approach has been described in [172, 173]. As can be seen the sub-threshold region is strongly influenced by ${{L}}_{{\rm{ext}}}$ and, in particular, for ${{L}}_{{\rm{ext}}}$ = 20 nm the device can attain a subthreshold swing as low as 11 mV/dec (average value for ${I}_{{\rm{DS}}}$ between 10 pA μm⁻¹ and 1 μA μm⁻¹), whereas for smaller ${{L}}_{{\rm{ext}}}$ the swing tends to degrade. This behavior has been explained by a close analysis of the LDOS in the off-state of the device. In particular, it was found that electrons in the valence band of ${\mathrm{WTe}}_{2}$ have energy values that, in the overlap region, correspond to states deep in the energy band gap of the ${\mathrm{MoS}}_{2}$ top layer. Consequently, the most favorable tunneling path from the bottom to the top layer was found to be at the right edge of the bottom layer and the tunneling distance increases by enlarging ${{L}}_{{\rm{ext}}}$ [172, 173]. As a result, the current in the off-state is exponentially suppressed by increasing ${{L}}_{{\rm{ext}}}$, which explains the large sensitivity of ${I}_{{\rm{DS}}}$ to ${{L}}_{{\rm{ext}}}$ in figure 17, and the degradation of SS for small ${{L}}_{{\rm{ext}}}$ values.

Zoom In Zoom Out Reset image size

The edge transport discussed above is an undesired effect in van der Waals TFETs, and it can lead to an ${I}_{{\rm{DS}}}$ that does not scale with the length of the overlap region [172, 173]. While a number of design aspects have been recently analyzed by means of numerical simulations, experiments for van der Waals TFETs are still at a quite embryonal stage. However new experimental data for Esaki diodes or transistors consisting of hetero-structure van der Waals systems is appearing with a steady pace [176–178], and a transistor based of an ${\mathrm{MoS}}_{2}$–Ge hetero-structure has recently shown sub-60 mV/dec subthreshold swing over several orders of magnitude of ${I}_{{\rm{DS}}}$ variation [155].