Electronic Band Structure of Transition Metal Dichalcogenides from Ab Initio and Slater-Koster Tight-Binding Model

Semiconducting transition metal dichalcogenides present a complex electronic band structure with a rich orbital contribution to their valence and conduction bands. The possibility to consider the electronic states from a tight-binding model is highly useful for the calculation of many physical properties, for which first principle calculations are more demanding in computational terms when having a large number of atoms. Here, we present a set of Slater-Koster parameters for a tight-binding model that accurately reproduce the structure and the orbital character of the valence and conduction bands of single layer MX$_2$, where M = Mo,Wand X = S, Se. The fit of the analytical tight-binding Hamiltonian is done based on band structure from ab initio calculations. The model is used to calculate the optical conductivity of the different compounds from the Kubo formula.

Semiconducting transition metal dichalcogenides present a complex electronic band structure with a rich orbital contribution to their valence and conduction bands. The possibility to consider the electronic states from a tight-binding model is highly useful for the calculation of many physical properties, for which rst principle calculations are more demanding in computational terms when having a large number of atoms.
Here we present a set of Slater-Koster parameters for a tight-binding model that accurately reproduce the structure and the orbital character of the valence and conduction bands of single layer MX2, where M = Mo, W and X = S, Se. The t of the analytical tight-binding Hamiltonian is done based on band structure from ab initio calculations. The model is used to calculate the optical conductivity of the di erent compounds from the Kubo formula.

I. INTRODUCTION
Soon after the discovery of graphene by mechanical exfoliation, this technique was applied to the isolation of other families of van der Waals materials. 1 Among them, semiconducting transition metal dichalcogenides (TMD) are of special interest because they have a gap in the optical range of the energy spectrum, what makes them candidates for applications in photonics and optoelectronics. [2][3][4] The electronic properties of these materials are highly sensitive to the external conditions such as strain, pressure or temperature. For instance, a direct-to-indirect gap and even a semiconducting-to-metal transition can be induced under speci c conditions. [5][6][7][8][9][10] They also present a strong spin-orbit coupling (SOC) that, due to the absence of inversion symmetry in single layer samples, lifts the spin degeneracy of the energy bands. 11 In time reversal-symmetric situations, inequivalent valleys have opposite spin splitting, leading to the so called spin-valley coupling, [12][13][14] which has been observed experimentally. [15][16][17][18][19][20] The coupling of the spin and valley degrees of freedom in semiconducting TMDs opens the possibility to manipulate them for future applications in spintronics and valleytronics nano-devices. 15,[21][22][23][24] On the other hand, TMDs present a high stretchability. Moreover, external strain can be used to e ciently manipulate their electronic and optical properties. 25 Non-uniform strain pro les can be used to create funnels of excitons, that allows to capture a broad light spectrum, concentrating carriers in speci c regions of the samples. 5,26 Strain engineering can be also used to exploit the piezoelectric properties of atomically thin layers of TMDs, converting mechanical to electrical energy. 27 The rich orbital structure of the valence and conduction bands of semiconducting TMDs 28 complicates the construction of a tight-binding (TB) model for these systems. Such a TB model must be precise enough as to include all the pertinent orbitals of the relevant bands, but at the same time, simple enough as to be used without too much effort in calculations of optical and transport properties of these materials. The advantage of a tight-binding description with respect to rst-principles methods is that it provides FIG. 1. a) Sketch of the atomic structure of MX2. The bulk compound has a 2H-MX2 structure with two MX2 layers per unit cell, each of them being built up from a trigonal prism coordination unit. The value of the lattice constants for each family is given in Table  I. b) Top view of monolayer M X2 lattice. Green (Blue) circles indicate M (X) atoms. The nearest neighbors (δi) and the next nearest neighbors (ai) vector are shown in the gure. a simple starting point for the further inclusion of manybody electron-electron interaction, external strains, as well as of the dynamical e ects of the electron-lattice interaction. Tight-binding approaches are often more convenient than ab initio methods for investigating systems involving a very large number of atoms, 26 disordered and inhomogeneous samples, 29 strained and/or bent samples, 30,31 materials nanostructured in large scales (nanoribbons, 32,33 ripples 34 ) or in twisted multilayer materials. The aim of the present paper is twofold. Starting from the TB model for MoS 2 developed by Cappelluti et al., 35 we present a more accurate set of Slater-Koster parameters obtained from a more sophisticated tting procedure, and we further generalize it to the other families of semiconducting TMDs, namely WS 2 , MoSe 2 and WSe 2 . Finally, we apply the obtained tight-binding models to calculate the optical conductivity of the four compounds.

