Flat band properties of twisted transition metal dichalcogenide homo- and heterobilayers of MoS$_2$, MoSe$_2$, WS$_2$ and WSe$_2$

Twisted bilayers of two-dimensional materials, such as twisted bilayer graphene, often feature flat electronic bands that enable the observation of electron correlation effects. In this work, we study the electronic structure of twisted transition metal dichalcogenide (TMD) homo- and heterobilayers that are obtained by combining MoS$_2$, WS$_2$, MoSe$_2$ and WSe$_2$ monolayers, and show how flat band properties depend on the chemical composition of the bilayer as well as its twist angle. We determine the relaxed atomic structure of the twisted bilayers using classical force fields and calculate the electronic band structure using a tight-binding model parametrized from first-principles density-functional theory. We find that the highest valence bands in these systems can derive either from $\Gamma$-point or $K$/$K'$-point states of the constituent monolayers. For homobilayers, the two highest valence bands are composed of monolayer $\Gamma$-point states, exhibit a graphene-like dispersion and become flat as the twist angle is reduced. The situation is more complicated for heterobilayers where the ordering of $\Gamma$-derived and $K$/$K'$-derived states depends both on the material composition and also the twist angle. In all systems, qualitatively different band structures are obtained when atomic relaxations are neglected.


I. INTRODUCTION
Introducing a twist between two van der Waals stacked two-dimensional materials creates a moiré pattern which results in novel emergent properties. For example, a graphene bilayer with a twist of ∼1.1 degree exhibits flat bands, strong electron correlations and superconductivity which are absent in the constituent monolayers [1][2][3][4][5] . These findings have generated significant interest and established the new field of twistronics 6 .
Besides graphene, there exist many two-dimensional materials that can be used as building blocks of moiré materials 7 . In particular, the transition metal dichalcogenides (TMDs) with chemical formula MX 2 [with M being a transition metal atom such as tungsten (W), molybdenum (Mo), niobium (Nb) or tantalum (Ta) and X denoting a chalcogen atom such as sulphur (S), selenium (Se) or tellurium (Te)] are a promising class of candidate materials. In contrast to graphene, many monolayer TMDs are semiconductors with band gaps in the range of 1-2 eV which makes these materials promising for applications in nano-and optoelectronics [8][9][10][11] . Moreover, monolayer TMDs exhibit strong spin-orbit coupling and spin-valley locking as a consequence of their crystal structure and the presence of heavy transition metal atoms 10,11 .
Recently, several experimental groups have started to explore the properties of twisted TMD bilayers. For example, Wang and coworkers 12 fabricated bilayers of WSe 2 with different twist angles and observed a correlated insulator state when the lowest valence band was half filled with holes. In the same system, Huang et al. 13 measured a giant nonlinear Hall effect at small twist angles, and control of optical properties through twisting has been reported in MoS 2 bilayers 14 .
In addition to homobilayers consisting of two identical TMD monolayers, it is also possible to create heterobilayers consisting of two different TMD monolayers. For heterobilayers, a moiré pattern emerges even without a twist between the layers, as a consequence of the different lattice constants of the constituent monolayers. Tran and coworkers 15,16 studied the optical properties of twisted WSe 2 /MoSe 2 bilayers and observed signatures of interlayer excitons that are trapped by the moiré potential. A similar experiment but with an untwisted WSe 2 /MoSe 2 bilayer was performed by Gerardot and coworkers 17 , who also observed spin-layer locking of interlayer excitons in a 2H-MoSe 2 /MoSe 2 /WSe 2 trilayer 18 . Tang et al. 19 detected interactions between excitons and magnetically ordered holes in angle-aligned WSe 2 /WS 2 structures indicating that this system can be used to simulate the phase diagram of the triangular Hubbard model. The existence of stripe phases over a large doping range has recently been reported in untwisted WSe 2 /WS 2 bilayers by Mak et al. 20 . The same system also shows an abundance of correlated insulating states across a range of electron and hole doping levels 21 .
To understand these experimental findings, detailed knowledge of the electronic structure of twisted TMD bilayers is required. Several groups have carried out density-functional theory (DFT) calculations of twisted homobilayers. For example, Naik and Jain 22 have calculated the band structure of several homobilayers (neglecting the effect of spin-orbit coupling) at a twist angle of 3.5 degree and found flat valence bands. However, accessing smaller twist angles is challenging because of the unfavorable scaling of standard first-principles techniques with system size. To access smaller twist angles, Zhan et al. 23 used the ab initio tight-binding model developed by Fang and coworkers 24 for untwisted homobilayers and calculated the band structure of MoS 2 homobilayers for twist angles as small as 1.6 degree, including the effect of arXiv:2102.03259v2 [cond-mat.mtrl-sci] 17 May 2021 spin-orbit coupling. In a similar work, de Laissardìere et al. 25 used a Slater-Koster based tight-binding approach to study the evolution of flat bands in twisted homobilayer MoS 2 . As an alternative to atomistic methods, Wu and coworkers 26 employed a continuum effective mass approach to study the electronic structure of twisted heterobilayers. Similar work was carried out by Zhang, Yuan and Fu 27,28 and Vogl et al. 29 . A different approach based on generalised Wigner crystals, has been proposed by Phillips and coworkers 30 to explain the emergence of insulating states at fractional filling.
In this work, we systematically study the atomic and electronic structure of all 3R stacked (θ ∼ 0 • ) twisted homo-and heterobilayers that can be constructed by combining MoS 2 , MoSe 2 , WS 2 and WSe 2 monolayers. Specifically, we use classical force fields to calculate the relaxed atomic structure of these systems, which display significant in-plane and out-of-plane relaxations. For the relaxed structures, we use an atomistic tight-binding model derived from first-principles DFT calculations to calculate the electronic band structure including the effect of spin-orbit coupling. In all homobilayers, we find that for relatively small angles (θ < 4 • ), the two highest valence bands are composed of Γ-valley states of the constituent monolayers and become extremely flat as the twist angle approaches zero, reaching bandwidths of a few meV for twist angles near 1.5 • . In contrast, not all heterobilayers exhibit such Γ-derived valence bands. In some heterobilayers (most notably those containing a WSe 2 layer), the top valence bands derive from monolayer states at K and K . Such K/K -derived valence states are less affected by interlayer coupling and are found to be more dispersive compared to Γ-derived states. The different ordering of Γ-derived and K/K -derived valence states in the various twisted bilayer systems can be understood by comparing the energy scale for interlayer hopping with the energy difference between the valence band K-and Γ-states of the constituent monolayers. Importantly, the neglect of atomic relaxations leads to qualitatively different electronic properties.

