Comparative Analysis of Tight-Binding models for Transition Metal Dichalcogenides

We provide a comprehensive analysis of the prominent tight-binding (TB) models for transition metal dichalcogenides (TMDs) available in the literature. We inspect the construction of these TB models, discuss their parameterization used and conduct a thorough comparison of their effectiveness in capturing important electronic properties. Based on these insights, we propose a novel TB model for TMDs designed for enhanced computational efficiency. Utilizing $MoS_2$ as a representative case, we explain why specific models offer a more accurate description. Our primary aim is to assist researchers in choosing the most appropriate TB model for their calculations on TMDs.


Introduction
Monolayers of transition metal dichalcogenides (TMDs) can be regarded as the semiconductor equivalent to graphene [1,2].They are two dimensional (2D) materials with a hexagonal symmetry, similar to graphene.In the unit cell they have one metal (M ) atom and two chalcogen atoms (X ), written as M X 2 with M = (M o, W ) and X = (S, Se, Te).The TMDs have a direct band gap around the two inequivalent K and K ′ points at the corners of the Brillouin zone [3,4].This gap is in the range of visible light.The spin and valley degrees of freedom are coupled what gives rise to interesting phenomena such as spin and valley dependent optical transitions [5].This makes TMDs interesting candidates for spintronic, electronic, optoelectronic and valleytronic applications.
More insights can be gained by studying these phenomena using theoretical methods.Density functional theory (DFT) gives an ab initio description of these materials, but can only be used for smaller systems.Alternatively, simple two-band effective models are constructed using symmetry considerations.These accurately describe the behaviour of massive Dirac fermions seen at the band-edge [6,7].However, these simple models are only valid close to the band-edge.To get a better understanding of the electronic structure of the material, different tight-binding (TB) models were proposed in the literature to give an improved description of the electronic structure.
Compared with graphene, the TB models for TMDs are more involved.In graphene, a simple model describing the p z -orbitals with the out-of-plane π-bond already gives reasonable results, even with just nearest neighbour (NN) hoppings [8].In TMDs on the other hand, the conduction and valence band have rich orbital contributions due to the hybridization of the d-orbitals from the metal atom with the two p-orbitals of the chalcogen atoms.This makes it harder to construct a basic model that can accurately describe the electronic structure of the system, which explains the large number of proposed TB models for TMDs.
(a) This work presents an overview of the relevant TB models for TMDs in literature and proposes a new TB model.The different steps in the construction of the Hamiltonian for the TB model are discussed and how these are connected to symmetries in the TMDs.The TB models are divided in two categories: TB models that take the Slater Koster (SK) [9] two center approximations as basis and ones that take the symmetry group (SG) as a starting point when constructing the Hamiltonian.Furthermore, these TB models are divided according to the orbitals they include and the amount of hoppings they take into account.We give a clear comparison between all these TB models to indicate their strengths and weaknesses, and correspondingly their possible applications.Based on this, we propose a new SG TB model which has a good balance between accuracy and computational efficiency.
The TB models from literature were constructed with various objectives, such as getting a good approximation around certain points in the Brillouin zone, a better fit for the highest valence and lowest conduction bands or for multiple bands.To better compare these models, we refit the parameters to DFT results for MoS 2 .In this work the models are describing the electronic properties of MoS 2 , other materials follow the same trends due to their similar structures.
The uniform construction of the Hamiltonian allows for an elegant implementation of all the discussed TB models.They are all implemented within a new Python package TMDybinding [10] that builds upon the Pybinding software [11].This package simplifies the construction of the Hamiltonian, independent from the TB model that is used in the calculation.
The paper is further organized as follows.In section 2 we present the construction of the Hamiltonian for the TB models considering the relevant orbitals, symmetries, hopping and spin-orbit coupling (SOC) terms.The construction procedures and properties of these TB models are discussed and the new SG model is introduced in this section.We then explain the methods used to fit the parameters in the TB models to the DFT results in section 3 and we compare these fitted models against the DFT results considering different properties: the DOS and its orbital contribution, the band structure, the effect of SOC and the Berry phase.Finally, we give a summary of the models and identify which model to use for specific situations and the possible applications in large-scale systems in section 4.