II. ELECTRONIC BAND STRUCTURE
The crystal structure of MX 2 is schematically shown in metal M atoms ordered on a triangular lattice, which is sandwiched between two layers of chalcogen X atoms placed on the triangular lattice of alternating hollow sites. We use a notation such that a corresponds to the distance between nearest neighbor in-plane M − M and X − X atoms, b is the nearest neighbor M − X separation and u is the distance between the M and X planes. The MX 2 crystal forms an almost perfect trigonal prism structure with b 7/12a and u a/2. The lattice parameters of the bulk compounds corresponding to the more commonly studied TMDs are given in Table I. [36][37][38] The in-plane Brillouin zone is an hexagon, and it is shown in Fig. 2. It contains the high-symmetry points Γ = (0, 0), K= 4π/3a(1, 0) and M= 4π/3a(0, √ 3/2). The six Q points correspond to the approximate position of a local minimum of the conduction band.
In order to study the electronic band structure of singlelayer TMDs, we use the Density Functional Theory (DFT) calculations presented by some of the authors in Ref. 28, in which the intrinsic spin-orbit interaction term for all atoms is included. Fig. 3 shows the band structures for single-layer M X 2 (black lines) together with the TB bands that will be discussed later (red lines). 39 One of the main characteristics of TMDs is that, contrary to what happens in other 2D crystals like graphene or phosphorene, the valence and conduction bands of M X 2 present a very rich orbital contribution. As explained in detail in Ref. 35, they are made by hybridization of the d orbitals of the transition metal, and the p orbitals of the chalcogen. More speci cally, the analysis of the orbital content of the set of bands containing the rst four conduction bands and the rst seven valence bands, which a u c MoS2 3.160 1.586 6.140 MoSe2 3.288 1.664 6.451 WS2 3.153 1.571 6.160 cover an energy window from -7 to 5 eV around the Fermi level, approximately, reveals that these bands are dominated by the ve 4d (5d) orbitals of the metal Mo (W) and the six (three for each layer) 3p (4p) orbitals of the chalcogen S (Se), summing up to the 93 % of the total orbital weight. 35 Single-layer TMDs are direct gap semiconductors, with the gap located at the two inequivalent K and K' points of the Brillouin zone (Fig. 3). The main orbital character at the edge of the valence band is due to a combination of d x 2 −y 2 and d xy orbitals of the metal M , which hybridize with p x and p y orbitals of the chalcogen X. On the other hand, the edge of the conduction band is formed by d 3z 2 −r 2 orbital of M , plus some contribution of p x and p y orbitals of X. 35 Contrary to single-layer samples, multi-layer compounds are indirect gap semiconductors. The edge of the valence band lies at the Γ point of the BZ, having a major contribution from d 3z 2 −r 2 and p z orbitals of M and X atoms, respectively. The edge of the conduction band in multi-layer samples is placed at the Q point of the BZ. It is important to notice that the Q point is not a high symmetry point of the Brillouin zone, and therefore its exact position depends on the number of layers and on the speci c compound. The orbital character at the Q point originates mainly from the d xy and d x 2 −y 2 orbitals of the metal M , plus p x and p y orbitals of the chalcogen X, plus a non negligible contribution of p z and d 3z 2 −r 2 of X and M FIG. 4. Band structure and orbital character of single-layer MoS2. The thickness of the bands represents the orbital weight, where the d character (d2 = d x 2 −y 2 , dxy and d0 = d 3z 2 −r 2 ) refers to the Mo atom 4d orbitals, while the p character (where pxy = px, py) refers to 3p orbitals of sulfur. Top panels correspond to orbital weight from DFT calculations, whereas bottom panels correspond to orbital weight from TB results. The black lines in the bottom panels are the DFT bands. Notice that spin-orbit coupling is not included in these plots. atoms, respectively. Figs. 4-8 represent these relative orbital weights in detail for the di erent compounds. The extremely rich orbital contribution to the relevant bands that occur in semiconducting TMDs complicates the derivation of a minimal TB model, valid in the whole Brillouin zone. Another important feature of TMDs is that they present a strong SOC, that leads to a large splitting of the valence band at the K and K' points of the BZ (see Fig. 3). The splitting is bigger for W than for Mo compounds, due to the heavier mass of the former. SOC also leads to a splitting of the conduction band at the K point, as well as at the minimum at the Q point. 37,40