A. Atomic structure
As a first step, we generate structures of flat (i.e. unrelaxed) twisted TMD homo-and heterobilayers (tBL-TMDs). We start from 3R stacked bilayers, where metal and chalcogen atoms of the top layer are directly above corresponding metal and chalcogen atoms of the bottom layer, and rotate the top layer by an angle θ around the axis perpendicular to the plane of the bilayer and going through the metal atoms. For homobilayers, a commensurate structure is obtained when the moiré cell vectors t 1 and t 2 can be expressed as 25 t 1 = na 1 + ma 2 , t 2 = −ma 1 + (n + m)a 2 , (1) where a 1 = a 2 ( √ 3, 1, 0) and a 2 = a 2 ( √ 3, −1, 0) are primitive lattice vectors of the monolayer (with a being the lattice constant) and m and n are integers. The twist angle is given by cos θ = n 2 +4nm+m 2 2(n 2 +nm+m 2 ) and the number of atoms in the cell is N at = 6(n 2 + nm + m 2 ).
For heterobilayers, we first consider systems whose constituent monolayers contain the same species of chalcogen atom. In this case, the lattice constants of both monolayers differ by less than 1% and we generate a commensurate moiré cell for the twisted heterobilayers by increasing the lattice constant of the monolayer with the smaller lattice constant to the value of the larger lattice constant and then use the same approach described above for homobilayers.
In contrast, for heterobilayers whose constituent monolayers contain different species of chalcogen atom, the lattice constants of the monolayers differ by several percent. To generate moiré cells for these systems, we follow the approach of Zeller and Günther 31 . In their work the moiré vectors t 1 and t 2 are defined as (2) where a 1 = a(1, 0, 0) and a 2 = a 2 (−1, √ 3, 0) are the primitive lattice vectors of the monolayer with smaller equilibrium lattice constant a (a denoting the lattice constant of the other layer). We use DFT equilibrium lattice constants from Ref. 32 . The integers n and m are determined from the numerical solution of a diophantine equation (see Appendix of Ref. 31 for details). Here, we only consider so-called first-order moiré structures 31 . Importantly, to generate a commensurate moiré cell near a desired target twist angle a certain level of strain must be applied. In this work, we only study systems with an overall strain of less than 3%. The strain (which can be either tensile or compressive) is always applied to the layer with the larger equilibrium lattice constant.
For both homo-and hetero-bilayers, using the flat twisted bilayers as starting points, we determine the relaxed equilibrium atomic structure via classical force fields as implemented in the LAMMPS software package 33,34 .
In particular, we employ the force fields developed by Naik and coworkers, based on the Kolmogorov-Crespi potential for the interlayer interaction 35 . For intralayer interactions, a Stillinger-Weber type force field is used 36 . The relaxed structures of most twisted TMD bilayers are of a breathing-mode type, i.e., the two layers have out-of-plane displacements in opposite directions (see Sec. III A). The only exception are heterobilayers with different chalcogen atoms. These systems also exhibit breathing-mode relaxations for twist angles ≥ 4.5°, but for smaller twist angles qualitatively different relaxed structures are found in which the two layers have out-of-plane displacements in the same direction and the amplitude of these displacements is larger compared to those in breathing mode structures. In an experimental setting, we expect such structures to be less likely to occur because the twisted TMD bilayers are placed on a substrate. Therefore, we do not present results for these systems and focus our attention on breathing mode structures.