Space group and symmetries
A single T M D − H layer has a hexagonal structure as shown in figure 1.The unit cell consists of a metal atom (M ) and two out-of-plane chalcogen atoms (X ).For MoS 2 , the metal and chalcogen atoms are M o and S respectively.The in-plane distance from the metal atom to the chalcogen atom is a M X , the distance between bottom (X b ) and top (X t ) the chalcogen atoms is d X −X and the angle between the plane of metal atoms and the line connecting the metal and chalcogen atom is θ .The lattice vectors are given by a 1 = a x and a 2 = a − The outer orbitals of the TMD are the dand p-orbitals for the metal and chalcogen atom respectively, which makes a total of eleven orbitals where p t and p b denote the orbitals from the top and bottom chalcogen respectively.The T M D − H structure belongs to the D 3h group and has a mirror symmetry in the x yplane σh , three mirror symmetries σv along the armchair direction, two threefold rotation axes Ĉ3 through each atomic position, three twofold rotation axes Ĉ′ 2 along the armchair direction and two Ŝ3 -symmetries.The irreducible representations for the orbitals in the Mulliken notation is Under the σh -symmetry, the p-orbitals transform into each other.These p-orbitals are combined to give new orbitals that are even (e) or odd (o) under this symmetry, given by [12][13][14] p e x/ y = With these new p-orbitals, the atomic basis is divided into four parts: the parts that are even and odd under σh and the parts given by the the metal (M ) and chalcogen (X ) atoms, It is assumed that this atomic basis is orthonormal, there is no overlap between the different orbitals.
The hexagonal structure of the TMD with the metal atoms at the blue larger circles and the chalcogen atoms at the smaller yellow circles.The unit cell is indicated in the upper left corner in red.The in-plane hopping vectors r δ l are grouped per n-th NN from left to right.For each of these groups, the hopping matrix T µν,ρ 1,δ associated with r δ 1 is given in appendix B. The hoppings 4a and 4b can be considered as one hopping along the armchair direction, that is rotated along small angles.

Tight-binding Hamiltonian
Using the atomic basis from (3), we write down a general description for the tight-binding Hamiltonian as with and ρ = {e, o}, the Hamiltonian is block diagonal in the even and odd parts under σh .The onsite and hopping energies are indicated as ϵ µ,ρ and T µν,ρ , with {µ, ν} = {M , X }.

Hopping energies and rotations
The hopping energies T µν,ρ connect the different orbitals.For the TB models with multiple orbitals considered here, these hopping energies are represented by matrices.The values of these hopping matrices are different for the different TB models.
In general, they consist of several different n-th nearest neighbour (NN) hopping matrices between the metals and chalcogen atoms.Written in reciprocal space, these matrices also have an additional phase dependence on the distance between the orbitals, with l the direction of this hopping and δ denoting the n-th NN-hopping as given in figure 2, δ = 1, 3, 4a, 4b , when µν = X M , 2µ, 5µ, 6µ , when µ = ν .
Depending on the TB model, a set of different hoppings δ is taken, like δ ∈ {1, 2M , 2X }.The furthest hopping from this set is δ ma x .
The three matrices associated with l ∈ {1, 2, 3} per hopping δ are related by the Ĉ3symmetry.Only one direction for each hopping δ from ν to µ must be specified, the other two are obtained from l = 1 as with the rotation matrices given by with γ = 120 • .The matrix R M ,e transforms the orbitals over double the rotation angle γ because of the |m l | = 2 azimuthal quantum number for the d x 2 − y 2 , d x y -orbitals.For example, the even (ρ = e) first nearest neighbour hoppings (δ = 1) from the metal atom (ν = M ) to the chalcogen atoms (µ = X ) are given by the three vectors r 1 l , where the three matrices given by l ∈ {1, 2, 3} are related by

