# Single Particle Transport in Two-dimensional Heterojunction Interlayer Tunneling Field Effect Transistor

###### Abstract

The single particle tunneling in a vertical stack consisting of monolayers of two-dimensional semiconductors is studied theoretically and its application to a novel Two-dimensional Heterojunction Interlayer Tunneling Field Effect Transistor (Thin-TFET) is proposed and described. The tunneling current is calculated by using a formalism based on the Bardeen’s transfer Hamiltonian, and including a semi-classical treatment of scattering and energy broadening effects. The misalignment between the two 2D materials is also studied and found to influence the magnitude of the tunneling current, but have a modest impact on its gate voltage dependence. Our simulation results suggest that the Thin-TFETs can achieve very steep subthreshold swing, whose lower limit is ultimately set by the band tails in the energy gaps of the 2D materials produced by energy broadening. The Thin-TFET is thus very promising as a low voltage, low energy solid state electronic switch.

###### pacs:

## I Introduction

The electronic integrated circuits are the hardware backbone of todays information society and the power dissipation has recently become the greatest challenge, affecting the lifetime of existing portable equipments, the sustainability of large and growing in number data centers, and the feasibility of energy autonomous systems for ambience intelligence Rabaey et al. (2002); Amirtharajah and Chandrakasan (1998), and of sensor networks for implanted monitoring and actuation medical devices Dreslinski et al. (2010). While the scaling of the supply voltage, V, is recognized as the most effective measure to reduce switching power in digital circuits, the performance loss and increased device to device variability are a serious hindrance to the V scaling down to 0.5 V or below.

The voltage scalability of VLSI systems may be significantly improved by resorting to innovations in the transistor technology and, in this regard, the ITRS has singled out Tunnel filed effect transistors (FETs) as the most promising transistors to reduce the sub-threshold swing, SS, below the 60 mV/dec limit of MOSFETs (at room temperature), and thus to enable a further V scaling Group et al. (2011); Seabaugh and Zhang (2010). Several device architectures and materials are being investigated to develop Tunnel FETs offering both an attractive on current and a small SS, including III-V based transistors possibly employing staggered or broken bandgap heterojunctions Zhou et al. (2012); Tomioka et al. (2013); Knoll et al. (2013); Mohata et al. (2011), or strain engineering Conzatti et al. (2012). Even if encouraging experimental results have been reported for the on-current in III-V Tunnel FETs, to achieve a sub 60 mV/dec subthreshold swing is still a real challenge in these devices, probably due to the detrimental effects of interface states Zhou et al. (2012); Pala and Esseni (2013); Esseni and Pala (2013). Therefore, as of today the investigation of new material systems and innovative device architectures for high performance Tunnel FETs is a timely research field in both the applied physics and the electron device community.

In such a contest, two-dimensional (2D) crystals attract increasingly more attention primarily due to their scalability, step-like density of states and absence of broken bonds at interface. They can be stacked to form a new class of tunneling transistors based on an interlayer tunneling occurring in the direction normal to the plane of the 2D materials. In fact tunneling and resonant tunneling devices have been recently proposed Feenstra et al. (2012), as well as experimentally demonstrated for graphene-based transistors Britnell et al. (2013); Britnell et al. (2012a). Furthermore, monolayers of group-VIB transition metal dichalcogenides MX (M = Mo, W; X = S, Se, Te) have recently attracted remarkable attention for their electronic and optical properties Radisavljevic et al. (2011); Wang et al. (2012). Monolayers of transition-metal dichalcogenides (TMDs) have a bandgap varying from almost zero to 2 eV with a sub-nanometer thickness such that these materials can be considered approximately as two-dimensional crystals Gong et al. (2013). The sub-nanometer thickness of TMDs can provide excellent electrostatic control in a vertically stacked heterojunction. Furthermore, the 2D nature of such materials make them essentially immune to the energy bandgap increase produced by the vertical quantization when conventional 3D semiconductors are thinned to a nanoscale thickness, and thus immune to the corresponding degradation of the tunneling current density Jena (2013). Moreover, the lack of dangling bonds at the surface of TMDs may allow for the fabrication of material stacks with low densities of interface defects Jena (2013), which is another potential advantage of TMDs materials for Tunnel FETs applications.

In this paper we propose a two-dimensional heterojunction interlayer tunneling field effect transistor (Thin-TFET) based on 2D semiconductors and develop a transport model based on the transfer-Hamiltonian method to describe the current voltage characteristics and discuss, in particular, the subthreshold swing. In Section II we first present the device concept and illustrate examples of the vertical electrostatic control, then we develop a formalism to calculate the tunneling current. Upon realizing that the subthreshold swing of the Thin-TFET is ultimately determined by the energy broadening, in Sec.II.3 we show how this important physical factor has been included in our calculations. In Sec.II.4 we address the effect of a possible misalignment between the two 2D semiconductor layers, while in Sec.II.5 we derive some approximated, analytical expressions for the tunneling current density, which are useful to gain insight in the transistor operation and to guide the device design. In Sec.III we present the results of numerically calculated current voltage characteristics for the Thin-TFET, and finally in Sec.IV we draw some concluding remarks about the modeling approach developed in the paper and about the design perspectives for the Thin-TFET.