B. Electronic structure
To calculate the electronic properties of twisted TMD homo-and heterobilayers, we use an atomistic tightbinding approach based on the work of Fang and coworkers 24 , who studied untwisted homobilayers. The atomic orbital basis consists of 5 d-like orbitals for the metal atoms and 3 p-like orbitals for each chalcogen atom (which doubles to 10 d-like orbitals and 6 p-like orbitals, respectively, if spin-orbit coupling is included). In a first step, we construct a symmetry-adapted tight-binding model for the monolayers including on-site, first, second and selected third nearest neighbor hoppings. The required hopping parameters are determined from a Wannier transformation 37,38 of the DFT Hamiltonian. To model bilayers, Fang and coworkers describe interlayer hoppings between the p-orbitals of the chalcogen atoms at the interface between the two layers, which we refer to as inner chalcogens, using the Slater-Koster approach 24 . The Slater-Koster parameters are fitted to a set of DFT calculations of untwisted bilayers in which the top layer is translated horizontally while the bottom layer is kept fixed. Finally, spin-orbit coupling is introduced via an on-site atomic term λ SO M/X L · S (with L and S denoting orbital and spin angular momentum operators, respectively, and λ SO M/X is the spin-orbit coupling strength of M or X atoms, whose value for each atom is given in Ref. 24 .) To model twisted homo-and heterobilayers, we have extended the tight-binding model of Fang et al. 24 in several ways. In particular, we have included interlayer hoppings from inner chalcogen p z -like orbitals on one layer to metal d z 2 -like orbitals on the other layer using a Slater-Koster approach. Moreover, to better capture the effect of out-of-plane displacements of the atoms, we improve the description of interlayer hoppings (both pp and p z -d z 2 ) by using a different set of Slater-Koster parameters for different values of the interlayer separation. All Slater-Koster parameters for the interlayer interactions as well as all intralayer hoppings were obtained from a Wannier transformation 37,38 of the DFT Hamiltonian. For heterobilayers, additional care must be taken to ensure that the on-site energies are referenced to the vacuum level.
To determine the interlayer hoppings in a twisted bilayer, the orbital basis of the rotated monolayer must be transformed. As described above, only p x -, p y -, p zand d z 2 -like orbitals are involved in interlayer hoppings. Since d z 2 -like and p z -like orbitals are unaffected by rotations around the z-axis, we only need to transform the p x -like and p y -like orbitals and the rotated orbitals are given by (p x , p y ) T = R(θ)(p x , p y ) T with R(θ) denoting a two-dimensional rotation matrix. Of course, interlayer hoppings involving p x -like and p y -like orbitals transform in a similar fashion when a twist is introduced.
Additional details about the interlayer tight-binding model, the determination of the hopping parameters and a full list of the parameters for all systems can be found in the Appendix A and in the Supplementary Information (Sec. S5). We have compared the band structures (without spin-orbit coupling) from this tight-binding model to results from explicit DFT calculations for different twisted bilayers, see Sec. S2 in the Supplementary Information, and find good agreement between the two methods, in particular for the valence bands.
Besides modulating the interlayer hopping, the introduction of a twist also gives rise to significant in-plane atomic relaxations which in turn induce changes in the intralayer hoppings. Such changes, however, are not captured by our model as we assume that intralayer hoppings of the twisted bilayer are the same as those in a monolayer. Recently, it was shown that such twist-induced changes to the intralayer hoppings are responsible for the flattening of K/K'-derived valence band states in WSe 2 /WS 2 superlattices 39 . To capture this effect, a fully position-dependent intralayer tight-binding Hamiltonian for the TMDs should be developed in the future.

A. Atomic structure
Introducing a twist between two 3R aligned TMD layers results in the creation of a moiré pattern consisting of regions with different stacking arrangements. Highsymmetry stackings include AA regions, where the metal (chalcogen) atoms of the one layer are directly above the metal (chalcogen) atoms of the other layer, as well as two types of Bernal-like regions, where in one case (denoted B M/X ) the chalcogen atom (X) in one layer lies directly above the metal atom (M) in the other layer, or vice versa (denoted B X/M ). The high-symmetry stackings are shown in Fig. 1(a).

