Ab-initio simulation of dissipative transport in tunnel devices based on heterostructures of 2D materials

We present a first-principles model to study tunnel transistors based on van der Waals heterojunctions of 2D materials in the presence of dissipative mechanisms due to the electron–phonon interaction. To this purpose, we employed a reduced basis set composed of unit-cell restricted Bloch functions computed with a plane wave ab-initio solver and performed self-consistent quantum transport simulations within the non-equilibrium Green’s functions formalism. Phonon scattering was included with specific self-energies making use of the deformation potential approximation for the electron–phonon coupling. Our simulations identify the van der Waals tunnel FET as a promising option to attain high on-state currents at low supply voltages, but also show a strong impact of the phonon scattering on the transport properties of such device in the sub-threshold regime.


Introduction
Since the discovery of the family of exfoliable 2D materials such as the transition metal dischalcogenides (TMDs), an intensive research has been dedicated to investigate their possible use for nanoelectronics [1]. Thanks to their atomic thickness, they represent the ultimate limit of geometrical scaling and are candidate to replace silicon in future nanodevices. The ability to provide heterostructures with clean interfaces is a key factor to design the next generation of electronic [2] and optical devices [3,4]. Moreover, in the quest of more efficient low-power systems, the possibility to use type II or type III heterostructures opens the way to design tunnel transistors supporting high on-state currents and, thanks to the lows trap densities, with sub-threshold swings (SS) smaller than the thermionic limit of 60 mV/dec at room temperature [5].
Heterostructures of 2D materials can be either lateral or vertical. The former have already been fabricated experimentally [6] and present the advantage to be free from mismatched alignment [7]. Using them may be an option to improve the performance of tunnel devices where on-state currents remain relatively small for most materials because of large band energy gaps and effective masses [8]. This problem can be solved by adopting van der Waals (vdW) heterostructures, where an inversion of the valence and conduction band can be induced by the electric field occurring at high voltages, thus allowing high on-current values [8][9][10]. However, a possible drawback of such a vdW heterojunction (HJ) is that it presents a smaller bandgap even in the off state, which boosts the degradation of the off-current with phonon-assisted tunneling [11].
Due to the possible combination of many 2D materials [12] and the lack of extensive experimental results, it is convenient to rely on simulations in order to identify the best candidates for such a device. State-of-theart ab-initio models like density functional theory (DFT) [13] have been extensively used to study the electronic properties of these materials and recently adopted also in quantum transport simulations of electron devices within the non-equilibrium Green's functions (NEGF) framework [8,10,14]. Even if we suppose the absence of interfacial traps that can degrade the off-state current ( I off ) of 1 3 these devices [15], in order to provide realistic simulations of these devices at room temperature, it is important to include dissipative effects due to the electron-phonon (el-ph) interaction, whose coupling can be also assessed by means of ab-initio approaches such as the density functional perturbation theory (DFPT) [16].
In this work, we propose to study vdW tunnel FETs with an NEGF solver based on DFT-based Hamiltonians, including a simplified description of electron-phonon (el-ph) interaction based on the deformation potential approximation. To perform this study, we have identified the heterojunction composed of the monolayers of 1T-SnS 2 and 1T-HfSe 2 as an optimal candidate as channel material. It is indeed a type-II HJ with a bandgap around 0.24 eV [17], which is small enough to provide the inversion of valence and conduction bands with relatively low gate voltages.
This work is organized as follows: Sect. 2 describes the adopted methodology, and Sec. 3 details the simulation results with an emphasis on the role played by the el-ph interaction in the physics of the vdW TFETs. Finally, Sec. 4 addresses the conclusion.

Methodology
The electronic properties of the materials under investigation were computed with plane wave (PW) DFT calculations by means of the Quantum ESPRESSO suite [13]. Since the hexagonal primitive cell of the monolayers of 1T-SnS 2 and 1T-HfSe 2 is not suited for transport calculations, we have considered orthorhombic unit cells as sketched in Fig. 1. Moreover, we adopted norm-conserving pseudopotentials based on the Perdew-Burke-Ernzerhof (PBE) [18] exchange-correlation functional and considered the non-local van der Waals functional vdW-DF3 [19] for the vdW interactions. To obtain the convergence we used a cutoff energy of 60 Ry and a 20 × 20 × 1 Monkhorst-Pack k-points grid.