## Ii Modeling of the tunneling transistor

### ii.1 Device concept and electrostatics

The device structure and the corresponding band diagram are sketched in Fig.1, where the 2D materials are assumed to be semiconductors with sizable energy bandgap, for example, transition-metal dichalcogenide (TMD) semiconductors without losing generality Mak et al. (2010); Wang et al. (2012). Both the top 2D and the bottom 2D material is a monolayer and the thickness of the 2D layers is neglected in the modeling of the electrostatics.

The working principle of the tunneling transistor sketched in Fig.1(a) can be explained as follows. When the conduction band edge E of the top 2D layer is higher than the valence band edge E of the bottom 2D layer (see Fig.2(a)), there are no states in the top layer to which the electrons of the bottom layer can tunnel into. This corresponds to the off state of the device. When E is pulled below E (see Fig.2(b)), a tunneling window is formed and consequently an interlayer tunneling can flow from the bottom to the top 2D material. The crossing and uncrossing between the top layer conduction band and the bottom layer valence band is governed by the gate voltages and it is described by the electrostatics of the device.

To calculate the band alignment along the vertical direction of the intrinsic device in Fig.1 we write the Gauss law linking the sheet charge in the 2D materials to the electric fields in the surrounding insulating layers, which leads to

(1) |

where is the capacitance per unit area of top oxide (interlayer, bottom oxide) and is the potential drop across top oxide (interlayer, bottom oxide). The potential drop across the oxides can be written in terms of the external voltages , , and of the energy and defined in Fig.1(b) as

(2) |

where , are fermi levels of majority carriers in the top and bottom layer. , are the electron and hole concentration in the top layer, , the concentrations in bottom layer, , are the electron affinities of the 2D materials, , the workfunctions of the top and back gate and E is the energy gap in the bottom layer. Eq. 2 implicitly assumes that the majority carriers of the two 2D materials are at thermodynamic equilibrium with their Fermi levels, with the split of the Fermi levels set by the external voltages (i.e. EeV), and the electrostatic potential essentially constant in the 2D layers.

Since in our numerical calculations we shall employ a parabolic effective mass approximation for the energy dispersion of the 2D materials, as discussed more thoroughly in Sec.III, the carrier densities can be readily expressed as an analytic function of and Esseni et al. (2011)

(3) |

where is the valley degeneracy.

### ii.2 Transport model

In this section we develop a formalism to calculate the tunneling current based on the transfer-Hamiltonian method Bardeen (1961); Harrison (1961); Duke (1969), as also revisited recently for resonant tunneling in graphene transistors Zhao et al. (2013); Feenstra et al. (2012); Britnell et al. (2013). We start by writing the single particle elastic tunneling current as

(4) |

where is the elementary charge, k, k are the wave-vectors respectively in the bottom and top 2D material, denote the corresponding energies, and are the Fermi occupation functions in the bottom and top layer (depending respectively on E and E, see Fig.1), and is the valley degeneracy. The matrix element expresses the transfer of electrons between the two 2D layers is given by Britnell et al. (2013)

(5) |

where () is the electron wave-function in the bottom (top) 2D layer and is the perturbation potential in the interlayer region.

Eq.5 acknowledges the fact that in real devices several physical mechanisms occurring in the interlayer region can result in a relaxed conservation of the in plane wave-vector k in the tunneling process. We will return to the discussion of in this section.

To proceed in the calculation of we write the electron wave-function in the Bloch function form as

(6) |

where is a periodic function of r and is the number of unit cells in the overlapping area of the two 2D materials. Eq.6 assumes the following normalization condition:

(7) |

where is the in-plane abscissa in the unit cell area and =.

The wave-function is assumed to decay exponentially in the interlayer region with a decay constant Feenstra et al. (2012); Britnell et al. (2013); such a dependence is absorbed in and we do not need to make it explicit in our derivations. It should be noticed that absorbing the exponential decay in recognizes the fact that in the interlayer region the r dependence of the wave-function changes with . In fact, as already discussed Feenstra et al. (2012), while the are localized around the basis atoms in the two 2D layers, these functions are expected to spread out while they decay in the interlayer region, so that the r dependence becomes weaker when moving farther from the 2D layers.

To continue in the calculation of we let the scattering potential in the interlayer region be separable in the form Britnell et al. (2013)

(8) |

where is the in-plane fluctuation of the scattering potential, which is essentially responsible for the relaxation of momentum conservation in the tunneling process.

By substituting Eqs.6 and 8 in Eq.5 and writing , where is a direct lattice vector and is the in-plane position inside each unit cell, we obtain

(9) | |||||

We now assume that F(r) corresponds to relatively long range fluctuations so that it can be taken as approximately constant inside a unit cell, and that, furthermore, the top and bottom 2D layer have the same lattice constant, hence the Bloch functions and have the same periodicity in the r plane. Moreover, for the time being we consider that the conduction band minimum in the top layer and the valence band maximum in the bottom layer are at the same point of the 2D Brillouin zone, so that is small compared to the size of the Brillouin zone and is approximately 1.0 inside a unit cell. These considerations and approximations allow us to rewrite Eq.9 as

