A generic tight-binding model for monolayer, bilayer and bulk MoS2

Molybdenum disulfide (MoS2) is a layered semiconductor which has become very important recently as an emerging electronic device material. Being an intrinsic semiconductor the two-dimensional MoS2 has major advantages as the channel material in field-effect transistors. In this work we determine the electronic structure of MoS2 with the highly accurate screened hybrid functional within the density functional theory (DFT) including the spin-orbit coupling. Using the DFT electronic structures as target, we have developed a single generic tight-binding (TB) model that accurately produces the electronic structures for three different forms of MoS2 - bulk, bilayer and monolayer. Our TB model is based on the Slater-Koster method with non-orthogonal sp3d5 orbitals, nearest-neighbor interactions and spin-orbit coupling. The TB model is useful for atomistic modeling of quantum transport in MoS2 based electronic devices.

Molybdenum disulfide (MoS 2 ) belongs to a family of layered transition metal dichalcogenides (TMDC) in which the layers are held together by weak van der Waals forces, and it can be exfoliated mechanically to a single layer thickness. In its bulk form MoS 2 is an indirect band gap semiconductor which turns into a direct band gap semiconductor for monolayer structure. 1 This intrinsic semiconducting nature of MoS 2 is a major advantage over graphene (which has no intrinsic band gap) as a two-dimensional (2D) channel material in field-effect transistors (FET). Indeed, FET devices based on MoS 2 monolayer and bilayer have already been fabricated in the experimental labs and demonstrated to have useful device performances. 2,3 Several theoretical studies of MoS 2 FET devices 4-6 based on simplified description of the electronic structures within an effective mass model have also been reported recently.
Due to the complex nature of the layered TMDC materials and the great interests in the electronics community for its applications in emerging devices, a more reliable, accurate, and atomistic treatment of the electronic structures of TMDC is desired. Since ab initio models are computationally expensive and often intractable for realistic device structures having large number of atoms, a widely applied intermediary option is the tight-binding (TB) model. TB models are atomistic full-band yet computationally efficient due to the small number of basis orbitals and the rigid nature of the Hamiltonian. Properly parameterized TB models have been extensively applied to predict reliable results for a diverse range of devices including the resonant tunneling diodes, 7 quantum dots, 8 Si nanowireFET, 9 and compound semiconductor heterostructures. 10 However, for the very exciting 2D MoS 2 and other TMDC materials, accurate TB models for electronic structures has not yet been devela) Electronic mail: fzahid@hku.hk oped. It is the purpose of this work to fill this void, and the resulting TB model should be extremely useful for predicting quantum transport properties of MoS 2 based devices.
In particular, we have developed a generic TB model with non-orthogonal sp 3 d 5 orbitals including the nearestneighbour interactions and spin-orbit coupling, that accurately determines the electronic structures of bulk, monolayer and bilayer MoS 2 . A reliable and accurate target band structure is the primary requirement for a successful tight-binding modeling. In this study, for the calculations of the target band structures we employ density functional theory (DFT) with the screened hybrid functional of Heyd, Scuseria, and Ernzerhof 11 (HSE) that has been shown in the literature [12][13][14] to produce accurate band gaps and reasonable effective masses for a wide range of semiconductors. Our calculated band structures of MoS 2 from the self-consistent DFT-HSE shows excellent agreement with the available experimental data and compare well with other theoretical studies. [15][16][17] Even though an accurate target band structure is obtained by DFT-HSE, it is nontrivial to develop a generic TB model which can accurately capture electronic structures and effective masses of not only bulk MoS 2 but also its doublelayer and single-layer forms. For the layered TMDC structures there is a thickness dependent interlayer interaction that significantly affects the electronic properties, the TB model must capture these effects so as to be applicable to devices made of different layer thicknesses. As shown below, our generic TB model accurately produces the band structures for three different forms of MoS 2 : bulk, monolayer and bilayer. The reported TB model can be easily extended to other MoS 2 structures of higher number of layers as well as other layered TMDC materials such as WS 2 , WSe 2 , and MoSe 2 . The accuracy of our TB model is validated by comparing the band gaps and effective masses with our calculated ab initio results. In our study, all the DFT-HSE calculations are performed using the projector augmented-wave (PAW) pseudopotential plane-wave method as implemented in the VASP software package. 18 A Monkhorst-Pack scheme is adopted for the sampling of the Brillouin zone with kpoint meshes of 7 × 7 × 2 for the bulk and 9 × 9 × 1 for the two-dimensional (2D) structures. An energy cutoff of 280 eV is used in a plane wave basis set. We also included the spin-orbit coupling self-consistently. Due to the presence of van der Waals interaction which determines the interlayer distance in MoS 2 , geometry optimization becomes quite tricky and complicated. Without the van der Waals intereaction the lattice parameter c (vertical length of the unit cell) is overestimated by a large margin. However, it is not yet possible to include van der Waals interaction in the HSE calculations, and it is an area of active research. We therefore optimize the bulk crystal structure using the Perdew-Burke-Ernzerhof (PBE) functional 19 for which the van der Waals interaction can be included explicitly by adopting the DFT-D2 method. 20 We then use this optimized structure for the HSE calculations. The values of our optimized lattice parameters for bulk MoS 2 are: a = 3.179Å and c = 12.729 A with a layer thickness of 3.135Å (S-S vertical distance). Although the value of c parameter improves significantly with the DFT-D2 method, it is still overestimated by around 3.6% compared to the experimental value.
It is a well-known fact that DFT within the local density approximation (LDA) or the generalized gradient approximation (GGA) underestimates band gaps of semiconductors due to factors such as the self interactions. This error can be minimized in the HSE scheme of hybrid functional method 11 by incorporating 25% short-range (SR) exact Hartree-Fock (HF) exchange with the PBE functional. A screening parameter µ defines the range separation, and is empirically set to 0.2Å −1 (HSE06 scheme) for both the HF and PBE parts. In our calculations for MoS 2 we observe that the band gap energies vary significantly with the screening parameter for both the bulk and the 2D structures. We therefore decided to adjust the value of µ to fit the experimental band gap, and an optimize value for µ is found to be 0.4Å −1 which corresponds to an interaction range of π/µ = 7.85Å for the SR nonlocal exchange. Note that a single value of µ is used for all the three structures of MoS 2 : bulk, doublelayer and single-layer, namely the HSE results presented here were all obtained using µ = 0.4Å −1 . Our calculated band structure results for MoS 2 are plotted in Fig. 1. We notice that the monolayer MoS 2 has a direct band gap at the K point. However, with the addition of just one extra layer (i.e. for a bilayer structure) an indirect band gap opens up due to the presence of interlayer interactions. Although both bulk and bilayer MoS 2 have indirect band gaps the conduction band minima (CBM) are located at different points: for bulk it is at Σ point (midpoint of Γ and K points) and for bilayer it is at K point. Similar behavior is observed in a previous theoretical study. 15 For monolayer MoS 2 the splitting of the valence band maximum (VBM) at K point is solely due to the spin-orbit coupling (SOC) whereas a combination of SOC and interlayer interactions is responsible for the VBM splitting in bulk and bilayer MoS 2 . The splitting of CBM due to SOC is minimal (around 5 meV). The calculated band gap energies for different transitions are presented in Table I. Our DFT-HSE results show good agreement with the experimental data. The values of effective masses at different band edges along different directions as presented in Table II also compare well with other available theoretical studies. 15,17 With the calculated DFT-HSE band structures as target, we employed the Nanoskif 22 software package to ob- tain a set of TB parameters for the on-site energies, the Slater-Koster energy integrals, overlap integrals, 23 and the spin-orbit splitting. 24 In our TB model we adopted a non-orthogonal basis set of sp 3 d 5 orbitals considering only the nearest neighbour interactions. The effects of spin-orbit coupling is included through a split-off energy term. A set of 96 parameters are optimized through the automated process built in Nanoskif. 22 The root-meansquare deviation of the fitting is within 25 meV. The optimized TB parameters are listed in Table III. Note that this single set of TB parameters is capable of producing accurate band structures of all three different forms of MoS 2 . The band structures obtained from our TB model are plotted in Fig. 1 along with those obtained from DFT-HSE. The agreement in energies around both VBM and CBM for all three structures of MoS 2 is excellent. Fig. 2 shows the energy contours around the valence band maximum (VBM) and the conduction band minimum (CBM) for the monolayer MoS 2 . Again, we observe excellent agreement in the shape of the energy contours from DFT-HSE and TB fitting. The band gap energies for different transitions and the effective mass values at different symmetry points obtained from the TB model are presented in Table I and Table II respectively. For the band gap energies a fitting accuracy of less than 2.5% is achieved. On the other hand, for the effective masses, in most cases the deviations between the DFT-HSE and the TB values are less than 10% which is quite acceptable. In four cases (excluding the monolayer structure) the deviations are higher with the worst situation at ∼ 27.9%. Given that just a single set of TB parameters is used to produce the band structurues of all three different forms of the material, this level of quantitative consistency is rather satisfactory. In conclusion, we have developed a generic TB model for accurately calculating the band structures of bulk, monolayer and bilayer MoS 2 . Our TB model is based on the Slater-Koster method. For the optimization of the TB parameters, accurate target band structures are obtained using the screened hybrid DFT method. The ac-curacy of our TB model is verified by comparing the band gaps for different transitions and the effective masses at different band edges against the ab initio band structures. One main feature of our TB model is that with only one set of parameters it can reproduce the band structures of MoS 2 of different structural configurations: from bulk TABLE III. Tight-binding parameters for MoS2 using non-orthogonal model with sp 3 d 5 orbitals, nearest-neighbour interactions, and spin-orbit coupling: on-site energies (E), spin-orbit splitting (λ), Slater-Koster energy integrals (E1 for intra-layer and E2 for inter-layer interaction) and overlap integrals (O1 for intra-layer and O2 for inter-layer interaction). The energies are in the unit of eV.

Es
Ep E d λSO to 2D structures. Our tight-binding model can be easily extended to other TMDC materials that show electronic characteristics similar to MoS 2 . The TB model reported in this paper is useful for simulations of quantum transport in nanoelectronic devices based on the TMDC materials.