Homobilayers.
All twisted homobilayers exhibit similar out-of-plane and in-plane displacement patterns upon relaxation. For example, Figs. 1(a)-(d) show results for twisted MoSe 2 /MoSe 2 at a twist angle of θ = 2.6°. The interlayer separation (ILS), defined as the distance between the two surfaces on which the inner chalcogen atoms lie, is large in the AA regions (which form a triangular lattice), but smaller in the triangle-shaped B M/X and B X/M regions (which form a honeycomb lattice), see Fig. 1(a). Fig. 1(b) shows the out-of-plane displacement along the diagonal of the moiré unit cell for three different twist angles. As the twist angle decreases, the size of the AA regions shrinks, whereas B M/X and B X/M regions expand. This allows the system to reduce its energy as AA regions are energetically unfavorable because of their large steric repulsion. It can further be observed that the maximum ILS increases, while the minimum ILS decreases as the twist angle is reduced. Again, this reduces the energy cost associated with steric repulsion.   Figure 2 shows the in-plane and out-of-plane relaxations of twisted MoSe 2 /MoS 2 at a twist angle of θ = 4.5°. For the set of angles studied in this work, we find that heterobilayers exhibit similar relaxation patterns as homobilayers: large ILSs are found in the AA regions, which form a triangular lattice. The relative size of the AA regions shrinks as the twist angle is decreased while B M/X and B X/M regions grow. In contrast to the homobilayers, the in-plane and out-of-plane displacements of the two layers are not symmetric, as can be seen in Figs. 2(b),(c) and (d). The difference of the out-of-plane displacements in the AA and B M/X regions is about four times larger in the MoS 2 layer than in the MoSe 2 layer (Fig. 2(b)). As we show in the next section, this asymmetry is less pronounced in heterobilayers    Fig. 3(a)) decreases monotonically as the twist angle is reduced and converges to the ILS of the untwisted bilayers. In contrast, the ILS in the AA regions (top panel of Fig. 3(a)) increases with decreasing twist angle, but does not converge to the value of the untwisted bilayer in the case of homobilayers and heterobilayers with same chalcogen atoms. This discontinuity of the maximum ILS at θ = 0 • is a consequence of the structural relaxations which result in a growth of the B M/X and B X/M regions and a shrinkage of the AA regions at small twist angles. At the center of the large B M/X /B X/M regions the twisted bilayer has a similar structure as the untwisted B M/X /B X/M bilayer while the small size of the AA restricts the atoms from reaching the same structure as the untwisted AA bilayer.
Comparing the ILS of different bilayers, we observe that bilayers where both constituent monolayers contain S atoms (MoS 2 /MoS 2 , WS 2 /WS 2 and WS 2 /MoS 2 ) exhibit the smallest interlayer distances (both in AA and B M/X /B X/M regions), whereas bilayers containing Se atoms (MoSe 2 /MoSe 2 , WSe 2 /WSe 2 and WSe 2 /MoSe 2 ) in both layers exhibit the largest ILSs. Bilayers with S atoms in one layer and Se atoms in The out-of-plane displacements for all homobilayers and heterobilayers with same chalcogen species at θ = 5.1°and for all heterobilayers with different chalcogens at θ = 5.4°are shown in Fig. 3(b). As shown in the left panel of Fig. 3(b), out-of-plane displacements in homobilayers are layer-symmetric and the shape of the dis-placement patterns is similar for all systems. In contrast, out-of-plane displacements in heterobilayers ( Fig. 3(b) right panel) are layer-asymmetric. In these systems, the amplitude of the displacement pattern of the bottom layers (which are unstrained) is similar to that found in the homobilayers, while the amplitudes of the strained top layer are somewhat smaller. In this section we study the evolution of the band structure of the twisted homobilayers MoS 2 /MoS 2 , MoSe 2 /MoSe 2 , WS 2 /WS 2 and WSe 2 /WSe 2 as function of twist angle. All calculations were carried out for the relaxed structures and include spin-orbit coupling. For all twist angles, the homobilayers exhibit a semiconducting band structure with a band gap separating the valence and conduction bands. Moreover, the two highest valence bands (each of which is spin degenerate) are separated from all other "remote" valence bands by energy gaps when θ < 4 • . We refer to these two highest valence bands as VB1 and VB2, respectively.
As an example, we focus on the valence band structure of twisted MoS 2 /MoS 2 . Fig. 4  structure at θ = 5.1°. It can be seen that the valence band maximum occurs at the Γ-point of the moiré Brillouin zone and that VB1 and VB2 touch at the K-point forming a Dirac cone. This is also reflected in the Vshaped density of states ( Fig. 4 (b)). The graphene-like dispersion of VB1 and VB2 can be understood by analyzing the wavefunctions of these states. Fig. 4(d) shows that the wavefunctions at Γ are localized in the B M/X and B X/M regions which form a honeycomb lattice. Importantly, the total bandwidth of the two highest valence bands is less than 30 meV demonstrating the formation of flat bands upon twisting.
To understand the chemical origin of the flat bands, we analyze their projections onto atomic orbitals. Fig. 4(c) shows that these states are localized symmetrically on the inner layers of chalcogen atoms and also on the two metal layers. These states are mostly composed of inner chalcogen p z -like orbitals and metal d z 2 -like orbitals, as show in Fig. S5 of the Supplementary Information. This suggests that VB1 and VB2 originate from Γ-states of the constituent monolayers: in all TMD monolayers, the top valence band states at Γ have large projections onto chalcogen p z -like orbitals and metal d z 2 -like orbitals, whereas the top valence band states at K and K have large projections onto metal d xy -like and d x 2 −y 2 -like orbitals 24 .
When atomic relaxations are not taken into account, a qualitatively different valence band structure is obtained. Figs. 5(a) and 5(b) compare the band structures of unrelaxed and relaxed MoS 2 /MoS 2 at θ = 5.1 • . In contrast to the relaxed system, the unrelaxed system does not exhibit a Dirac-like dispersion and exhibits an energy gap between VB1 and VB2. This is a consequence of the different spatial structure structure of the corresponding wavefunction, see Figs. 5(c) and 5(d): in the relaxed system VB1 and VB2 localize in the B M/X and B X/M regions, while in the unrelaxed system the top valence band states are localized in the AA regions which form a triangular lattice 26 . We find a similar effect of atomic relaxations in all homobilayers, see Supplementary Information (S1). Figure 6 compares the band structures of all homobilayers at three twist angles (θ = 5.1°, 2.6°and 1.6°). At θ = 5.1°, the top valence band in twisted MoS 2 /MoS 2 and MoSe 2 /MoSe 2 are Γ-derived and exhibit a Dirac-like dispersion with a valence band maximum at Γ, as discussed above. In contrast, for WS 2 /WS 2 and WSe 2 /WSe 2 the Γ-derived valence states are intersected by dispersive bands which are derived from monolayer K/K -states (see discussion below). In WSe 2 /WSe 2 , the highest valence band is K/K -derived and the valence band maximum is found at the K-point.
When the twist angle is reduced to θ = 2.6°, the Γ-derived bands become flatter and the K/K -derived bands are shifted to lower energies such that they no longer intersect the flat Γ-derived top valence bands. Decreasing the twist angle further to θ = 1.6°, we observe that the highest four remote Γ-derived valence bands become isolated in energy from all other remote bands, see Figs. 6(c),(f),(i),(l). The two middle bands of this set also exhibit a Dirac-like dispersion near the K-point, while the highest and lowest bands are very flat. Our results are in good agreement with DFT calculations performed by Xian and coworkers 42 , who also analyzed the origin of these bands and proposed that they can be described by a set of p x -like and p y -like orbitals on a honeycomb lattice.
Figure 7(a) shows the bandwidth w, computed as the energy difference between states at Γ and K, of the top valence band (denoted as VB1 as in Fig. 4(b)) as function of twist angle for relaxed and unrelaxed (flat) homobilayers. As the twist angle decreases, the bandwidths approach zero. For relaxed homobilayers, the magnitude of the bandwidths in the different systems are relatively similar with MoSe 2 /MoSe 2 exhibiting the smallest one (reaching ≈ 0.5 meV at θ = 1.6°). When relaxations are neglected, the bandwidths are smaller. For example, a bandwidth of only 0.2 meV is found in unrelaxed MoSe 2 /MoSe 2 at θ = 2.6°. Figure 7(b) shows the band gap ∆ between valence and conduction states as function of twist angle for both relaxed and unrelaxed homobilayers. For all relaxed homobilayers ∆ decreases linearly with a slope of ≈ 20 meV/degree as the twist angle is reduced. for three twist angles θ = 5.1°, 2.6°and 1.6°. The high-symmetry path Γ-M -K-Γ is shown in Fig. 4(a).
B X/M ) bilayers, shown on the left panel of Fig. 7(b) at θ = 0°. This is expected as the B M/X and B X/M regions are energetically favorable (compared to AA regions) and their relative size grows as the twist angle is reduced (see Sec. II A).
Interestingly, the nature of the band gap of WSe 2 /WSe 2 changes from direct (K → K) to indirect (Γ → K) around θ = 5.1°. This is a consequence of the change in ordering of Γ-derived and K/K -derived valence states, see Fig. 6. All other systems exhibit indirect band gaps. In particular, for MoS 2 /MoS 2 and MoSe 2 /MoSe 2 the valence band maximum is at Γ and the conduction band minimum at K as in the untwisted case, while for WS 2 /WS 2 the conduction band minimum is half-way between the Γ-point and the M -point (referred to as the X-point), which explains the deviation from the untwisted case.
Without relaxations, the band gaps are almost constant and do not depend sensitively on the twist angle, see Fig. 7(b). Also, the nature of the band gap for WSe 2 /WSe 2 is different compared to the relaxed systems for θ < 5.1°, as is for WS 2 /WS 2 at θ = 2.6°. Results are presented for three twist angles: θ = 5.1°, 2.6°and 1.6°.