(10) |

where the integral in the unit cell has been written for r0 because it is independent of the unit cell.

Consistently with the assumption that k and k are small compared to the size of the Brillouin zone, in Eq.10 we neglect the k (k) dependence of u (u) and simply set uu, uu, where u and u are the periodic parts of the Bloch function at the band edges, which is the simplification typically employed in the effective mass approximation approach Esseni et al. (2011). By recalling that the u and u retain the exponential decay of the wave-functions in the interlayer region with a decay constant , we now write

(11) |

where T is the interlayer thickness and is a k independent matrix element that will remain a prefactor in the final expression for the tunneling current. Since F has been assumed a slowly varying function over a unit cell, then the sum over the unit cells in Eq.10 can be rewritten as a normalized integral over the tunneling area

(12) |

By introducing Eq.11 and 12 in Eq.10 we can finally express the squared matrix element as

(13) |

where and S(q) is the power spectrum of the random fluctuation described by , which is defined as Esseni et al. (2011)

(14) |

Before we proceed with some important integrations of the basic model that will be discussed in Secs.II.3 and II.4, a few comments about the results obtained so far are in order below.

According to Eq.15 the current is proportional to the squared matrix element defined in Eq.11 and decreases exponentially with the thickness interlayer according to the decay constant of the wave-functions. Attempting to derive a quantitative expression for is admittedly very difficult, in fact it is difficult to determine how the periodic functions u and u spread out when they decay in the barrier region and, furthermore, it is not even perfectly clear what potential energy or Hamiltonian should be used to describe the barrier region itself, which is an issue already recognized and thoroughly discussed in the literature since a long time Duke (1969). Our model essentially circumvents these difficulties by resorting to the semi-empirical formulation of the matrix element given by Eq.11, where is left as a parameter to be determined and discussed by comparing to experiments.

It is also worth noting that in our calculations we have not explicitly discussed the effect of spin-orbit interaction in the bandstructure of 2D materials, even if giant spin-orbit couplings have been reported in 2D transition-metal dichalcogenides Zhu et al. (2011). If the energy separations between the spin-up and spin-down bands are large, then the spin degeneracy in current calculations should be one instead of two, which would affect the current magnitude but not its dependence on the gate bias. Our calculations neglected also the possible modifications of band structure in the TMD materials produced by the vertical electrical field, in fact we believe that in our device the electrical field in the 2D layers is not strong enough to make such effects significant Ramasubramaniam et al. (2011).

The decay constant in the interlayer region may be estimated from the electron affinity difference between the 2D layers and the interlayer material Feenstra et al. (2012). Moreover, according to Eq.15 the constant determines the dependence of the current on , so that may be extracted by comparing to experiments discussing such a dependence, which, for example, have been recently reported for the interlayer tunneling current in a graphene-BN system Britnell et al. (2012a).

As for the spectrum S(q) of the scattering potential, in our calculations we utilize

(16) |

where is the correlation length, which in our derivations has been assumed large compared to the size of a unit cell. Eq.16 is consistent with an exponential form for the autocorrelation function of F(r) Esseni et al. (2011), and a similar dependence has been recently employed to reproduce the experimentally observed line-width of the resonance region in graphene interlayer tunneling transistors Britnell et al. (2013). Such a functional form can be representative of phonon assisted tunneling, short-range disorderLi et al. (2011), charged impuritiesYan and Fuhrer (2011) or Moiré patterns that have been observed, for instance, at the graphene-BN interfaceYankowitz et al. (2012); Xue et al. (2011); Decker et al. (2011). We will see in Sec.II.5 that the has an influence on the gate voltage dependent current, which has a neat physical interpretation, hence a comparison to experimental data will be very informative for an estimate of .

### ii.3 Effects of energy broadening

According to Eq.4 and Eq.15 the tunneling current is simply zero when there is no energy overlap between the conduction band in the top layer and the valence band in the bottom layer, that is for . In a real device, however, the 2D materials will inevitably have phonons, disorder, host impurities in the 2D layer and be affected by the background impurities in the surrounding materials, so that a finite broadening of the energy levels is expected to occur because of the statistical potential fluctuations superimposed to the ideal crystal structure Van Mieghem et al. (1991). The energy broadening in 3D semiconductors is known to lead to a tail of the density of states (DoS) in the gap region, that has been also observed in optical absorption measurements and denoted Urbach tail Urbach (1953); Cody (1992). It is thus expected that the finite energy broadening will be a fundamental limit to the abruptness of the turn on characteristic attainable with the devices of this work, hence it is important to include this effect in our model.

Energy broadening in the 2D systems can stem from the interaction with randomly distributed impurities and disorder in the 2D layer or in the surrounding materials Kane (1963); Van Mieghem et al. (1991); Sarma and Vinter (1982), by scattering events induced by the interfaces Knabchen (1995), as well as by other scattering sources. We recognize the fact that a detailed description of the energy broadening is exceedingly complicated due to the many-body and statistical fluctuation effects Ghazali and Serre (1985), and thus resort to a relatively simple semi-classical treatmentKane (1963) Van Mieghem et al. (1991). We start by recalling that the density of states for a 2D layer with no energy broadening is