III. TIGHT-BINDING MODEL
In this section we consider the electronic band structure of TMDs, in the whole BZ, from a Slater-Koster tight-binding approximation. 41 We use the model developed by Cappelluti et al., 35 which contains 11 orbitals per layer. In particular, the model contains the ve d orbitals of the metal M atom and the six p orbitals of the two chalcogen X atoms in the unit cell. The used scheme has been recently used in other works studying the electronic band structure of TMDs from a tight-binding perspective. 42,43 The corresponding base can be expressed as 35 (1) where the indices t and b label the top and bottom chalcogen planes, respectively. The model is de ned by the hopping integrals between the di erent orbitals, which are described in terms of σ, π and δ ligands. In the following we reproduce the most important results, and we refer the reader to Refs. 35 and 40 for details of the model. The Slater-Koster parameters that account for the relevant hopping processes of the model are V pdσ and V pdπ for M − X bonds, V ddσ , V ddπ and V ddδ for M − M bonds, and V ppσ and V ppπ for X − X bonds. Additional parameters of the theory are the crystal elds ∆ 0 , ∆ 1 , ∆ 2 , ∆ p , ∆ z , describing respectively the atomic level of the l = 0 (d 3z 2 −r 2 ), the l = 1 (d xz , d yz ), the l = 2 (d x 2 −y 2 , d xy ) M orbitals, the in-plane (p x , p y ) X orbitals and of the out-of-plane p z X orbitals.
This model can be simpli ed by performing an unitary transformation that takes the p orbitals of the top and bottom layers of the X atoms into their symmetric and antisymmetric combinations with respect to the z-axis. This way the 11bands model is decoupled into a 6 × 6 block with even (odd) symmetry of the p x , p y (p z ) orbitals with respect to z → −z inversion, and a 5 × 5 bands block with opposite combination. Since low energy excitations belong exclusively to the rst block, the t to DFT that we will present later will be performed within this sector. Therefore the relevant bands above and below the gap are well accounted by the reduced Hilbert space: where the S and A superscripts stand for the symmetric and antisymmetric combinations of the p-orbitals of the X atom, The tight-binding Hamiltonian de ned by the base (2), including local spin-orbit coupling, can be expressed in real space as where c † i,µ (c i,µ ) creates (annihilates) an electron in the unit cell i in the atomic orbital µ = 1, . . . , 6, belonging to the Hilbert space de ned by (2). The Hamiltonian in k-space can be expressed as: 30,35,40 where the vectors (δ i ) and (a i ) are shown in Fig. 1(c). The analytical expressions for the TB model are given in Appendix A.