Heterobilayers.
We first consider heterobilayers with the same chalcogen species in each layer, i.e., WS 2 /MoS 2 and WSe 2 /MoSe 2 . As discussed in Sec. II, it is possible to generate commensurate moiré structures with very little strain for these systems. The band structures of these systems at three different twist angles are shown in Fig. 8. Again, we find that these systems feature both flat Γderived valence bands as well as dispersive K/K -derived states.
In WSe 2 /MoSe 2 at θ = 5.1 • , the valence band maximum corresponds to a K/K'-derived state. At smaller twist angles, the K/K -derived states are shifted to lower energies and the top valence bands are Γ-derived and very flat. Similar to the homobilayers, the Γ-derived top two valence bands are separated from all other remote valence bands at small twist angles. However, these bands no longer have a Dirac-like dispersion, but are separated by an energy gap at K. This energy gap caused by the "chemical asymmetry" of the two constituent layers.
In contrast, the highest valence bands in WS 2 /MoS 2 are derived from monolayer K/K -states at all twist angles.
Next, we first consider the heterobilayers MoSe 2 /MoS 2 and MoSe 2 /WS 2 , i.e. hetero-bilayers containing different species of chalcogens but no WSe 2 . Information. In contrast to the homobilayers and the heterobilayers containing a single chalcogen species, the top valence bands in these systems are not spin-degenerate. Fig. 10(a) shows the spin-resolved dispersion of top valence bands in MoSe 2 /MoS 2 at θ = 4.5°which exhibits spin splittings with magnitudes up to 13 meV. Also, an energy gap between VB1 and VB2 of 8 meV is found at the K-point. Interestingly, these bands are partially spinpolarized: it can be observed that the top valence bands are only spin-polarized in the vicinity of the K-point even though spin splitting occurs along the whole band structure path. This scenario was recently discussed by Liu and coworkers 43 who demonstrated that spin-orbit cou- pling can lead to spin splittings without spin polarization in non-magnetic materials without inversion symmetry. Finally, Figs. 9(e)-(h) show the band structures of twisted WSe 2 /MoS 2 and WSe 2 /WS 2 bilayers, i.e., the heterobilayers with different chalcogen atoms that contain WSe 2 . The top valence bands in these systems are dispersive and the flat bands are observed at lower energies, see discussion in Sec. III C. Figure 12 compares the bandwidths and band gaps of the different heterobilayers as function of twist angle. Note that we only show results for twisted bilayers whose top valence bands are flat, i.e., no results are shown for WS 2 /MoS 2 , WSe 2 /MoS 2 and WSe 2 /WS 2 . For WSe 2 /MoSe 2 flat bands are only found for θ < 5.1°. Similar to the case of the homobilayers, the bandwidths of the heterobilayers decrease monotonically as the twist angle approaches zero and the value of the bandwidth at a fixed twist angle is roughly the same for the different homo-and heterobilayers.
The band gap of the twisted heterobilayers does not depend sensitively on twist angle, see Fig. 12(b). For most systems, a mild reduction of the gap can be observed as the twist angle decreases. The smallest bands gaps (≈ 0.8 eV) are found for WSe 2 /MoS 2 , while MoSe 2 /WS 2 and WS 2 /MoS 2 exhibit the largest gaps (≈ 1.3 − 1.4 eV). Interestingly, the nature of the band gap depends sensitively on the twist angle and many systems exhibit a change from a direct to an indirect gap as the twist an-