(17) |

where E(k) denotes the energy relation with no broadening and , are spin and valley degeneracy. In the presence of a randomly fluctuating potential V(r), instead, the DoS can be written as Kane (1963); Van Mieghem et al. (1991)

(18) |

where P(v) is the distribution function for V(r) (to be further discussed below), and we have used the definition in Eq.17 to go from the first to the second equality.

Comparing Eq.18 to Eq.17, we see that the (E) of the system in the presence of broadening can be calculated by substituting the Dirac function in Eq.17 with a finite width function P(v), which is the distribution function of V(r) and it is thus normalized to one.

In order to include the energy broadening in our current calculations, we rewrite the tunneling rate in Eq.4 as

(19) |

and note that, consistently with Eq.18, the energy broadening can be included in the current calculation by substituting with . By doing so the tunneling rate becomes

(20) |

where we have introduced an energy broadening spectrum that is defined as

(21) |

where and is the potential distribution function due to the presence of randomly fluctuating potential in the top and the bottom layer, respectively.

On the basis of Eq.20, in our model for the tunneling current we accounted for the energy broadening by using in all numerical calculations the broadening spectrum defined in Eq.21 in place of . More precisely we used a Gaussian potential distribution for both the top and the bottom layer

(22) |

which has been derived by Evan O.Kane for a broadening induced by randomly distributed impurities Kane (1963), in which case can be expressed in terms of the average impurity concentration.

Quite interestingly, for the Gaussian spectrum in Eq.22 the overall broadening spectrum defined in Eq.21 can be calculated analytically and reads

(23) |

Hence also S has a Gaussian spectrum, where and are the broadening energies for the top and bottom 2D layer, respectively.

### ii.4 Rotational misalignment and tunneling between inequivalent extrema

The derivations in Sec.II.2 assumed that there is a perfect rotational alignment between the top and the bottom layer and that the tunneling occurs between equivalent extrema in the Brillouin zone, that is tunneling from a to a extremum (or from to extremum). We now denote by the angle expressing a possible rotational misalignment between the two 2D layers (see Fig.3), and still assume that the top 2D crystal has the same lattice constant as the bottom 2D crystal. The principal coordinate system is taken as the crystal coordinate system in the bottom layer, and we denote with , the position and wave vectors in the crystal coordinate system of the top layer (with r, k being the vectors in the principal coordinate system). The wave-function in the top layer has the form given in Eq.6 in terms of , , hence in order to calculate the matrix element in the principal coordinate system we start by writing r, k, where is the rotation matrix from the bottom to the top coordinate system, with being the matrix going from the top to the bottom coordinate system and denoting the transpose of the matrix . The rotation matrix can be written as

(24) |

in terms of the rotational misalignment angle .

Consistently with Sec.II.2 we set uu, uu, where u, u are the periodic part of the Bloch function respectively at the band edge in the top and bottom layer. We then denote with the wave-vector at the conduction band edge in the top layer (expressed in the top layer coordinate system), and with the wave-vector at the valence band edge in the bottom layer (expressed in the principal coordinate system); the derivations in this section account for the fact that and may be inequivalent extrema (i.e. ).

By expressing and in the principal coordinate system we can essentially follow the derivations in Sec.II.2 and write the matrix element as

(25) | |||||

where qkk and we have introduced the vector

(26) |

Eq.25 is an extension of Eq.10 that accounts for a possible rotational misalignment between the 2D layers and describes also the tunneling between inequivalent extrema. The vector is zero only for tunneling between equivalent extrema (i.e. ) and for a perfect rotational alignment (i.e. 0). Considering a case where all extrema are at the point, we have 43, then for the magnitude of is simply given by Feenstra et al. (2012).

One significant difference in Eq.25 compared to Eq.10 is that, in the presence of rotational misalignment, the top layer Bloch function has a different periodicity in the principal coordinate system from the bottom layer . Consequently the integral over the unit cells of the bottom 2D layer is not the same in all unit cells, so that the derivations going from Eq.10 to Eq.15 should be rewritten accounting for a matrix element depending on the unit cell . Such an could be formally included in the calculations by defining a new scattering spectrum that includes not only the inherently random fluctuations of the potential , but also the cell to cell variations of the matrix element . A second important difference of Eq.25 compared to Eq.10 lies in the presence of Q in the exponential term multiplying F.

For the case of tunneling between inequivalent extrema and with a negligible rotational misalignment (i.e. 0), Eq.26 gives QKK and the current can be expressed as in Eq.15 but with the scattering spectrum evaluated at qQ. Since in this case the magnitude of Q is comparable to the size of the Brillouin zone, the tunneling between inequivalent extrema is expected to be substantially suppressed if the correlation length of the scattering spectrum S(q) is much larger than the lattice constant, as it has been assumed in all the derivations.

Quite interestingly, the derivations in this section suggest that a possible rotational misalignment is expected to affect the absolute value of the tunneling current but not to change significantly its dependence on the terminal voltages.