Hopping energies and symmetries
The exact formulation of the hopping matrices depends on the orbitals that are included in the TB model φ and the set of nearest neighbour hoppings δ.All the symmetries present in the symmetry group for the TMD must be considered by the TB model and its hopping matrices.The σh -symmetry is already satisfied with the even (e) and odd (o) formulation of the orbitals and the Ĉ3 -symmetry by rotating the hopping matrices with the expressions in equation (8).
The other important symmetry operation is the σv symmetry.This last symmetry will further reduce the amount of parameters to a minimal set, giving the final shape of the hopping matrices.Here, we will consider the mirror plane through the yz-axes ( σ1 v ).The two other mirror planes can be considered by a rotation as given by the matrices in equation (9).The basis φ splits in even (g) and odd (u) parts under σ1 v , given by the orbitals When the Hamiltonian is written in this basis, it has parts that are even and odd under σ1 v .Here, we will only write down the hopping matrix for a general hopping between the even (g) and odd (u) sectors, with T g g and T uu denoting the hoppings between orbitals that are even and odd with σ1 v respectively, T gu and T ug denote the hoppings that take a hopping from an even to an odd sector and vice versa.
The hoppings along the armchair directions (δ ∈ {1, 3, 4a, 4b, 5M , 5X }) are parallel to the mirror plane of the σ1 v -symmetry.Under this σ1 v symmetry transformation, the matrix is invariant the sections that connect the even and odd parts will disappear for hoppings along the armchair direction.
For the hoppings along the zigzag directions (δ ∈ {2M , 2X , 6M , 6X }) the mirror plane will be perpendicular to the hopping.The symmetry thus flips the hopping from going to the right to going to the left giving the Hermitian conjugate of the hopping, so that T g g = T † g g , T uu = T † uu and T gu = −T † ug .The remaining Ĉ′ 2 -and Ŝ3 -symmetries are a combination of σh -, σv -and/or Ĉ3 -symmetries that are already satisfied by the hopping matrices.For example, Ĉ′ 2 can be obtained by combining the two mirror symmetries σh and σv .Therefore, a TB model constructed based on the rules related to the σh (even/odd), Ĉ3 (rotations) and σv (shape of matrix) symmetries agrees with all the symmetries in the D 3h symmetry group of the TMD.
The hopping matrices will have different shapes depending on the hopping δ considered and the orbitals φ that are included in the TB model.The values for the parameters in these hopping matrices depend on the used TB model.We call TB models that use all the parameters allowed by the symmetry group (SG), the SG TB models.The matrices ϵ µ,ρ and T µν,ρ 1,δ for the onsite and hopping energies respectively, are given in appendices A and B.

Broken rotation symmetry
When constructing the Hamiltonian, care has to be taken to determine if the hopping is taken with the correct angle and direction.When the hopping is calculated along the wrong angle, say a rotation of 180 • (−r δ l ) from the correct hopping (r δ l ), this corresponds to applying a Ĉ2 -symmetry.This Ĉ2 rotation can be expressed as a combination of σv and σd .The σdsymmetry is a mirror symmetry along the zigzag direction.This σd -symmetry is not a member of the symmetry group for the TMDs as the positions of the metal and chalcogen atoms are not interchangeable.The Ĉ2 is thus not in the symmetry group.
It is possible to transform the TB model to a different basis, connected with a σd operation.Some models in literature [12,15] use this different basis, where the position of the metal and chalcogen atoms are swapped.The parameters given in appendix B and the TMDybinding package always use the same convention for the position of the metal and chalcogen atoms.

Slater-Koster parameterization
The hopping matrices formed by the SG approach still contain a considerable amount of parameters.The two-center integral approximation from Slater and Koster (SK) [9] gives a simplified generalized description for these parameters.Here, only hoppings between pure, hydrogen-like, non-hybridized por d-orbitals and bond types between the orbitals are considered.There are three different types of bonds for the TMDs considered here, σor π-bonds between pand/or d-orbitals and a δ-bond between the d-orbitals.For example, a hopping matrix between the even d-orbitals at different sites can have up to 6 different parameters according to the SG approach.However, with the SK approach, this reduces to only 3 parameters, V d dσ , V d dπ and V d dδ .These parameters are the same for hoppings r δ l with the same radial distance |r δ l |.In the SK approach, the hopping matrix between a metal atom and a chalcogen atom has a contribution for the hopping from the even p e x to the even d x y given by which is the element E x,x y from table 2 in Slater and Koster [9] with l, m and n the components of the unit vector from the p e x -to This approach can be used for all the parameters in the hopping matrices.
The SK approach always produces hopping parameters that are allowed by the symmetry group of a specific material, but does not take into account every broken symmetry in the material.For example, it does not consider the T gu hopping elements from equation (14).Due to the hybridization in the TMDs, the orbitals deviate from the simple, hydrogen-like orbitals.These hybridized orbitals do break certain symmetries, allowing for more parameters in the full SG hopping matrices that are not included in the SK approach [16].
In the SK term from equation ( 15) the parameters are already divided in the even and odd terms under σh .However, if the terms are constructed from the original top-and bottomp-orbitals, there will be an additional cross hopping term between these top-and bottom-porbitals for hoppings δ ∈ {2X , 5X , 6X }.This cross terms leads to additional non-zero elements, also in the T gu -part, that allow for a better description of the TB model.
In appendix D the most general parameterization of the SK approach is noted down, including the splitting in even and odd parameters and the cross-hopping term between the topand bottom-p-orbitals for the hoppings δ = {2X , 5X , 6X }.

