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

Abstract

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, W and $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.

1. 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, which is what makes them candidates for applications in photonics and optoelectronics [2,3,4]. The electronic properties of these materials are highly sensitive to 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 specific 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 creates the possibility to manipulate them for future applications in spintronics and valleytronics nanodevices [15,21,22,23,24].

On the other hand, TMDs present a high stretchability. Moreover, external strain can be used to efficiently manipulate their electronic and optical properties [25]. Non-uniform strain profiles can be used to create funnels of excitons, which allows for capturing a broad light spectrum, concentrating carriers in specific 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 first-principles methods is that it provides a simple starting point for the further inclusion of many-body electron–electron interaction, external strains, as well as of the dynamical effects 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 fitting 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.

2. Electronic Band Structure

The crystal structure of MX${}_2$ is schematically shown in Figure 1. A single layer is composed by an inner layer of 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\simeq \sqrt{7/12}a$ and $u\simeq a/2$. The lattice parameters of the bulk compounds corresponding to the more commonly studied TMDs are given in

Table 1. Lattice parameters for the TMDs (transition metal dichalcogenides) considered in the text. a represents the M–M atomic distance, u the internal vertical separation between the M and the X planes, and $c^{\prime }$ the distance between the metal M layers. In bulk systems, $c=2c^{\prime }$ accounts for the z-axis lattice parameter. All values are in Å units.

	a	u	${\mathit{c}}^{\prime }$
MoS${}_2$	$3.160$	$1.586$	$6.140$
MoSe${}_2$	$3.288$	$1.664$	$6.451$
WS${}_2$	$3.153$	$1.571$	$6.160$
WSe${}_2$	$3.260$	$1.657$	$6.422$