C. Physical origin of the dispersive valence bands
We now focus on the set of heterobilayers that exhibit dispersive bands at the valence band edge and discuss the origin of these bands. As described above, such bands are observed near the top of the valence band in WS 2 /WS 2 and WSe 2 /WSe 2 homobilayers at large twist angles and form the top valence states at all twist angles in WS 2 /MoS 2 , WSe 2 /MoS 2 and WSe 2 /WS 2 , see Figs. 9(e)-(h). Compared to the flat bands, the width of these bands decreases much less as the twist angle is reduced. Projecting the corresponding states onto atomic orbitals reveals large contributions from W d xy -like and d x 2 −y 2 -like orbitals, see Fig. S8 of Supplementary Information, suggesting that they originate from states at the K-point of WS 2 monolayer in the case of WS 2 /MoS 2 and similarly from K-point states of the WSe 2 monolayer in the WSe 2 /MoS 2 and WSe 2 /WS 2 systems 24 . As these states have very small projections onto the inner chalcogen atoms (see Fig. S9 of the Supplementary Information), they are only weakly affected by interlayer coupling which explains why introduction of a twist does not result in a significant reduction of their band width in our model. However, it has been established that such K/K -derived bands can become flat as a consequence of modulations in the intralayer hopping induced by inplane relaxations 39 . As explained in Sec. II, our current tight-binding approach does not capture such variations of the intralayer hopping and therefore does not capture this additional band flattening mechanism.
The top valence band at the K-point of the WS 2 and WSe 2 monolayers is spin-polarized and the dispersive valence bands of the twisted bilayers inherit this property as shown for WSe 2 /MoS 2 in Fig. 10(b). The valence band maximum in bilayers with dispersive bands is located at the K-point. In the WS 2 /MoS 2 bilayer the band gap is direct, whereas it is indirect in the WSe 2 /MoS 2 and WSe 2 /WS 2 bilayers as the conduction band minimum is at Γ. We have found that the ordering of the flat and dispersive valence bands depends sensitively on the atomic structure and the twist angle. For example, Fig. 11(a) shows that neglecting atomic relaxations results in a different ordering of dispersive and flat valence bands. Moreover, we have found that it is possible to switch the order of flat and dispersive bands when the interlayer separation of the relaxed structures is reduced suggesting that the electronic properties of these materials can be easily tuned by applying pressure, as show in Fig. S10 of the Supplementary Information.
To understand the ordering of flat and dispersive valence bands in the different heterobilayers, we propose a simple model in which the bilayer states originating from K/K -or Γ-valleys of the constituent monolayers are obtained from a two-level system. Specifically, the Hamiltonians for the coupled valleys are given by where (1) Γ(K) denotes the energy at Γ (K) of the monolayer with the higher-lying valence band (and (2) Γ(K) denotes the corresponding energies for the monolayer with the lowerlying valence band). Also, ∆ Γ(K) describes the interlayer coupling. As the wavefunctions in the K/K -valley are predominantly localized on the metal atoms, we assume that they are not affected by interlayer coupling and use ∆ K = 0. In contrast, the wavefunctions of Γ-valley states have projections onto chalcogen atoms and these states are pushed to higher energies by interlayer coupling 32 . To calculate ∆ Γ , we carry out tight-binding calculations for untwisted bilayers with B M/X stacking. Then, ∆ Γ is cho-sen such that the largest eigenvalue max Γ of H Γ is equal to the energy of the highest tight-binding valence band state at Γ of the bilayer. For heterobilayers with different chalcogen atoms, we compute the value of max Γ and ∆ Γ by averaging the results from two calculations: one in which the lattice constants of both layers are set equal to max(a 1 , a 2 ), and one in which the lattice constants of both layers are set to min(a 1 , a 2 ), where a 1 and a 2 are the equilibrium lattice constants of the two monolayers, respectively. Table I shows the results from this analysis. For the MoS 2 /MoS 2 and MoSe 2 /MoSe 2 homobilayers, we find that max Γ is significantly larger than (1) K indicating that the flat bands have higher energies than the dispersive bands. This is in agreement with our explicit band structure calculations, see Fig 6. For the WS 2 /WS 2 and WSe 2 /WSe 2 homobilayers, max Γ is only slightly larger than (1) K and we expect that flat and dispersive bands have similar energies. Again, this is consistent with our explicit calculations which show that the ordering of flat and K/K -derived bands can depend on the twist angle for these systems, see Fig. 6. Considering the heterobilayers, we find that Γ-states are predicted to lie above the K/K -derived states in WSe 2 /MoSe 2 , MoSe 2 /MoS 2 and MoSe 2 /WS 2 , while the opposite ordering is predicted for WS 2 /MoS 2 , WSe 2 /MoS 2 and WSe 2 /WS 2 . These predictions are again in agreement with our explicit band structure calculations.