From a technological viewpoint, if the stack of the 2D materials is obtained using a dry transfer method the rotational misalignment appears inevitable Britnell et al. (2013); Britnell et al. (2012b). Experimental results have shown that, when the stack of 2D materials is obtained by growing the one material on top of the other, the top 2D and bottom 2D layer can have a fairly good angular alignment Tiefenbacher et al. (2000); Koma (1999).

### ii.5 An analytical approximation for the tunneling current

The numerical calculations for the tunneling current obtained with the model derived in Secs.II.2 and II.3 will be presented in Sec.III, while in this section we discuss an analytical, approximated expression for the tunneling current which is mainly useful to gain an insight about the main physical and material parameters affecting the current versus voltage characteristic of the Thin-TFET. In order to derive an analytical current expression we start by assuming a parabolic energy relation and write

(27) |

where , are the energy relation respectively in the bottom layer valence band and top layer conduction band and , the corresponding effective masses.

In the analytical derivations we neglect the energy broadening and start from Eq.15, so that the model is essentially valid only in the on-state of the device, that is for .

We now focus on the integral over k and k in Eq.15 and first introduce the polar coordinates kk,, k k,, and then use Eq.27 to convert the integrals over k, k to integrals over respectively E, E, which leads to

where the spectrum (q) is given by Eq.16 and thus depends only on the magnitude of . Assuming , the Dirac function reduces one of the integrals over the energies and sets , furthermore the magnitude of depends only on the angle , so that Eq.II.5 simplifies to

(29) |

In the on-state condition (i.e. for ), the zero Kelvin approximation for the Fermi-Dirac occupation functions f, f can be introduced to further simplify Eq.29 to

(30) |

where max{}, min{} define the tunneling window .

The evaluation of Eq.30 requires to express as a function of the energy inside the tunneling window and of the angle between k and k. By recalling , we can use Eq.27 to write

(31) |

with . When Eq.31 is substituted in the spectrum S(q) the resulting integrals over and in Eq.30 cannot be evaluated analytically. Therefore to proceed further we now examine the maximum value taken by . The value leading to the largest is , and the resulting expression can be further maximized with respect to the energy varying in the tunneling window. The energy leading to maximum is

(32) |

and the corresponding is

(33) |

When neither the top nor the bottom layer are degenerately doped the tunneling window is given by EE and EE, in which case the E defined in Eq.32 belongs to the tunneling window and the maximum value of q is given by Eq.33. If either the top or the bottom layer is degenerately doped the Fermi levels become the edges of the tunneling window and the maximum value of q may be smaller than in Eq.33.

A drastic simplification in the evaluation of Eq.30 is obtained for 1, in which case Eq.16 returns to S(q), so that by substituting S(q) in Eq.29 and then in Eq.15 the expression for the current simplifies to

(34) |

where we recall that max{}, min{} define the tunneling window.

It should be noticed that Eq.34 is consistent with a complete loss of momentum conservation, so that the current is simply proportional to the integral over the tunneling window of the product of the density of states in the two 2D layers. Since for a parabolic effective mass approximation the density of states is energy independent, the current turns out to be simply proportional to the width of the tunneling window. In physical terms, Eq.34 corresponds to a situation where the scattering produces a complete momentum randomization during the tunneling process.

As can be seen, as long as the top layer is not degenerate we have and the tunneling window widens with the increase of the top gate voltage V, hence according to Eq.34 the current is expected to increase linearly with V. However, when the tunneling window increases to such an extent that becomes comparable to or larger than 1, then part of the values in the integration of Eq.30 belong to the tail of the spectrum S(q) defined in Eq.16, and so their contribution to the current becomes progressively vanishing. The corresponding physical picture is that, while the tunneling window increases, the magnitude of the wave-vectors in the two 2D layers also increases, and consequently the scattering can no longer provide momentum randomization for all the possible wave-vectors involved in the tunneling process. Under these circumstances the current is expected to first increase sub-linearly with V and eventually saturate for large enough V values.

## Iii Numerical results for the tunneling current

The 2D materials used for the tunneling current calculations reported in this paper are the hexagonal monolayer and . The band structure for and have been calculated by using a density functional theory (DFT) approach Gong et al. (2013); Liu et al. (2013), showing that these materials have a direct bandgap with the band edges for both the valence and the conduction band residing at the point in the 2D Brillouin zone. Fig.4 shows that in a range of about 0.4 eV from the band edges the DFT results can be fitted fairly well by using an energy relation based on a simple parabolic effective mass approximation (dashed lines). Hence the parabolic effective mass approximation appears adequate for the purposes of this work, which is focussed on a device concept for extremely small supply voltages ( 0.5 V). The values for the effective masses inferred from the fitting of the DFT calculations are tabulated in Tab.1 together with some other material parameters relevant for the tunneling current calculations.