Spin-orbit coupling
The spin-orbit coupling (SOC) in TMDs leads to strong modifications in the band structure [7].
To include this effect, the Hamiltonian is expanded in a spin-up (↑) and spin-down (↓) part, The spin-orbit coupling is included in the models by a L • S-term Table 1: Properties of the TB models.With the exception of TB models that include the s-orbital [17,18] and overlap matrices [18,19], all relevant models in literature can be reduced to a model in this table.The dand p-orbitals determine if a spin flip can happen.The full p-orbitals from the X atoms must be present to allow the model to be transformed into one with separate orbitals for the p-orbitals at the X -top and X -bottom atoms.The third NN hopping for Fang is between the even orbitals.(△) The original work [12] has one parameter less than a follow-up work [15].

Model
The summation goes over the different types of atoms µ, with matrices that connect both the spin-up and spin-down components.The exact formulation of the SOC matrices in the expanded orbital basis of (3) is given in appendix C. The main contribution to the SOC is given by the L µ z S µ z -term.When the other terms are neglected, the additional SOC hopping terms only couple bands with the same parity ρ under σh .The Hamiltonian stays block-diagonal.The broken inversion symmetry leads to a spin splitting at K and K ′ .The time reversal symmetry is preserved, so that spin-up states at valley K are degenerate with spin-down states at valley K ′ and vice versa.The same splitting occurs at the Q and Q ′ -points.
The other spin terms give rise to a spin-flip.This spin-flip coincides with a change in |m l | azimuthal quantum number.This SOC breaks the parity ρ, resulting in a denser Hamiltonian.The SOC is best described in the chiral basis, given by This chiral basis is a linear combination of orbitals in the Cartesian basis with the same azimuthal |m l | and angular l quantum number.The orbitals in this basis have a well defined m l -quantum number, which simplifies the description of the L µ ± ladder operator.
The SG TB models from Liu et al [21] are the simplest TB models and only include the even d-orbitals, φ M ,e , resulting in 3 band TB models.These d-orbitals are effective orbitals originating from the hybridization between the dand p-orbitals [30].Though they are categorised with δ ma x = 2 (Liu 2 ) or δ ma x = 6 (Liu 6 ), the models only have metal atoms, such that they are effective TB models with only δ max = 1 or δ max = 3 as maximal nearest neighbour hopping lengths in practice.An extension to the Liu 6 model was given by Wu et al [22] where they also included the odd d-orbitals giving a total of 5 bands, with the basis given by φ M ,e/o .
Fang et al [12,15] reported a SG model with all the orbitals φ M /X ,e/o up to δ max = 3.Compared to Liu 6 /Wu, it also includes the deeper energy-levels in the valence band.
The models from Rostami et al [13], Cappelluti et al [20] and Dias et al [14] are all SK models.These models neglect the cross-hopping between p t -and p b -orbitals on different sites.Dias et al splits the model in even and odd parameters.

A new TB model
We present a novel SG TB model utilizing even orbitals as the basis (φ M /X ,e ) and incorporating hoppings up to second nearest neighbours.For the TMDs investigated in this work, the energy bands around the band gap are all given by contributions from the even orbitals.By only including these even orbitals, the model reduces the number of bands from 11 to 6 while preserving the important features around the band gap.The limited number of hoppings in the model (δ ma x = 2) results in fewer off-diagonal elements yielding an overall sparser matrix which is easier to solve.The chosen combination of hoppings and orbitals makes this model ideally suited for the study of larger systems.The inclusion of the chalcogen atoms in the set of even orbitals ensures a good description for calculations that need these atoms, as for example required in the study of nanoribbons.
However, the models from Fang or Liu 6 will provide a more accurate description, but they are computationally more demanding or do not include the chalcogen atoms, respectively.The model presented here cannot account for terms introduced by out-of-plane deformations, impurities or spin flips as these will mix the even and odd orbitals in the Hamiltonian which the new model cannot describe.The decomposition in top-and bottom p-orbitals is not straightforward, making it less suited for calculations on heterostructures.
In order to provide a clear comparison, we summarize the properties of our new model alongside existing models from the literature in table 1.The parameters for this new model and the existing models are given in appendices B and D. All the models are implemented in TMDybinding, as detailed in appendix H.