IV. CONCLUSIONS
We have studied the electronic band structures of all twisted transition metal dichalcogenide (TMD) bilayers   and WSe 2 monolayers. Specifically, we have carried out tight-binding calculations taking into account the effect of atomic relaxations and also spin-orbit coupling. In all twisted homobilayers, the top valence bands are derived from monolayer states at Γ and become flat when the twist angle decreases. For twisted heterobilayers, we find two scenarios: either the highest valence band derives from Γ-states of the monolayer and becomes flat upon twisting or it derives from K-and K -states of the monolayer and remains dispersive even at small twist angles. Interestingly, the ordering of flat and dispersive bands depends sensitively on the atomic structure of the bilayer and can be changed by applying pressure. Our findings reveal that the chemical complexity of the twisted TMD bilayers can be harnessed to design flat band properties and pave the way to understanding electron-electron interaction effects in these materials.

ACKNOWLEDGMENTS
We wish to thank Arta Safari for useful discussion on the construction of heterobilayers moiré cells. All authors acknowledge funding from the EPSRC grant E/P77380. Via our membership of the UK's HEC Materials Chemistry Consortium, which is funded by EPSRC (EP/L000202, EP/R029431, EP/T022213), this work used the the ARCHER UK National Supercomputing Service. This work also used the Cirrus UK National Tier-2 HPC Service at EPCC (http://www.cirrus.ac.uk) funded by the University of Edinburgh and EPSRC (EP/P020267/1). Finally, we acknowledge the Imperial College London Research Computing Service (DOI:10.14469/hpc/2232) for the computational resources used in carrying out this work. .

Appendix A: Improvements on the Tight-binding Model
In this Appendix, we describe the modifications that were required to generalize the tight-binding model for untwisted TMD bilayers developed by Fang and coworkers 24 to twisted TMD bilayers.
In particular, we improved the description of interlayer hoppings and parametrized these hoppings for all combinations of homo-and heterobilayers composed of MoS 2 /MoS 2 , MoSe 2 /MoSe 2 , WS 2 /WS 2 and WSe 2 /WSe 2 monolayers.
To calculate the band structure of untwisted TMD bi- layers, Fang and coworkers used Slater-Koster expressions for the interlayer hoppings between chalcogen porbitals. They also included a term for the interlayer hopping between chalcogen p z -orbitals and transition metal d z 2 -orbitals, but did not describe this with a Slater-Koster expression. Instead, they only calculated the value of this hopping for the specific geometry of an untwisted 2H bilayer.
To generalize the description of p z to d z 2 hoppings to twisted bilayers, we used the Slater-Koster formula t pz,d z 2 (r) = n n 2 − 1 2 (l 2 + m 2 ) V pdσ (r) (A1) where the directional cosines are defined as l = r x /r, m = r y /r and n = r z /r. To determine the functions V pdσ (r) and V pdπ (r), we calculated t pz,d z 2 and also t pz,dxz and t pz,dyz for a set of untwisted bilayers with different stacking configurations and different interlayer separations using a Wannier transformation of the DFT Hamiltonian. Next, a least square fitting process was used to extract V pdσ and V pdπ at different interatomic distances and the results were fitted to functions of the type with V, α, β and γ denoting fitting parameters, and d = 3.5Å is an average interlayer distance. Fig. 13 (d) demonstrates that this yields an accurate description of the calculated t pz,d z 2 hopping parameters. We have also tested the influence of other hoppings between chalco-gen p-orbitals and transition metal d-orbitals, but found that the most important contribution arises from p z to d z 2 hoppings.
The accuracy of the Slater-Koster approximation, in describing the orientation dependence of the interlayer hopping integrals for p-p and p-d orbitals, is demonstrated in Fig. 13 for hoppings extracted from displaced untwisted bilayers. The main error for p-p hoppings (see Fig. 13a-c) arises from approximating the orthogonal Wannier basis as non-orthogonal atomic-like orbitals to make use of the Slater-Koster rules. This discrepancy, however, is very small, i.e. less than 20 meV on average for interlayer hopping matrix elements which suggests that it is an appropriate model to describe various configurations seen in twisted bilayers. Figure 14 compares the band structures of untwisted and twisted bilayer MoS 2 at a twist angle of 7.3°from tight-binding with and without the p z to d z 2 hopping to a first-principles density-functional theory result. We find that inclusion of p z to d z 2 hoppings improves the agreement with the first principles result significantly. In particular, the valence band states near Γ are pushed to higher energies which reduces the band gap by approximately 200 meV. A similar shift of the highest valence bands is found also in the twisted bilayers.
Aside from the inclusion of interlayer p z to d z 2 hoppings, we also discovered that the description of interlayer hoppings between chalcogen p-orbitals developed by Fang and coworkers required improvements to obtain an accurate description of twisted bilayers. To parametrize the corresponding Slater-Koster expressions, Fang et al. carried out first-principles calculations of untwisted bilayers with different stacking configurations, but using   We have found that these changes in the interlayer separation are not well captured by the simple Slater-Koster expressions used by Fang and coworkers. Fig. 15(a) shows the Slater-Koster parameters V ppσ and V ppπ as function of the interatomic distance for different interlayer separations: for small interatomic distances, the Slater-Koster parameters depend sensitively on the interlayer separation.
To account for the dependence of the Slater-Koster parameters on the interlayer separation, the following procedure is used: for a given pair of atoms, we first calculate the interlayer separation as the difference of their z-coordinates as well as the interatomic distance. Then, the Slater-Koster parameters for the specific interlayer separation are used to obtain the desired hopping matrix element.
The DFT calculations for monolayers and untwisted bilayers were performed with Quantum Espresso 44 within the optB88 generalized gradient approximation for the exchange-correlation potential and a plane-wave cutoff value of 70 Ry (≈950 eV). Monolayer calculations were performed with a 25 × 25 × 1 Monkhorst-Pack k-point grid whereas a 12 × 12 × 1 k-point grid was used for the untwisted bilayer calculations. The DFT Hamitonian was transformed into the basis of 11 (22) Wannier functions consisting of atomic-like p and d orbitals with average spreads around 2.2Å for the monolayer (bilayer) calculations. For twisted bilayers DFT calculations have been carried out with Onetep 45,46 , a linear-scaling DFT code.
We use the Perdew-Burke-Ernzerhof exchange-correlation functional 47 with projector-augmented-wave pseudopotentials 48,49 , generated from ultra-soft pseudopotentials 50 , and a kinetic energy cutoff of 800 eV. A basis consisting of 9 nonorthogonal generalized Wannier functions (NGWFs) for calcogen atoms and 13 NGWFs for metal atoms is employed. The NGWFs' radii are set to 9.0 a 0 .  . S1 shows the band structures of unrelaxed (flat) and relaxed homobilayers for θ = 5.1°. The main effect of relaxation is to shift up in energy the flat bands (compared to the K/K'-derived bands) and to close the gap at the K-point, giving rise to Dirac cones. Moreover, in WS 2 /WS 2 and WSe 2 /WSe 2 flat and K/K'-derived bands intersect. Similarly, unrelaxed (flat) and relaxed band structures for heterobilayers with same chalcogen species for θ = 5.1°are shown in Fig. S2. The ordering of K/K'-derived vs flat bands is reversed by relaxation. In particular, in unrelaxed WS 2 /MoS 2 , K/K'-derived bands and flat bands intersect and the valence band maximum (VBM) is at the K/K'-point, wherease in the relaxed structure the VBM is at the Γ-point. In WSe 2 /MoSe 2 the situation is reversed. In both cases, the gap between the flat bands at the K-point vanishes. Fig. S3 shows the band structures for unrelaxed and relaxed heterobilayers with different chalcogen species for θ = 4.5°. In MoSe 2 /WS 2 and WSe 2 /MoS 2 the effect of relaxation is to shift up in energy the flat bands and make them separated from all remote bands. For these systems, the gap at the K point remains finite and each band is not doubly degenerate. In WSe 2 /MoS 2 and WSe 2 /WS 2 bilayers the K/K'-derived bands are always on top. Relaxation has a minor effect for these systems, as K/K'-derived bands and flat bands are separated by several hundreds meV.  Fig. S4 shows the comparison between DFT and tight-binding (TB) band structures without spin-orbit coupling (SOC) for relaxed homobilayers at two twist angles, θ = 7.3°and θ = 5.1°, respectively. DFT band structure have been computed with Onetep[1, 2], a linear-scaling DFT code. We use the Perdew-Burke-Ernzerhof exchange-correlation functional [3] with projector-augmented-wave pseudopotentials [4,5], generated from ultra-soft pseudopotentials [6], and a kinetic energy cutoff of 800 eV. A basis consisting of 9 non-orthogonal generalized Wannier functions (NGWFs) for chalcogen atoms and 13 NGWFs for metal atoms is employed. We set the NGWFs' radii to 9.0 a 0 . Agreement between DFT and TB band structures improves at smaller angles.