In all current calculations we assume a top gate work function of 4.17 eV (Aluminium) and back gate work function of 5.17 eV (p++ Silicon) and the top and bottom oxide have an effective oxide thickness (EOT) of 1 nm (see Fig.1). The top 2D layer consists of hexagonal monolayer MoS while the bottom 2D layer is hexagonal monolayer WTe. An -type and -type doping density of by impurities and full ionization are assumed respectively in the top and bottom 2D layer and the relative dielectric constant of the interlayer material is set to 4.2 (e.g. boron nitride). The voltage between the drain and the source is set to 0.3 V and the back gate is grounded for all calculations, unless otherwise stated.

In Fig. 5, the results of numerical calculations are shown for the band alignment and the current density versus the top gate voltage . Figure 5(a) shows that the top gate voltage can effectively govern the band alignment in the device and, in particular, the crossing and uncrossing between the conduction band minimum E in the top layer and the valence band maximum E in the bottom layer, which discriminates between the on and off state of the transistor.

The I versus V characteristic in Fig.5(b) can be roughly divided into three different regions: sub-threshold region, linear region and saturation region. The sub-threshold region corresponds to the condition EE (see also Fig.5(a)), where the very steep current dependence on V is illustrated better in Fig.6 and will be discussed below.

In the second region I exhibits an approximately linear dependence on V, in fact the current is roughly proportional to the energy tunneling window, as discussed in Sec.II.5 and predicted by Eq.34, because the tunneling window is small enough that the condition 1 is fulfilled. In this region I is proportional to the long-wavelength part of scattering spectrum (i.e. small values), hence the current increases with L, as expected from Eq.34. The super-linear behavior of I at small V values observed in Fig.5(b) is due to the tail of the Fermi occupation function in the top layer. When V is increased above approximately 0.5V, the current in Fig.5(b) enters the saturation region, where I increasing with V slows down because of the decay of the scattering spectrum S for values larger than (see Eq.16).

In Fig.6 we analyze the I-V curves for different interlayer thicknesses T and broadening energies ; in all cases an average inverse sub-threshold slope (SS) is extracted in the I range from 10 and 1 [A]. Figure 6(a) shows that the tunneling current increases exponentially by decreasing T, and the decay constant 3.8 nm employed in our calculations results in a dependence on T that is roughly consistent with the dependence experimentally reported in graphene based interlayer tunneling devices Britnell et al. (2012a). The threshold voltages are also shifted to lower values by increasing T. It can be seen that the T impact on SS is overall quite modest and for all the T values the simulations indicate a very steep I-V curve in the sub-threshold region ( 20 mV/dec).

Figure 6(b) shows that according to the model employed in our calculations SS is mainly governed by the parameter of the energy broadening (Eq.22). This result is expected, as already mentioned in Sec.II.3, since in our model the energy broadening is the physical factor setting the minimum value for SS and the I versus V approaches a step-like curve when is zero due to the step-like DoS of these 2D semiconductors Agarwal and Yablonovitch (2011). These results suggest that the energy broadening in the 2D materials plays a very critical role in achieving experimentally low SS values in the proposed Thin-TFETs.

## Iv Discussion and conclusions

This paper proposed a new steep slope transistor based on the interlayer tunneling between two 2D semiconductor materials and presented a detailed model to discuss the physical mechanisms governing the device operation and to gain an insight about the tradeoffs implied in the design of the transistor.

The tunnel transistor based on 2D semiconductors has the potential for a very steep subthreshold region and the subthreshold swing is ultimately limited by the energy broadening in the two 2D materials. The energy broadening can have different physical origins such as disorder, charged impurities in the 2D layers or in the surrounding materials Ghazali and Serre (1985) Sarma and Vinter (1982), phonon scattering Bockelmann and Bastard (1990) and microscopic roughness at interfaces Knabchen (1995). In our calculations we accounted for the energy broadening by assuming a simple gaussian energy spectrum with no explicit reference to a specific physical mechanism. However, a more detailed and quantitative description of the energy broadening is instrumental in physical modeling of the device and its design.

Quite interestingly, our analysis suggests that, while a possible rotational misalignment between the two 2D layers can affect the absolute value of the tunneling current, the misalignment is not expected to significantly degrade the steep subthreshold slope, which is the crucial figure of merit for a steep slope transistor.

An optimal operation of the device demands a good electrostatic control of the top gate voltage V on the band alignments in the material stack, as shown for example in Fig.5(a), which may become problematic if the electric filed in the interlayer is effectively screened by the high electron concentration in the top 2D layer. Consequently, since high carrier concentrations in the 2D layers are essential to reduce the layer resistivities, a tradeoff exists between the gate control and layer resistivities; as a result, doping concentrations in these 2D layers are important design parameters in addition to tuning the threshold voltage. In this respect, chemical doping of TMD materials have been recently demonstrated Fang et al. (2013, 2012), however these doping technologies are still far less mature than they are for 3D semiconductors, and improvements in in-situ doping will be very important for optimization of the device performance. Since our model does not include the lateral transport in the 2D materials, an exploration of the above design tradeoffs goes beyond the scope of the present paper and demands the development of more complete transport models.