Fitting the parameters
In the literature, each model was fitted or constructed with a different objective.This makes a direct comparison difficult.Therefore, we refitted all the models on results obtained from a DFT-calculation.
The energy bands from the DFT calculation were disentangled, i.e. the overlap of the wavefunction at neighbouring points along the path in the reciprocal space 〈u m 〉 was used to differentiate between a crossing and an avoided crossing.We used the function linear sum assignment from Scipy [31] to obtain the ordering of indices n, m where the overlap of the wavefunction was the largest.Then, we sorted these disentangled bands according to energy at the K-point.We neglect the energy bands from DFT that are not taken into account by the Table 2: Ranking of the weights in the cost function from large to small contribution to fit the parameters in the TB models to the DFT results.The bang gap (BG) energy levels are the highest and lowest energy levels in the valence and conduction band respectively.The high-symmetry path is the path that connects the Γ , M , K and Q points.
i Property 1 Effective mass at the K point 2 Effective mass at the Γ point 3 BG energy levels at the Γ , M , K and Q points 4 BG energy levels along the high-symmetry path 5 All energy levels along the high-symmetry path TB models.The bands from the fitted TB model are also disentangled and sorted in the same manner.
The fitting procedure uses a cost function where the various objectives have different weights.The weights were selected to enhance the precision of fitting specific properties within the TB model, a summary is given in table 2. Each objective has a cost function that is based on the sum of the squared error between the function that is fitted and the objective function derived from the DFT calculation, resulting in the total cost function: with C the cost, w i the different weights per objective i and the values for the objective functions for objective i and the different contributions to that objective j, such as the values for the energy at different levels or points in the BZ.This cost function is minimized using the minimize function of SciPy [31] for starting configurations given by the original parameters from the literature.
We found that the SG TB models always converged to a reasonable minimum, where a small change in weight for a specific objective leads to a small change in the parameters.However, we never obtained a parameterization for the TB models using the SK approach that fitted well to the DFT results like we found in the case of the SG approach.The SK models also had difficulties in obtaining the correct orbital character for the energy bands.

DFT calculations
We use the results from DFT-PBE [32] calculations in the GGA-approach as a reference for refitting the models and comparing the results obtained from TB calculations.Although more precise methods are available, the aim is to compare the models in the ability to capture the precise details from a given structure.
The DFT calculations were preformed with Abinit [33] using pseudopotentials obtained from PsueodoDojo [34].Norm-conserving full-relativistic [33] pseudopotentials were used for the calculations with SOC and PAW [34] pseudopotentials for the spin degenerate calculations.The plane-wave energy cutoff was 40H a on a 10 × 10 × 1 Γ -centered Monkhorst-Pack grid of k-points.The monolayers were separated with a 25Å vacuum region.The vacuum level was chosen as reference for the energy.
Abipy [35] was used to analyze the DFT wavefunction to calculate the Berry phase, spinpolarisation, band disentangling and projected DOS.

TB calculations
All the models are implemented in the Python software library Pybinding [11].We created an extension, TMDybinding, that allows the user to easily select which model they want to use and preform the calculations in Pybinding.The installation procedure and an example are included in appendix H.

Band structure
The energy band structures for the different TB models with a comparison to the DFT results are presented in figure 3.Although features at the Γ and K points in the conduction and valence bands are reasonably well described by most models, only the model from Liu 6 seems to give a good description around the Q-point.This model is the only one with a δ = 6 nearest neighbour hopping.This hopping seems to give an important contribution to the shift of the bands at the Q-point as can be seen in appendix G.A long range TB model is thus needed to accurately fit the bandstructure at the Q-point.
Though the deeper energy bands, far from the band gap are not important in most calculations, they still play a role in determining the orbital contribution of all the other bands.As can be seen in the SK TB models (Rostami, Cappelluti and Dias), these deeper energy bands often have wrong energy values compared to the DFT results.By placing wrong energy values at these deeper energies for the TB model, the model can obtain a better fit for the higher energies.The inexact energy value allows the model to get a correct orbital contribution for bands close to the band gap.

DOS
As can be seen in figure 4 from the DFT result, the orbital contribution around the band gap is mainly given by the d-orbitals, φ M ,ρ .All TB models have a higher contribution from the d-orbitals in the conduction band and the higher valence bands, like as seen in the DFT results.Though the shape of the DOS varies between the models, they agree reasonably well around the valence and conduction band edges.
In our novel TB model and the TB models from Liu 2 , Liu 6 and Wu there is a gap in the valence band that mainly originates from the absence of the odd-bands in these models.The DOS for the conduction band from Wu matches the DFT results well, even without the porbitals.1, except for the model from Rostami.The DFT results also include the sorbitals.1, except for the model from Rostami.Between the M and Γ points, the bands are degenerate due to the time reversal symmetry.
In the DOS obtained from the DFT calculations, there are small contributions from the s-orbitals.These orbitals are not described by the TB models, they are relatively small and far from the valence and conduction band edges.