[36,37,38]. The in-plane Brillouin zone is a hexagon, and it is shown in Figure 2. It contains the high-symmetry points $\Gamma =(0,0)$, K$=4\mathsf{\pi }/3a(1,0)$ and M$=4\mathsf{\pi }/3a(0,\sqrt{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 single-layer TMDs, we use the Density Functional Theory (DFT) calculations presented by some of the authors in Reference [28], in which the intrinsic spin–orbit interaction term for all atoms is included. Figure 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) (DFT calculations were done using the Siesta code [39,40], with the exchange–correlation potential of Ceperly–Alder [41] as parametrized by Perdew and Zunger [42]. A split-valence double-$\mathsf{\zeta }$ basis set including polarization functions was considered [43]. The energy cutoff and the BZ sampling were chosen to converge the total energy with a value of 300 Ry and $30\times 30\times 1$. The energy cutoff was set to ${10}^{-6}$ eV. Spin-orbit interaction of the different compounds was considered by following the method developed in Reference [44]. The lattice parameters used in the calculation are set to their experimental values [45,46,47] and are given in Table 1). 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 Reference [35], they are made by hybridization of the d orbitals of the transition metal, and the p orbitals of the chalcogen. More specifically, the analysis of the orbital content of the set of bands containing the first four conduction bands and the first seven valence bands, which cover an energy window from −7 to 5 eV around the Fermi level, approximately, reveals that these bands are dominated by the five 4d (5d) orbitals of the metal Mo (W) and the six (three for each layer) 3p (4p) orbitals of the chalcogen S (Se), totaling to up to 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 (Figure 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 specific 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 atoms, respectively. Figure 4 and Figure 5 and Figure B1, Figure B2 and Figure B3 represent these relative orbital weights in detail for the different 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, which leads to a large splitting of the valence band at the K and K’ points of the BZ (see Figure 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,48].

3. 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 [49]. We use the model developed by Cappelluti et al. [35], which contains 11 orbitals per layer. In particular, the model contains the five 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 [50,51]. The corresponding base can be expressed as [35]

$$\left(p_x^t,p_y^t,p_z^t,d_{3z^2-r^2},d_{xz},d_{yz},d_{x^2-y^2},d_{xy},p_x^b,p_y^b,p_z^b,\right)$$

(1)

where the indices t and b label the top and bottom chalcogen planes, respectively. The model is defined by the hopping integrals between the different orbitals, which are described in terms of $\mathsf{\sigma }$, $\mathsf{\pi }$ and $\mathsf{\delta }$ ligands. In the following, we reproduce the most important results, and we refer the reader to Reference [35,48] for details of the model. The Slater–Koster parameters that account for the relevant hopping processes of the model are $V_{pd\mathsf{\sigma }}$ and $V_{pd\mathsf{\pi }}$ for M–X bonds, $V_{dd\mathsf{\sigma }}$, $V_{dd\mathsf{\pi }}$ and $V_{dd\mathsf{\delta }}$ for M–M bonds, and $V_{pp\mathsf{\sigma }}$ and $V_{pp\mathsf{\pi }}$ for X–X bonds. Additional parameters of the theory are the crystal fields ${\Delta }_0$, ${\Delta }_1$, ${\Delta }_2$, ${\Delta }_p$, ${\Delta }_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 the out-of-plane $p_z$ X orbitals.

This model can be simplified by performing a 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 11-band model is decoupled into a $6\times 6$ block with even (odd) symmetry of the $p_x,p_y$ ($p_z$) orbitals with respect to $z\to -z$ inversion, and a $5\times 5$ band block with opposite combination. Since low energy excitations belong exclusively to the first block, the fit 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 for by the reduced Hilbert space:

$$\mathsf{\psi }=\left(d_{3z^2-r^2},d_{x^2-y^2},d_{xy},p_x^S,p_y^S,p_z^A\right),$$

(2)

where the S and A superscripts stand for the symmetric and antisymmetric combinations of the p-orbitals of the X atom, $p_i^S=1/\sqrt{2}(p_i^t+p_i^b)$ and $p_i^A=1/\sqrt{2}(p_i^t-p_i^b)$, with $i=x,y,z$. The tight-binding Hamiltonian defined by base (2), including local spin–orbit coupling, can be expressed in real space as

$$\begin{array}{ccc}\hfill H& =& \sum _{i,\mathsf{\mu }\nu }{\mathsf{ϵ}}_{\mathsf{\mu },\mathsf{\nu }}{c}_{i,\mathsf{\mu }}^{†}{c}_{i,\mathsf{\nu }}+\sum _{ij,\mathsf{\mu }\mathsf{\nu }}[t_{ij,\mathsf{\mu }\mathsf{\nu }}{c}_{i,\mathsf{\mu }}^{†}{c}_{j,\nu }+\mathrm{H}.\mathrm{c}.],\hfill \end{array}$$

(3)

where $c_{i,\mu }^{†}(c_{i,\mathsf{\mu }})$ creates (annihilates) an electron in the unit cell i in the atomic orbital $\mathsf{\mu }=1,\dots ,6$, belonging to the Hilbert space defined by base (2). The Hamiltonian in k-space can be expressed as [30,35,48]:

$$\begin{array}{c}\mathcal{H}=\left(\begin{array}{c}{\mathcal{H}}_{MM}{\mathcal{H}}_{MX}\\ {{\mathcal{H}}_{MX}}^{†}{\mathcal{H}}_{XX}\end{array}\right),\\ {\mathcal{H}}_{MM}={\mathsf{ϵ}}_M+2\sum _{i=1,2,3}{t}_i^{MM}cos\left(\mathbf{k}·{\mathbf{a}}_i\right),\\ {\mathcal{H}}_{XX}={\mathsf{ϵ}}_X+2\sum _{i=1,2,3}{t}_i^{XX}cos\left(\mathbf{k}·{\mathbf{a}}_i\right),\\ {\mathcal{H}}_{MX}=\sum _{i=1,2,3}{t}_i^{MX}{e}^{-i\mathbf{k}·{\mathsf{\delta }}_i},\end{array}$$

(4)

where the vectors (${\mathsf{\delta }}_i$) and (${\mathbf{a}}_i$) are shown in Figure 1b. The analytical expressions for the TB model are given in Appendix A.

4. 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 difficult task. Two main requirements must be satisfied: 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 strains do not affect 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 considers the effect of vacancies, adatoms, etc.

In this work, we have obtained the Slater–Koster parameters for each compound from a minimization procedure that has the possibility to consider a band/momentum resolved weight that allows us to resolve more accurately particular $\mathbf{k}$ regions of selected bands (e.g., edges of the valence and conduction bands)(The calculation of tight-binding band structure and fitting to DFT have been performed using the MathQ package, developed by P. San-Jose. (Source: http://www.icmm.csic.es/sanjose/MathQ/MathQ.html)). Furthermore, we can apply constrictions for the orbital contribution at specific 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 2. Spin–orbit coupling ${\mathsf{\lambda }}_{\mathsf{\alpha }}$ and tight-binding parameters for single-layer $M{X}_2$, where the metal M is Mo or W and X is S or Se. All of the Slater–Koster parameters are in units of eV. Spin-Orbit (SO) coupling parameters are taken from Reference [37].

		MoS${}_2$	MoSe${}_2$	WS${}_2$	WSe${}_2$
SOC	${\mathsf{\lambda }}_M$	0.086	0.089	0.271	0.251
${\mathsf{\lambda }}_X$	0.052	0.256	0.057	0.439
Crystal Fields	${\Delta }_0$	−1.094	−1.144	−1.155	−0.935
${\Delta }_1$	−0.050	−0.250	−0.650	−1.250
${\Delta }_2$	−1.511	−1.488	−2.279	−2.321
${\Delta }_p$	−3.559	−4.931	−3.864	−5.629
${\Delta }_z$	−6.886	−7.503	−7.327	−6.759
M–X	$V_{pd\mathsf{\sigma }}$	3.689	3.728	7.911	5.803
$V_{pd\mathsf{\pi }}$	−1.241	−1.222	−1.220	−1.081
M–M	$V_{dd\mathsf{\sigma }}$	−0.895	−0.823	−1.328	−1.129
$V_{dd\mathsf{\pi }}$	0.252	0.215	0.121	0.094
$V_{dd\delta }$	0.228	0.192	0.442	0.317
X–X	$V_{pp\mathsf{\sigma }}$	1.225	1.256	1.178	1.530
$V_{pp\mathsf{\pi }}$	−0.467	−0.205	−0.273	−0.123

. During the fitting, we have used only the $6\times 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, ${\Delta }_1$, which is the crystal field corresponding to $d_2=d_{xz,yz}$ orbitals, whose contribution is absent in the $6\times 6$ block. What we have done to estimate the value of ${\Delta }_1$ is to impose that the edges of the TB and DFT bands coincide for the lower energy band of the $5\times 5$ block at the K point of the Brillouin zone. The band structure calculated with the $6\times 6$ block of this model is plotted in Figure 3, as compared to DFT calculations. In Figure 4 and Figure 5 and Figure B1, Figure B2 and Figure B3, we compare the orbital contribution of the TB model, using the Slater–Koster parameters of Table 2, 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 presents an acceptable fit 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 3. Comparison of the orbital contribution at band edges at K and Γ points obtained from DFT (Density Functional Theory) and TB (tight-binding) models. K${}^v$ (K${}^c$) refers to the edge of the valence (conduction) band at K point, and ${\Gamma }^v$ refers to the edge of the valence band at the Γ point.

		K${}^{\mathit{v}}$		K${}^{\mathit{c}}$		${\mathbf{\Gamma }}^{\mathit{v}}$
		DFT	TB	DFT	TB	DFT	TB
MoS${}_2$	$d_0$	0.0	0.0	0.82	0.77	0.66	0.96
	$d_2$	0.76	1.0	0.0	0.0	0.0	0.0
	$p_{xy}$	0.20	0.0	0.12	0.23	0.0	0.0
	$p_z$	0.0	0.0	0.0	0.00	0.28	0.04
MoSe${}_2$	$d_0$	0.0	0.0	0.83	0.83	0.71	0.96
	$d_2$	0.78	1.0	0.0	0.0	0.0	0.0
	$p_{xy}$	0.18	0.0	0.10	0.17	0.0	0.0
	$p_z$	0.0	0.0	0.0	0.0	0.23	0.04
WS${}_2$	$d_0$	0.0	0.0	0.80	0.76	0.64	0.98
	$d_2$	0.74	0.94	0.0	0.0	0.0	0.0
	$p_{xy}$	0.21	0.06	0.07	0.24	0.0	0.0
	$p_z$	0.0	0.0	0.0	0.0	0.28	0.02
WSe${}_2$	$d_0$	0.0	0.0	0.82	0.86	0.69	0.99
	$d_2$	0.73	0.95	0.0	0.0	0.0	0.0
	$p_{xy}$	0.20	0.05	0.05	0.14	0.0	0.0
	$p_z$	0.0	0.0	0.0	0.0	0.23	0.01

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 fit to the DFT bands cannot be perfect. More sophisticated methods such 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 [52]. For the case presented here, and due to the automatised fitting procedure, we can conclude that the set of parameters presented in Table 2 must be close to the ideal solution.

5. 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 $\mathsf{\sigma }(\mathsf{\omega })$. For this aim, we use the Kubo formula

$$\begin{array}{c}\mathsf{\sigma }(\mathsf{\omega })=\frac{e^2}{h}\frac{\hslash }{A\mathsf{\omega }}\sum _{mn,\mathbf{k}}(f(E_n)-f(E_m))\times \hfill \\ {\left|〈{\Psi }_n(\mathbf{k}){\left|\widehat{v}\right|\Psi }_m(\mathbf{k})〉\right|}^2\delta [\hslash \mathsf{\omega }-(E_n(\mathbf{k})-E_m(\mathbf{k}))],\hfill \end{array}$$

(5)

where A is the area of the unit cell, ${\Psi }_n(\mathbf{k})$ is the eigenstate of energy $E_n$, $f(E_n)=1/(1+e^{\mathsf{\beta }{E}_n})$ is the Fermi–Dirac distribution function, in terms of the inverse temperature $\mathsf{\beta }=1/k_{B}T$ and considering that the Fermi energy lies in the gap, and $\widehat{v}=(1/\hslash )\partial \widehat{\mathcal{H}}/\partial \mathbf{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 Figure 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 $\mathsf{\sigma }(\mathsf{\omega })$ at low energies (see insets of each panel in Figure 5) is due to the SOC, which leads to two sets of optical transitions in the spectrum. Due to the stronger SOC of heavier W atoms, the effect is especially visible in WS${}_2$ and WSe${}_2$, with plateaus of ∼0.4 eV in $\mathsf{\sigma }(\mathsf{\omega })$, corresponding to the energy splitting of the valence band at the K and K’ points of the Brillouin zone (see Figure 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., Reference [53]). At higher energies, the optical conductivity shows a series of peaks that are associated with optical transitions between flat 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 Figure 3). We notice that disorder (vacancies, adatoms, etc.), not included here, can lead to the creation of midgap states that allow for additional optical transitions in the spectrum [29].

6. Conclusions

In summary, we have generalized the tight-binding model for MoS${}_2$ of Reference [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 2, which have been obtained from a fit 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 different compounds. This approximation can be straightforwardly generalized to multi-layer systems with arbitrary stacking orders, heterostructures made from the stacking of layers of different compounds, twisted multilayers, strained and/or disordered samples, etc.