S4 -Wavefunction densities
Fig . S9 shows the square of the wavefunction of VBM at the K-point (|ψ K | 2 ) in WSe 2 /MoS 2 for θ = 4.5°. This state is localized mainly on B M/X regions, which form a triangular lattice, and is more diffuse than flat bands states. Figure S9: Top view of |ψ K | 2 of VBM in WSe 2 /MoS 2 at θ = 4.5°. Colored dots refer to different stacking regions as described in Fig. 1(a) of main text and moiré cell indicated by grey lines.

S7 Effect of interlayer distance on flat bands
Fig. S10 shows the effect of interlayer separation (ILS) on the ordering of flat bands over K/K'-derived bands in WSe 2 /WS 2 for θ = 4.5°. To compute the band structures in Fig. S10, we start from a relaxed WSe 2 /WS 2 bilayer and rigidly translate the layers in the z-direction without further relaxing them. For ILS greater or equal to the equilibrium ILS we find K/K'-derived bands at the top of the valence band edge (top panels in Fig. S10). As the ILS is reduced (bottom panels in Fig. S10), flat bands emerge at the top of the valence band edge, which suggests that pressure could potentially be used to change the ordering between flat bands and K/K'-derived bands, and consequently, being able to modify flat-band properties in these systems. Starting from a relaxed structure, the ILS is increased/decreased by rigidly moving the layers without further relaxing the system. In panel (a) the ILS is increased by 0.5Å. In panel (b) the ILS is increased by 1.0Å. In panel (c) the ILS is decreased by 0.5Åand in panel (d) the ILS is decreased by 1.0Å. In all panels the band structure at equilibrium is shown with violet circles and that with modified ILS with solid orange lines. The ordering of flat vs K/K'-derived bands is reversed when the two layers are closer, compared to the equilibrium case (which exhibit K/K'-derived bands on top), as shown in bottom panels (which exhibit flat bands on top).

S5 -Tables of tight-binding parameters
In this section we present tables of tigh-binding hopping parameters, both intralayer and interlayer, for all homo-and hetero-bilayers in our work.