SOC
In figure 5 we give the band structure for the TB models where the SOC contribution is included.The orbital contribution determines the splitting of the bands at the K-point.It is clear that the inclusion of the small L • S-term gives good predictions for the splitting and avoided crossings in the band structure.As mentioned in table 1, the spin-flip term is only possible in models that have the orbitals with consecutive |m l | azimuthal quantum numbers.At the Q-point, the local minimum splits, similar to the K-point.
The TB models from Liu/Wu describe an effective d-orbital, the specific values for the SOC are thus expected to be slightly different than the one from DFT.The metal and chalcogen atoms have different contributions to the λ µ -parameter in equation (18), where as the effective

Berry phase
The Berry phase is a topological quantity that gives more insight in the accuracy of the wavefunctions produced by the TB models.The Berry phase calculation implemented in Pybinding follows the approach introduced in PythTB [36,37].
The Berry phase obtained from the DFT results looks similar to the one obtained from Fang and Liu 6 /Wu.The other TB models have local minima/maxima or even a sign change around the Q-point that is not present in the DFT results.It seems that the inclusion of the third nearest neighbour hopping improves the results for the Berry phase.On the other hand, the Berry phase in the SK models is less accurate compared to the DFT results than the SG models.

Discussion
In general, the TB models based on the SK approach give less accurate results compared to the results obtained from DFT than the SG TB models.It appears that the SK approach does not have enough parameters to capture the electronic properties of the TMDs.For the SKmodels studied here, only two parameters contribute to the hybridization between the pand d-orbitals, V pdσ and V pdπ .The limited set of parameters makes these SK TB models sensitive to small changes in these parameters.This makes it hard to find a parameterization that reproduces the specific properties from the DFT results, such as getting a good fit for energy levels around the band gap and the right orbital contributions for these bands.
The even and odd parameterization in Dias creates duplicates for the even and odd sectors.This allows the parameters to describe the finer details in the TB model for the even and odd parts separately.The longer distance hopping in this TB model makes the Hamiltonian overall denser, which makes it unsuitable for a large scale calculation.
Even when all the different kind of hoppings up to δ max = 6 nearest neighbours are taken into account, we could never obtain a good fit for SK models.If the same long-range model is fitted with a SG TB model, an almost perfect fit with the DFT results is obtained as given in appendix F.
The less accurate results obtained from the SK approach likely originate from the hybridization between the pand d-orbitals in the valence band.The SK orbitals are expected to be hydrogen-like atomic orbitals.Instead, SG orbitals can describe the effective orbitals from the hybridization.
This discrepancy between the SK and SG TB models is illustrated below with an example.In the Liu 2 SG TB model, only the even d-orbitals (φ M ,e ) are taken into account with only the effective first nearest neighbour hoppings, δ = 2M .In the SK model from Rostami et al , this hopping δ = 2M is also present.Turning off all the interactions between the p-orbitals, it results in the same 3-band model.
In appendix E, the Hamiltonian at the Γ and K points is described, and the eigenvalues for the TB models from Liu are given.When the SK parameters are included in this description, the following eigenvalues are obtains for the SK and SG approaches at the Γ -point: 3u There is a difference between the SK and SG TB models for the values of ε K  When the (effective) 1NN hopping is needed the value of u 2,e 4 breaks the symmetry between the two d ±2 -orbitals.However, in the SK model this parameter is set to zero so that these energy levels are always degenerate at the K-point, the SK model can never describe the effective d-orbitals.
The inability to construct a SK model with only d-orbitals was discussed by Cappelluti et al [20].The model from Liu et al [21] was published almost simultaneously.Furthermore, Cappelluti et al mentioned that the fit was not good, it only represented the required orbital characteristics for the TB model.
For the 5th nearest neighbour hopping along the armchair direction, the equivalent u 5,e 4parameter is zero because of the σv -symmetry.The 6th NN hopping goes along the zigzag direction, allowing for another term that breaks the degeneracy of the d ±2 -orbitals at the K point, as mentioned in appendix E. This inclusion of the additional parameters with the 6th NN hopping can contribute to the better description of the local minima in the conduction band at the Q-point.Looking at variations for the u 6M ,e 4 -parameter indeed reveales that it has a strong influence for the values at the Q-point, as discussed in appendix G.
Concluding the previous discussion, we can give a ranking of the different TB models.It is clear that the SK models should be avoided, except when the parameters need to be written in terms of the types of bonds between the orbitals, to describe certain properties like strain [13,38].The good agreement with DFT and the overall simplicity of Liu 6 model means that it is a good start for numerical calculations.Our model also adds the p-orbitals, which makes it better suited for investigating small variations in the structure and large scale calculations.The novel model only includes short range hoppings so that the Hamiltonian stays relatively sparse.The model from Fang has all the orbitals, including the full top-and bottom-p-orbitals, which makes it useful to investigate heterostructures where there are different hoppings between the top-and bottom-p-orbitals.The model from Wu should be used when also the odd d-orbitals are needed, for example when the spin-flip needs to be taken into account.The Liu 2 model should be used when only the effective first nearest neighbour hopping δ max = 1 is needed.These conclusions are summarized in table 3.