Quantum transport simulation
By extracting the DFT Hamiltonian matrices computed by the PW DFT solver, we were able to perform self-consistent quantum transport simulations within the NEGF framework, by using the Green Tea code, as in the model presented in Ref. [14]. Essentially, this approach consists in the following steps. First, the PW Hamiltonian matrix is rewritten into a hybrid basis composed by specific unitary transformation [20], which makes it suitable for transport calculations since in this basis set the Hamiltonian matrix has a block tridiagonal shape, which allows us to use recursive algorithms [21] to efficiently compute the Green's functions of interest. Second, the order of each block matrix is strongly decreased by projecting them into a reduced basis composed of unitcell restricted Bloch functions (URBF) [14]. For each unit cell along the transport direction, the URBFs are in turn the solutions of eigenvalue problem being k = (k x , k yz ) and H 0,0 , H 0,1 the blocks of Hamiltonian matrix in the hybrid basis. The reduced basis set consists in a subset of the {Φ n,k } corresponding to a few k x in the reduced zone and to some tens of bands. For a given lateral wave vector k y , we typically choose the four wave vectors k x = (−0.25, 0, 0.25, 0.5) × 2 ∕a x to define the basis set. Hence, the Hamiltonian matrix in the URBF has a much lower order than the one in the PW basis set, but it provides the same electronic bands, as shown in the example of Fig. 1.
The adoption of this basis set to reduce the matrix order is an alternative to the approach based on the maximally localized Wannier functions (MLWF) [8,[22][23][24]. Even if the reduction in the size problem with respect to the MLWF method is comparable, a possible advantage of the URBF method is to avoid the extraction of the MLWF as a postprocessing of first-principle calculations, which in some cases can be quite delicate and computationally demanding.
Once that Hamiltonian matrix is obtained, the Green's functions matrices in the URBF basis are computed by solving their kinetic equations where G r(a) denote the retarded (advanced) Green's function matrices, G <(>) lesser-than (greater-than) Green's functions, whereas Σ <(>) represent the lesser-than (greater-than) selfenergies which take into account the interaction with the leads and the phonon bath. Here, we assumed ideal ohmic contact at the metal/semiconductor interfaces, therefore neglecting the contact resistance which can be induced by non-ideal interfaces characterized by Fermi-level pinning and metal-induced gap states. Such an approximation can be justified by recent theoretical and experimental works showing that ohmic contacts between metals and 2D materials can be obtained with the adoption of semimetals such as Bi and Sb [25][26][27].
The carrier concentration extracted from the lesser-than and greater-than Green's functions are then self-consistently coupled with the 3D nonlinear Poisson equation to obtain an accurate electrostatic simulation of the system as described in [14]. The self-consistent solution was obtained by setting the Fermi levels of the contacts in accordance with the applied bias and then varying the electrostatic potential in the doped regions until charge neutrality was attained. We assumed Neumann conditions for ungated boundary nodes and Dirichlet conditions for the gated ones.
To include el-ph interaction in our simulations, we opted to use simplified self-energies within the self-consistent Born approximation (SCBA) by assuming the el-ph coupling as local and the phonon bath at equilibrium, meaning that the phonon statistics is expressed in terms of the Bose-Einstein distribution N = [exp (ℏ ∕k B T) − 1] −1 . Further, we also separated the problem for acoustic and optical phonon branches and considered a linear dispersion for the acoustic phonon bands and a dispersion-less band for the optical ones.
Hence, by adopting the deformation potential approximation for the el-ph matrix elements, the elements of the lesser-than self-energy for acoustic phonons read where D ac is the acoustic phonon deformation potential, T is the temperature, is the material density and v s the sound velocity.
Similarly, the elements of the lesser-than self-energy for the optical phonons read where D op,j is the optical phonon deformation potential and j the frequency for the j-th phonon mode. In Eqs. (3) and (4), we have defined the generalized form factor I as where Ψ n,k are the Bloch functions, G the reciprocal space wave vectors and G 0 ≠ 0 accounts for Umklapp scattering processes.

Evaluation of deformation potentials
To compute the self-energy, it is still needed to determine the deformation potentials in Eqs. (3 and 4), which can be assessed by ab-initio methods. The optical phonon frequencies were determined by the calculation of the phonon dispersion within the density functional perturbation theory (DFPT) as implemented in Quantum ESPRESSO [13]. As mentioned before, the electron-phonon coupling was approximated in the form of a deformation potential. For acoustic phonons, it was calculated as in Ref. [28]: where E C∕V is the conduction or valence band energy for electrons and holes, respectively, and is the biaxial strain applied to the unit cell. For optical phonons, the calculation of el-ph matrix elements with the DFPT has shown that the el-ph coupling in materials under investigation is dominated by the polar LO mode. Therefore, we decided to include only this polar phonon mode, whose deformation potential was evaluated from the scattering rate −1 op calculated with the methodology proposed in Refs. [29,30], including the impact of the dielectric environment on the polar phonons. Hence, the optical phonon deformation potential was approximated with the following equation: where is the mass density of the 2D material and m x and m y are the effective masses along x and y directions, respectively. Table 1 shows the deformation potentials computed with this methodology for the monolayers of 1T-SnS 2 and 1T-HfSe 2 .   Figure 2 shows the bandstructures of 1T-SnS 2 , 1T-HfSe 2 , and their vdW heterostructure along the high-symmetry points of the orthorhombic unit cell used in transport calculations. A small hybridization between the orbitals of the two monolayers is reported resulting in a close superposition of the bands of the vdW HJ with those of the two monolayers. This allows us to easily identify the VB of the HJ as mainly composed of HfSe 2 states and the CB of SnS 2 states. Hence, we expect that a sufficiently high vertical electric field can induce an energetic inversion of the minimum of the CB and the maximum of the VB of the HJ at the Γ point. From the PW DFT calculations, we were able to extract the Hamiltonian matrix in the URBF basis, having a reduced order of 32 for the SnS 2 monolayer, 48 for the HfSe 2 monolayer and 74 for the vdW HJ. In transport calculations, we adopted a lateral wave vector discretization of Δk y = 0.1 × 2 ∕a y . The architecture of the vdW tunnel FET is sketched in Fig. 3. The source and drain regions are 8 nm long, the overlap region length can vary from 5 to 20 nm and the gate extension length is of 15 nm. Bottom and top gates cover the overlap region of length L ov and are prolonged for an extension region of length L ext in order to suppress the tunneling between the VB of the bottom layer material and the CB of the top layer one that can occur at low gate voltages [11]. The HfSe 2 monolayer is supposed to be p-type doped with an acceptor concentration of N A = 10 13 cm −2 , while the SnS 2 is undoped in the gated region, but it is n-type doped in the drain region with a donor a concentration of N D = 10 13 cm −2 .

Simulation of vdW TFETs
In order to study the scalability of this device, we simulated the transfer characteristics of the vdW TFET for overlap region lengths varying from 5 to 40 nm in the presence of el-ph scattering, as shown in Fig. 4. Here, the applied source-to-drain voltage is V DS = 0.35 V and the off-state current is chosen as I off = 10 −1 A/ m, whereas the on-state current I on is computed as the current at V GS = V GS,off + V DS , where V GS,off is the gate voltage in the off state. The phonon parameters used for this simulation were those of the SnS 2 given in Tab. 1, this material having the largest scattering rate. In this figure, we can first observe an increase in the current at V GS = 0 V and a significant degradation of the SS when increasing L ov , with SS = 28, 35, 45 and 58 mV/dec for L ov = 5, 10, 20 and 40 nm, respectively. On the other hand, we can notice an increase in the drain current at large gate overdrives, but not as large as expected for a vertical tunnel FET, suggesting that this device behaves more like a point-tunneling (tunneling direction perpendicular to the electric field induced by the gate) than like a line-tunneling FET (tunneling direction parallel to the electric field induced by the gate) [31]. In this regard, we notice that the high local current densities occurring at the overlap region edge could probably induce relevant self-heating effects with a consequent decrease in the I on due to the enhanced phonon scattering. The accurate analysis of such a phenomenon within Fig. 2 Bandstructures of the monolayers of (red symbols) 1T-SnS 2 and (blue symbols) 1T-HfSe 2 , as well as of (black lines) their vdW heterostructure along the high-symmetry points of the orthorhombic unit cell used for transport calculations with cell parameters a x = 3.726 Å, a y = 6.450 Å a full quantum approach [32] is beyond the scope of this paper, but it undoubtedly deserves further studies. A qualitatively similar behavior is reported in Fig. 5, where, in order to study the impact of polar optical phonon inelastic scattering, we considered a scattering rate 10 times larger with respect to that in Fig. 4. In this case, the increase in D op induces a larger degradation of the off-state current with SS values of 45, 61, 78 and 81 mV/dec for L ov = 5, 10, 20 and 40 nm, respectively. On the other hand, the drain current at large gate overdrives is almost independent of the vdW HJ length.
To shed light on the L ov dependence of the off-state current, in Fig. 6, we show the spectral current densities at V GS = 0 V of the devices with (a) L ov = 5 nm, (b) L ov = 10 nm, (c) L ov = 20 nm and (d) L ov = 40 nm. In all Fig. 6a-d, phonon absorption mechanisms are clearly visible in the vdW heterojunction region, where they promote electrons from the valence band of the HfSe 2 to the conduction band of the SnS 2 . This phenomenon is however more pronounced for larger L ov , as a consequence of the longer tunneling region allowing the occurrence of more phonon absorption events.
Similarly, in Fig. 7, we show the spectral current along the transport direction at V GS = 0.6 V for the devices with L ov varying from 5 to 40 nm. In this figure, we can observe that the tunneling current remains almost constant along the overlap region and then experiences a crunching at the edge between the vdW HJ and the SnS 2 , whereas phonon emission events induce electrons to dissipate their energy only in the drain region. The inter-layer current concentration at the channel edge can also be deduced by the shape of the CB and VB, which indicates an accumulated charge concentration induced by this phenomenon. As seen also for the device simulated in Ref. [33], this suggests that the vertical tunneling between the two layers is mainly localized at the edge of the vdW region, such a behavior being the responsible for the weak dependence of the on-state current on the overlap region length.
The analysis of the transfer characteristics in Figs. 4 and 5 suggests that a sub-thermionic behavior can be attained only at relatively short lengths of the vdW HJ and that no significant improvement of the on-state current can be expected by enlarging the overlap region between the two 2D materials.

Conclusion
We have presented a full ab-initio simulation of a vertical tunnel FET based on 2D materials in the presence of inelastic el-ph scattering and studied its dependence on the length of the vdW HJ. The vertical tunnel FET based on the van der Waals heterostructure of 1T-HfSe 2 and 1T-SnS 2 has been identified as a promising option to attain high on-state currents of several hundreds of A/ m for V DS = 0.35 V, thanks to its ability to achieve inversion of valence and conduction bands at high gate voltages. However, this device also proved to be very sensitive to inelastic phononabsorption mechanisms, thus demanding a special care in minimizing of the impact of polar optical modes by means of a suited electrostatic design in order to avoid the loss of the sub-thermionic behavior in the sub-threshold regime.

Data availability
The data of this study are available upon reasonable request. Fig. 7 (colormap) Spectral current density along the transport direction together with the VB and CB profiles at V GS = 0.6 V and V DS = 0.35 V for the vdW TFETs with (a) L ov = 5 nm, (b) L ov = 10 nm, (c) L ov = 20 nm and (d) L ov = 40 nm. Simulation parameters: D ac = 5.47 eV, −1 op = 6.7 ×10 12 s −1 , T = 300 K. The color scale is also normalized with respect to the maximum spectral current density value and ranges from 0 (white) to 1 (red)