The transport model proposed in this work does not account for possible traps or defects assisted tunneling, which have been recently recognized as a serious hindrance to the experimental realization of Tunnel-FETs exhibiting a sub-threshold swing better than 60 mV/dec Pala and Esseni (2013); Esseni and Pala (2013). A large density of states in the gap of the 2D materials may even lead to a Fermi level pinning that would drastically degrade the gate control on the band alignment and undermine the overall device operation. In this respect, from a fundamental viewpoint the 2D crystals may offer advantages over their 3D counterparts because they are inherently free of broken/dangling bonds at the interfaces Jena (2013). However, the fabrication technologies for 2D crystals are still in an embryonal stage compared to technologies for conventional semiconductors, hence the control of defects in the 2D materials will be a challenge for the development of the proposed tunneling transistor.

The simulation results reported in this paper indicate that the newly proposed transistor based on interlayer tunneling between two 2D materials has the potential for a very steep turn-on characteristic, because the vertical stack of 2D materials having an energy gap is probably the device structure that allows for the most effective, gate controlled crossing and uncrossing between the edges of the bands involved in the tunneling process. Our modeling approach based on the Bardeen’s transfer Hamiltonian is by no means a complete device model but instead a starting point to gain insight about its working principle and its design. At the present time an experimental demonstration of the device appears of crucial importance, first of all to validate the device concept, and then to help estimate the numerical value of a few parameters in the transport model that can be determined only by comparing to experiments.

Acknowledgments: This work was supported in part by the Center for Low Energy Systems Technology (LEAST), one of six SRC STARnet Centers, sponsored by MARCO and DARPA, by the Air Force Office of Scientific Research (FA9550-12-1-0257), and by a Fulbright Fellowship for D. Esseni. The authors are also grateful for the helpful discussions with Profs. K. J. Cho, R. Feenstra and A. Seabaugh. The authors are especially thankful to Dr. Cheng Gong in Dr. K. J. Cho’s group for providing the calculated band structure data shown in Fig. 4.

## References