Conclusion
We gave an overview of the most important TB models from literature, and proposed a novel TB model.The properties of the TB models and the construction method used determines the precision of the description of the electronic structure.While we used MoS 2 as an example to study these models, the formalism can be expanded to similar TMDs because of the similar electronic structure.Overall, the TB model from Liu 6 gives impressive results around the band gap.The TB model from Fang gives even better overall results, but it has more bands and thus leads to a denser Hamiltonian.
Our newly proposed TB model is more efficient, suitable for large scale calculations where the presence of chalcogen atoms is desired.The models from Liu or Wu do not have the chalcogen atoms in the model, and thus it is hard to model a nanoribbon with specific edges.Our model has the explicit dependence of the chalcogen atoms, while maintaining a minimal amount of hoppings.
Finding a parameterization that matches the symmetry group of the TB model is important.We saw that excluding a specific parameter can lead to a degeneracy that is not expected from the DFT results.This degeneracy is broken if all the parameters allowed by the symmetry group are taken into account.Furthermore, this parameter also helps in obtaining a better fit around the Q point.Our newly proposed model incorporates this additional parameter.
All the models, their new parameters and the parameters from literature are included in a new Python package, TMDybinding [10].This package uses the same construction method for all the TB models, so that parameters can be compared easily between the different models.
It is clear from this analysis, that there is not one go-to model for MoS 2 like there is for graphene.However, one should avoid a TB model that uses the SK approach for any numeric calculations, besides when there is an explicit dependence on the energy for specific type of bonds.with m = (2/6), f = (2/5/6) and n = (1/3/4).For the 4-th NN hopping, the matrices must also be rotated to obtain with γ = arctan 3 5 .The parameters for the hopping matrices for the SG TB models are given in tables 4 and 5.

C SOC matrices
The matrices for the SOC are given by t ↑↑,ee with the following relations between the different parts

D Slater-Koster parameterization
For the Slater-Koster models, the values of the parameters are, with c = tan θ ϵ with r 1 = −1, r 2 = 3, r 3 = 2, r 4 = − 7, r 5 = 3 and r 6 = 2 3.The parameters for the hopping matrices for the SK TB models are given in table 6.

E Hamiltonian at high-symmetry points
At the high-symmetry points Γ and K, the Hamiltonian can be solved exactly by preforming suitable basis transformations [20].The terms on the diagonal u p and u d , get contributions from the onsite energies and the second, fifth and sixth nearest neighbour hoppings (given by f ), the pd terms only get contributions from the first, third and fourth nearest neighbour hoppings (given by n).

E.1 Γ -point
In the 3-band TB model with only even d-orbitals, the second fifth and sixth nearest neighbour hoppings are present.This results in u pd being zero.In this model, there are no p-orbitals, the u p disappears and the Hamiltonian becomes diagonal.The eigenvalues are easily obtained as ε Γ -orbitals are zero, the hoppings along the zigzag direction.The fifth nearest neighbour hopping, along the armchair direction, does not contribute to this splitting because it is protected by the σ v -symmetry.

F All hoppings
To further analyze all the results, we did a fit up to sixth nearest neighbours.This included all the hoppings for all the orbitals.The fit corresponds quite well with the DFT results, as seen in figure 7. The discrepancies of the energy levels align with the energy levels where the s-orbital has a higher orbital contribution, as seen in the DFT results for the DOS in figure 4. The parameters for this SG TB model are given in table 4 and 5.