IV. SLATER-KOSTER PARAMETERS FROM FITTING TO AB INITIO CALCULATIONS
Finding the optimal set of Slater-Koster parameters for the TB model considered here is a di cult task. Two main requirements must be satis ed: a good reproduction of the structure of the relevant electronic bands, and faithful representation of the orbital contribution along such bands. The last condition is especially relevant because, for example, different kinds of strain do not a ect all the hoppings in the same manner. Therefore capturing the proper orbital contribution is essential when using the TB model for calculations of physical properties of strained membranes. The same happens when one consider the e ect of vacancies, adatoms, etc.
In this work we have obtained the Slater-Koster parameters for each compound from a minimization procedure which has the possibility to consider a band/momentum resolved weight, that allows us to resolve more accurately particular k regions of selected bands (e.g. edges of the valence and conduction bands). 44 Furthermore, we can apply constrictions for the orbital contribution at speci c band regions, taking as a reference the information from the DFT  wave-functions. The sets of Slater-Koster parameters that we have obtained for the four compounds are given in Table II. During the tting we have used only the 6 × 6 block of the Hamiltonian because, as explained above, it accounts for the valence and conduction bands. Therefore, we obtain as output all the Slater-Koster parameters but one, ∆ 1 , which is the crystal eld corresponding to d 2 = d xz,yz orbitals, whose contribution is absent in the 6 × 6 block. What we have done to estimate the value of ∆ 1 is to impose that the edges of the TB and DFT bands coincide for the lower  Fig. 3, as compared to DFT calculations. In Figs. 4,6-8 we compare the orbital contribution of the TB model, using the Slater-Koster parameters of Table II, to the corresponding orbital contribution as obtained from DFT. We show the results for the most relevant orbitals (d 0 = d 3z 2 −r 2 , d 2 = d xy , d x 2 −y 2 , p xy = p x , p y and p z ), and we can conclude that the TB model not only present an acceptable t to the band structure, but importantly, the wave-functions also reproduce the DFT orbital contribution at the most important points of the band structure. Table III contains the main orbital contribution of each compound at the most relevant edges of the band structure, namely valence and conduction bands at K point, and valence band at Γ point of the Brillouin zone. We notice that the main restriction of the TB model considered here is that it only includes up to next-nearest-neighbor hopping terms, and this is why the t to the DFT bands cannot be perfect. More sophisticated methods as DFT based tight-binding Hamiltonians represented in the basis of maximally localized Wannier functions can lead to better agreements, at the cost of inclusion of longer range hopping terms. 45 For the case presented here, and due to the automatised tting procedure, we can conclude that the set of parameters presented in Table II must be close to the ideal solution.

V. OPTICAL CONDUCTIVITY
Once we have the tight-binding models for the four compounds, we can use them to calculate physical observables as, for example, the optical conductivity σ(ω). For this aim we use the Kubo formula where A is the area of the unit cell, Ψ n (k) is the eigenstate of energy E n , f (E n ) = 1/(1 + e βEn ) is the Fermi-Dirac distribution function, in terms of the inverse temperature β = 1/k B T and considering that the Fermi energy lies in the gap, andv = (1/h)∂Ĥ/∂k is the velocity operator. The results are shown, for the approximate range of validity of our TB models (∼ 1 eV above and below the band gap), in Fig. 5. We observe that, for all the compounds there is a threshold for the onset of optical transitions that is equal to the gap ∆. The steplike structure of σ(ω) at low energies (see insets of each panel in Fig. 5) is due to the SOC, that leads to two set of optical transitions in the spectrum. Due to the stronger SOC of heavier W atoms, the e ect is specially visible in WS 2 and WSe 2 , with plateaus of ∼ 0.4 eV in σ(ω), corresponding to the energy splitting of the valence band at the K and K' points of the Brillouin zone (see Fig. 3). These transitions lead to the well known A and B absorption peaks observed in photoluminescence. 16 We further notice that our results for the optical conductivity are in good agreement, even for the onset energy, with experimental measurements (see e.g. Ref. 46). At higher energies, the optical conductivity shows a series of peaks that are associated to optical transitions between at bands in the spectrum (van Hove singularities). Such van Hove singularities are clearly evident for the valence and conduction bands at the M point of the Brillouin zone (see Fig. 3). We notice that disorder (vacancies, adatoms, etc.), not included here, can lead to the creation of midgap states that allows for additional optical transitions in the spectrum. 29

VI. SUMMARY
In summary, we have generalized the tight-binding model for MoS 2 of Ref. 35 to the other families of semiconducting TMDs: MoSe 2 , WS 2 and WSe 2 . Our main result is the set of Slater-Koster parameters of Table II, which have been obtained from a t to DFT calculations in which special care was paid to capture the main orbital contribution of the TB bands at the relevant regions of the band structure. The obtained models have been used to calculate the optical conductivity of the di erent compounds. This approximation can be straightforwardly generalized to multi-layer systems with arbitrary stacking orders, heterostructures made from the stacking of layers of di erent compounds, twisted multilayers, strained and/or disordered samples, etc.