- Rabaey et al. (2002) J. M. Rabaey, J. Ammer, T. Karalar, S. Li, B. Otis, M. Sheets, and T. Tuan, in Solid-State Circuits Conference, 2002. Digest of Technical Papers. ISSCC. 2002 IEEE International (IEEE, 2002), vol. 1, pp. 200–201.
- Amirtharajah and Chandrakasan (1998) R. Amirtharajah and A. P. Chandrakasan, Solid-State Circuits, IEEE Journal of 33, 687 (1998).
- Dreslinski et al. (2010) R. G. Dreslinski, M. Wieckowski, D. Blaauw, D. Sylvester, and T. Mudge, Proceedings of the IEEE 98, 253 (2010).
- Group et al. (2011) I. W. Group et al., URL http://www. itrs. net (2011).
- Seabaugh and Zhang (2010) A. C. Seabaugh and Q. Zhang, Proceedings of the IEEE 98, 2095 (2010).
- Zhou et al. (2012) G. Zhou, R. Li, T. Vasen, M. Qi, S. Chae, Y. Lu, Q. Zhang, H. Zhu, J.-M. Kuo, T. Kosel, et al., in Electron Devices Meeting (IEDM), 2012 IEEE International (IEEE, 2012), pp. 32–6.
- Tomioka et al. (2013) K. Tomioka, M. Yoshimura, and T. Fukui, Nano letters (2013).
- Knoll et al. (2013) L. Knoll, Q.-T. Zhao, A. Nichau, S. Trellenkamp, S. Richter, A. Schafer, D. Esseni, L. Selmi, K. K. Bourdelle, and S. Mantl, Electron Device Letters, IEEE 34, 813 (2013).
- Mohata et al. (2011) D. Mohata, R. Bijesh, S. Mujumdar, C. Eaton, R. Engel-Herbert, T. Mayer, V. Narayanan, J. Fastenau, D. Loubychev, A. Liu, et al., in Electron Devices Meeting (IEDM), 2011 IEEE International (IEEE, 2011), pp. 33–5.
- Conzatti et al. (2012) F. Conzatti, M. Pala, D. Esseni, E. Bano, and L. Selmi, Electron Devices, IEEE Transactions on 59, 2085 (2012).
- Pala and Esseni (2013) M. Pala and D. Esseni, Electron Devices, IEEE Transactions on 60, 2795 (2013), ISSN 0018-9383.
- Esseni and Pala (2013) D. Esseni and M. G. Pala, Electron Devices, IEEE Transactions on 60, 2802 (2013).
- Feenstra et al. (2012) R. M. Feenstra, D. Jena, and G. Gu, Journal of Applied Physics 111, 043711 (2012).
- Britnell et al. (2013) L. Britnell, R. Gorbachev, A. Geim, L. Ponomarenko, A. Mishchenko, M. Greenaway, T. Fromhold, K. Novoselov, and L. Eaves, Nature communications 4, 1794 (2013).
- Britnell et al. (2012a) L. Britnell, R. V. Gorbachev, R. Jalil, B. D. Belle, F. Schedin, M. I. Katsnelson, L. Eaves, S. V. Morozov, A. S. Mayorov, N. M. Peres, et al., arXiv preprint arXiv:1202.0735 (2012a).
- Radisavljevic et al. (2011) B. Radisavljevic, A. Radenovic, J. Brivio, V. Giacometti, and A. Kis, Nature nanotechnology 6, 147 (2011).
- Wang et al. (2012) Q. H. Wang, K. Kalantar-Zadeh, A. Kis, J. N. Coleman, and M. S. Strano, Nature nanotechnology 7, 699 (2012).
- Gong et al. (2013) C. Gong, H. Zhang, W. Wang, L. Colombo, R. M. Wallace, and K. Cho, Applied Physics Letters 103, 053513 (2013).
- Jena (2013) D. Jena, Proceedings of the IEEE 101, 1585 (2013), ISSN 0018-9219.
- Mak et al. (2010) K. F. Mak, C. Lee, J. Hone, J. Shan, and T. F. Heinz, Physical Review Letters 105, 136805 (2010).
- Esseni et al. (2011) D. Esseni, P. Palestri, and L. Selmi, Nanoscale MOS transistors: Semi-classical transport and applications (Cambridge University Press, 2011).
- Bardeen (1961) J. Bardeen, Phys. Rev. Letters 6 (1961).
- Harrison (1961) W. A. Harrison, Physical Review 123, 85 (1961).
- Duke (1969) C. B. Duke, Tunneling in solids (Academic Press New York, 1969), vol. 1999.
- Zhao et al. (2013) P. Zhao, R. Feenstra, G. Gu, and D. Jena, Electron Devices, IEEE Transactions on 60, 951 (2013), ISSN 0018-9383.
- Zhu et al. (2011) Z. Zhu, Y. Cheng, and U. Schwingenschlögl, Physical Review B 84, 153402 (2011).
- Ramasubramaniam et al. (2011) A. Ramasubramaniam, D. Naveh, and E. Towe, Physical Review B 84, 205325 (2011).
- Li et al. (2011) Q. Li, E. Hwang, E. Rossi, and S. D. Sarma, Physical review letters 107, 156601 (2011).
- Yan and Fuhrer (2011) J. Yan and M. S. Fuhrer, Physical Review Letters 107, 206601 (2011).
- Yankowitz et al. (2012) M. Yankowitz, J. Xue, D. Cormode, J. D. Sanchez-Yamagishi, K. Watanabe, T. Taniguchi, P. Jarillo-Herrero, P. Jacquod, and B. J. LeRoy, Nature Physics 8, 382 (2012).
- Xue et al. (2011) J. Xue, J. Sanchez-Yamagishi, D. Bulmash, P. Jacquod, A. Deshpande, K. Watanabe, T. Taniguchi, P. Jarillo-Herrero, and B. J. LeRoy, Nature materials 10, 282 (2011).
- Decker et al. (2011) R. Decker, Y. Wang, V. W. Brar, W. Regan, H.-Z. Tsai, Q. Wu, W. Gannett, A. Zettl, and M. F. Crommie, Nano letters 11, 2291 (2011).
- Van Mieghem et al. (1991) P. Van Mieghem, G. Borghs, and R. Mertens, Physical Review B 44, 12822 (1991).
- Urbach (1953) F. Urbach, Physical Review 92, 1324 (1953).
- Cody (1992) G. Cody, Journal of non-crystalline solids 141, 3 (1992).
- Kane (1963) E. O. Kane, Physical Review 131, 79 (1963).
- Sarma and Vinter (1982) S. D. Sarma and B. Vinter, Surface Science 113, 176 (1982).
- Knabchen (1995) A. Knabchen, Journal of Physics: Condensed Matter 7, 5209 (1995).
- Ghazali and Serre (1985) A. Ghazali and J. Serre, Solid-State Electronics 28, 145 (1985).
- Britnell et al. (2012b) L. Britnell, R. Gorbachev, R. Jalil, B. Belle, F. Schedin, A. Mishchenko, T. Georgiou, M. Katsnelson, L. Eaves, S. Morozov, et al., Science 335, 947 (2012b).
- Tiefenbacher et al. (2000) S. Tiefenbacher, C. Pettenkofer, and W. Jaegermann, Surface science 450, 181 (2000).
- Koma (1999) A. Koma, Journal of crystal growth 201, 236 (1999).
- Liu et al. (2013) G.-B. Liu, W.-Y. Shan, Y. Yao, W. Yao, and D. Xiao, arXiv preprint arXiv:1305.6089 (2013).
- Agarwal and Yablonovitch (2011) S. Agarwal and E. Yablonovitch, arXiv preprint arXiv:1109.0096 (2011).
- Bockelmann and Bastard (1990) U. Bockelmann and G. Bastard, Physical Review B 42, 8947 (1990).
- Fang et al. (2013) H. Fang, M. Tosun, G. Seol, T. C. Chang, K. Takei, J. Guo, and A. Javey, Nano letters 13, 1991 (2013).
- Fang et al. (2012) H. Fang, S. Chuang, T. C. Chang, K. Takei, T. Takahashi, and A. Javey, Nano letters 12, 3788 (2012).

Bandgap (eV) | Electron affinity () | Conduction band | Valence band | |

effective mass (m) | effective mass (m) | |||

MoS | 1.8 | 4.30 | 0.378 | 0.461 |

WTe | 0.9 | 3.65 | 0.235 | 0.319 |