G Sensitivity of the Q-point
The Q-point is sensitive to longer range hoppings.From figure 3 we can see that the Q-point has a better agreement with the DFT results for TB models that include higher order hoppings.To quantify this better, we investigated the sensitivity of the Liu 6 SG TB model to its parameters.As the case with the band gap opening up, the u 4  6 -parameter seems to play a role.However, this parameter also shifts the location of the location of the energy of the valence band.This

Figure 1 :
Figure 1: (a) Perspective view of the T M D−H structure with the different X b −M −X t indicated in yellow-blue-yellow respectively, (b) side view of the structure, (c) the Brillouin Zone with the high-symmetry points with emphasis on the inequivalent K and K ′ points.
a = 3a M X .For MoS 2 , we obtained the values a = 3.19Å, d X −X = 3.12Å and θ = 0.703 using density function theory (DFT) calculations.The reciprocal-space vectors are b 1 = 2π a x + 1 3 ŷ and b 2 = 4π a 3 ŷ.The special points in the first Brillouin zone (BZ) are at Γ = 0 (the center), K = −K ′ = 1 3 (2b 1 − b 2 ) (the corners), and M = 1 2 b 1 (the center of the edges).The line connecting these points is a line of high-symmetry.Around the Q point, between the Γ and K points, the conduction band has a local minimum.

Figure 3 :
Figure 3: The band structure for the different TB models given in table 1.The DFT results on which the TB models were fitted are given as a reference.The additional odd bands in the Wu model are indicated with a dotted line.

Figure 4 :
Figure 4: The (projected) DOS for DFT results and the different TB models given in table 1, except for the model from Rostami.The DFT results also include the sorbitals.

Figure 5 :
Figure 5: The band structure with SOC, with the projected s z values for the DFT results and the TB models given in table 1, except for the model from Rostami.Between the M and Γ points, the bands are degenerate due to the time reversal symmetry.

Figure 6 :
Figure 6: The Berry phase for the different models as given in table 1 and the DFT results.The models from Liu 6 and Wu have the same valence band and give the same results.The first Brillouin zone is indicated in red.

d.
±2 The value of ε K d 0 is again given by the d z 2 -orbital, the values for ε K d ±2 are given by a linear combination of the d x 2 − y 2 , d x y -orbitals, the chiral basis.At the K-point the bands are flipped, and the valence bands is given by the linear combination d +2 instead of the d z 2 -orbital.In the SG TB model,

Table 5 :
The fitted parameters for the SG TB models, as used in this manuscript.

2 ,
ϵ d 0 originates from the d z 2 -orbital, and the ϵ d 2 is degenerate between the d x yand d x 2 − y 2 -orbitals.E.2 K -pointAt the K-point, the Hamiltonian becomes block-diagonal in the following chiral basis:φ = φ e p e−1 d −again making 2 × 2 blocks and one 1 × 1 block, like explained for the case at the Γ -point.The eigenvalues can also be obtained in the same manner.The same hoppings still contribute to the diagonal or off-diagonal terms.For the 3-band TB model the Hamiltonian again becomes diagonal, the eigenvalues are given by energy levels.The contributions to the last two energy levels are given by the d 2 -states, thus resulting in a superposition of d x 2 − y 2 and d x y .In theory, one can construct a TB model consisting of only two bands: a d z 2 and d +2 -band.However, this will only be correct for the energy levels close to the band gap.It would model an effective d-orbital for the conduction band, thus resulting in an effective TB model of effective d-orbitals.The d-orbitals are degenerate if the u m,M e 4

Figure 7 :
Figure 7: A SG TB model for all the orbitals and hoppings up to sixth nearest neighbours.

Figure 8 :e 4 .
Figure 8: Sensitivity of the Liu 6 SG TB model to changes in the parameters u 4,M e 2 and u 6,M e 4. If both values are scaled so that they compensate each other, only the energy levels around the Q-point change.

Table 3 :
Ranking of TB models, which to consider for specific calculations.

Table 4 :
The fitted parameters for the SG TB models, as used in this manuscript.

Table 6 :
The fitted parameters (in eV ) for the SK TB models, as used in this manuscript.The positions with a dash (−) have parameters that have the same value as the one with the different parity ρ.
At the Γ -point, the Hamiltonian becomes block-diagonal in the following Cartesian basis Each block has contributions for a metal d and a chalcogen p orbital.This allows us to denote these blocks as