Flat band properties of twisted transition metal dichalcogenide homo- and heterobilayers of MoS_2, MoSe_2, WS_2 and WSe_2
Valerio Vitale, Kemal Atalar, Arash A. Mostofi, Johannes Lischner
FFlat band properties of twisted transition metal dichalcogenide homo- andheterobilayers of MoS , MoSe , WS and WSe Valerio Vitale , Kemal Atalar , Arash A. Mostofi , and Johannes Lischner Departments of Materials and Physics and the Thomas Young Centre for Theory and Simulation of Materials,Imperial College London, London SW7 2AZ, UK
Twisted bilayers of two-dimensional materials, such as twisted bilayer graphene, often feature flatelectronic bands that enable the observation of electron correlation effects. In this work, we studythe electronic structure of twisted transition metal dichalcogenide (TMD) homo- and heterobilayersthat are obtained by combining MoS , WS , MoSe and WSe monolayers, and show how flatband properties depend on the chemical composition of the bilayer as well as its twist angle. Wedetermine the relaxed atomic structure of the twisted bilayers using classical force fields and calculatethe electronic band structure using a tight-binding model parametrized from first-principles density-functional theory. For homobilayers, we find that the two highest valence bands exhibit a graphene-like dispersion and become flat as the twist angle is reduced. In contrast, not all heterobilayers haveflat valence bands. Specifically, we find that those systems in which the highest valence band derivesfrom K or K’ states of the constituent monolayers do not exhibit flat bands, even at small twistangles. In all systems, qualitatively different band structures are obtained when atomic relaxationsare neglected. I. INTRODUCTION
Introducing a twist between two stacked two-dimensional materials creates a moir´e pattern which re-sults in novel emergent properties. For example, agraphene bilayer with a twist of ∼ . ? –4 .These findings have generated significant interest and es-tablished the new field of twistronics .Besides graphene, there exist many two-dimensionalmaterials that can be used as building blocks of moir´ematerials . In particular, the transition metal dichalco-genides (TMDs) with chemical formula MX (with M be-ing a transition metal atom such as tungsten (W), molyb-denum (Mo), niobium (Nb) or tantalum (Ta) and X de-noting a chalcogen atom such as sulphur (S), selenium(Se) or tellurium (Te)) are a promising class of candi-date materials. In contrast to graphene, many monolayerTMDs are semiconductors with band gaps in the range of1-2 eV which makes these materials promising for appli-cations in nano- and optoelectronics . Moreover, mono-layer TMDs exhibit strong spin-orbit coupling and spin-valley locking as a consequence of their crystal structureand the presence of heavy transition metal atoms.Recently, several experimental groups have started toexplore the properties of twisted TMD bilayers. For ex-ample, Wang and coworkers fabricated bilayers of WSe with different twist angles and observed a correlated in-sulator state when the lowest valence band was half filledwith holes. In the same system, Huang et al. measureda giant nonlinear Hall effect at small twist angles, andcontrol of optical properties through twisting has beenreported in MoS bilayers .In addition to homobilayers consisting of two identicalTMD monolayers, it is also possible to create heterobi-layers consisting of two different TMD monolayers. For heterobilayers, a moir´e pattern already emerges withouta twist as a consequence of the different lattice constantsof the constituent monolayers. Tran and coworkers studied the optical properties of twisted WSe /MoSe bi-layers and observed signatures of interlayer excitons thatare trapped by the moir´e potential. A similar experi-ment but with an untwisted WSe /MoSe bilayer wasperformed by Gerardot and coworkers , who also ob-served spin-layer locking of interlayer excitons in a 2H-WSe /MoSe trilayer . Tang et al. detected interac-tions between excitons and magnetically ordered holes inangle-aligned WSe /WS structures indicating that thissystem can be used to simulate the phase diagram ofthe triangular Hubbard model. The existence of stripephases over a large doping range has recently been re-ported in untwisted WSe /WS bilayers by Mak et al. .The same system also shows an abundance of correlatedinsulating states across a range of electron and hole dop-ing levels .To understand these experimental findings, detailedknowledge of the electronic structure of twisted TMDbilayers is required. Several groups have carried outdensity-functional theory (DFT) calculations of twistedhomobilayers. For example, Naik and Jain have cal-culated the band structure of several homobilayers (ne-glecting the effect of spin-orbit coupling) at a twist angleof 3 . etal. used the ab initio tight-binding model developedby Fang and coworkers for untwisted homobilayers andcalculated the band structure of MoS homobilayers fortwist angles as small as 1 . etal. used a Slater-Koster based tight-binding approachto study the evolution of flat bands in twisted homobi- a r X i v : . [ c ond - m a t . m t r l - s c i ] F e b layer MoS . As an alternative to atomistic methods, Wuand coworkers employed a continuum effective mass ap-proach to study the electronic structure of twisted heter-obilayers. Similar work was carried out by Zhang, Yuanand Fu .In this work, we systematically study the atomic andelectronic structure of all twisted homo- and heterobilay-ers that can be constructed by combining MoS , MoSe ,WS and WSe monolayers. Specifically, we use classicalforce fields to calculate the relaxed atomic structure ofthese systems, which display significant out-of-plane re-laxations. For the relaxed structures, we use an atomistictight-binding model derived from first-principles DFTcalculations to calculate the electronic band structure in-cluding the effect of spin-orbit coupling. In all homobi-layers, we find that the two highest valence bands becomeextremely flat as the twist angle approaches zero, reach-ing bandwidths of a few meV for twist angles near 2 ◦ . Incontrast, not all heterobilayers exhibit such flat valencebands. In particular, heterobilayers whose top valencebands derive from monolayer states at K and K’ exhibitdispersive top valence band states. These differences canbe understood by comparing the energy scale for inter-layer hopping with the energy difference between the va-lence band K- and Γ-states of the constituent monolay-ers. Importantly, the neglect of atomic relaxations leadsto qualitatively different electronic properties. II. METHODSA. Atomic structure
As a first step, we generate structures of flat (i.e. un-relaxed) twisted TMD homo- and heterobilayers (tBL-TMDs). We start from perfectly aligned bilayers, wheremetal and chalcogen atoms of the top layer are directlyabove corresponding metal and chalcogen atoms of thebottom layer, and rotate the top layer by an angle θ around the axis perpendicular to the plane of the bilayerand going through the metal atoms. For homobilayers, acommensurate structure is obtained when the moir´e cellvectors t and t can be expressed as t = n a + m a , t = − m a + ( n + m ) a , (1)where a = a ( √ , , a = a ( √ , − ,
0) are the prim-itive lattice vectors of the monolayer (with a being thelattice constant) and m and n are integers. The twistangle is given by cos θ = n +4 nm + m n + nm + m ) and the numberof atoms in the cell is N at = 6( n + nm + m ).For heterobilayers, we first consider systems whoseconstituent monolayers contain the same type of chalco-gen atoms. In this case, the lattice constants of bothmonolayers differ by less than 1% and we generate acommensurate moir´e cell for the twisted heterobilayersby increasing the lattice constant of the monolayer withthe smaller lattice constant to the value of the larger lattice constant and then use the same approach as forhomobilayers, described above.In contrast, for heterobilayers whose constitutentmonolayers contain different types of chalcogen atoms,the lattice constants of the monolayers differ by severalpercent. To generate moir´e cells for these systems, wefollow the approach of Zeller and G¨unther . In theirwork the moir´e vectors t and t are defined as t = n a + m a , t = − m a + ( n − m ) a , (2)where a = a (1 , , a = a ( − , √ ,
0) are the primitivelattice vectors of the monolayer with smaller equilibriumlattice vector a (with a denoting the larger lattice con-stant). We use DFT equilibrium lattice constants fromRef. . The integers n and m are determined from the nu-merical solution of a diophantine equation, see Appendixof Ref. . Here, we only consider so-called first-ordermoir´e structures . Importantly, to generate a commen-surate moir´e cell near a desired target twist angle a cer-tain level of strain must be applied. In this work, weonly 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 equilibriumlattice constant, i.e. a is reduced.For both homo- and hetero-bilayers, using the flattwisted bilayers as starting points, we determinethe relaxed equilibrium atomic structure via classi-cal force fields as implemented in the LAMMPS soft-ware package . In particular, we employ theforce fields developed by Naik and coworkers, basedon the Kolmogorov-Crespi potential for the interlayerinteraction . For intralayer interactions, a Stillinger-Weber type force field is used . The relaxed structuresof most twisted TMD bilayers are of a breathing-modetype, i.e. the two layers have out-of-plane displacementsin opposite directions (see Sec. III A). The only excep-tion are heterobilayers with different chalcogen atoms.These systems also exhibit breathing-mode relaxationsfor twist angles above 4 . ° degree, but for smaller twistangles qualitatively different relaxed structures are foundin which both layers move in the same direction andthe amplitude of the out-of-plane displacement is verylarge. In an experimental setting, we do not expect suchstructures to occur because the twisted TMD bilayersare placed on a substrate. Therefore, we do not presentresults for these systems. B. Electronic structure
To calculate the electronic properties of twisted TMDhomo- and heterobilayers, we use an atomistic tight-binding approach based on the work of Fang andcoworkers for untwisted homobilayers. The atomic or-bital basis consists of 5 d-orbitals for the metal atoms and3 p-orbitals for each chalcogen atom (which increases to10 d-orbitals and 6 p-orbitals, respectively, if spin-orbitcoupling is included). In a first step, we construct asymmetry-adapted tight-binding model for the monolay-ers including on-site, first, second and selected third near-est neighbor hoppings. The required hopping parametersare determined from a Wannier transformation of theDFT Hamiltonian. To model bilayers, Fang and cowork-ers describe interlayer hoppings between the p-orbitalsof the chalcogen atoms at the interface between the twolayers, which we refer to as “inner” chalcogens, using theSlater-Koster approach . The Slater-Koster parametersare fitted to a set of DFT calculations of untwisted bilay-ers in which the top layer is translated horizontally whilethe bottom layer is kept fixed. Finally, spin-orbit cou-pling is introduced via an on-site atomic term λ SOM / X L · S (with L and S denoting orbital and spin angular momen-tum operators, respectively, and λ SOM / X is the spin-orbitcoupling strength of M or X atoms).To model twisted homo- and heterobilayers, we haveextended the tight-binding model of Fang et al. in sev-eral ways. In particular, we have included interlayer hop-pings from inner chalcogen p z -orbitals on one layer tometal d z orbitals on the other layer using a Slater-Kosterapproach. Moreover, to better capture the effect of out-of-plane displacements present of the atoms, we improvethe description of interlayer hoppings (both p-p and p z -d z ) by using a different set of Slater-Koster parametersfor each value of the interlayer separation. All Slater-Koster parameters for the interlayer interactions as wellas all intralayer hoppings were obtained from a Wanniertransformation of the DFT Hamiltonian. For heter-obilayers, additional care must be taken to ensure thatthe on-site energies are referenced to the vacuum level.To determine the interlayer hoppings in a twisted bi-layer, the orbital basis of the rotated monolayer must betransformed. As described above, only p x , p y , p z and d z orbitals are involved in interlayer hoppings. Since d z andp z orbitals are unaffected by rotations around the z -axis,we only need to transform the p x and p y orbitals and therotated 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 and p y or-bitals transform in a similar fashion when a twist is in-troduced.Additional details about the interlayer tight-bindingmodel, the determination of the hopping parameters anda full list of the parameters for all systems can be found inthe Appendix A and in the Supplementary Information(Sec. S5). We have compared the band structures (with-out spin-orbit coupling) from this tight-binding modelto results from explicit DFT calculations for differenttwisted bilayers, see Sec. S2 in the Supplementary Infor-mation. Overall, we find very good agreement betweenthe two methods, in particular for the valence bands. III. RESULTSA. Atomic structure
Introducing a twist between two aligned TMD lay-ers results in the creation of a moir´e pattern consistingof regions with different stacking arrangements. High-symmetry stackings include AA regions, where the metal(chalcogen) atoms of the one layer are directly above themetal (chalcogen) atoms of the other layer, as well as twotypes of Bernal-like regions, where in one case (denotedB M / X ) the chalcogen atom (X) in one layer lies directlyabove the metal atom (M) in the other layer, or viceversa (denoted B X / M ). The high-symmetry stackings areshown in Fig. 1(a).
1. Homobilayers
All twisted homobilayers exhibit similar out-of-planeand in-plane displacement patterns upon relaxation.For example, Figs. 1(a)-(d) show results for twistedMoSe /MoSe at a twist angle of θ = 2 . ° . The inter-layer separation (ILS), defined as the distance betweenthe two surfaces on which the inner chalcogen atoms lie,is large in the AA regions (which form a triangular lat-tice), 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 alongthe diagonal of the moir´e unit cell for three different twistangles. As the twist angle decreases, the size of the AAregions shrinks, whereas B M / X and B X / M regions expand.This allows the system to reduce its energy as AA regionsare energetically unfavorable because of their large stericrepulsion. It can further be observed that the maximumILS increases, while the minimum ILS decreases as thetwist angle is reduced. Again, this reduces the energycost associated with steric repulsion.Figures 1(c) and (d) show the in-plane displacementsof twisted MoSe /MoSe . Similar to twisted bilayergraphene, the in-plane displacements in tBL-TMDs formvortices around the AA regions, with the atoms in the topand bottom layers rotating in opposite directions. Atomsin the B M / X and B X / M regions are almost unaffected byin-plane relaxations. The magnitude of in-plane atomicdisplacements around AA regions increases for small an-gles. This allows the system to reduce the size of theenergetically unfavourable AA regions.
2. Heterobilayers
Figure 2 shows the in-plane and out-of-plane re-laxations of twisted MoSe /MoS at a twist angle of θ = 4 . ° . For the set of angles studied in this work, wefind that heterobilayers exhibit similar relaxation pat-terns as homobilayers: large ILSs are found in the AAregions, which form a triangular lattice. The relative (a) B M/X B X/M (b)(c) MoSe top layer (d) MoSe bottom layer FIG. 1: Atomic relaxations in twisted MoSe /MoSe homobilayer. (a) Left: Inter-layer separation (ILS), defined asthe distance between the two surfaces generated by the inner chalcogen atoms, for θ = 2 . ° . Right: atomic stackingarrangements in AA, B M / X and B X / M regions of the moir´e cell. (b) Out-of-plane displacement ∆ z along thediagonal of a 2 × s = α ( t + t ) with α ranging from 0 to 2, shown as green dashed line in (a), forthree different twist angles. (c) and (d) show the in-plane displacements | ∆ u | of the top and bottom layer,respectively, for θ = 2 . ° . Arrows indicate the direction of the in-plane displacements with their lengths beingproportional to the magnitude of the atomic displacements.size of the AA regions shrinks as the twist angle is de-creased while B M / X and B X / M regions grow. In contrastto the homobilayers, the in-plane and out-of-plane dis-placements of the two layers are not symmetric, as canbe seen in Figs. 2(b),(c) and (d). The difference of theout-of-plane displacements in the AA and B M / X regionsis about four times larger in the MoS layer than in theMoSe layer (Fig. 2(b)). As we show in the next sec-tion, this asymmetry is less pronounced in heterobilayersthat have the same chalcogens. Similar to the out-of-plane displacements, the in-plane displacements are alsolarger in the MoS layer compared to the MoSe layer(Figs. 2(c) and 2(d)).
3. Chemical trends
Figure 3(a) shows the maximum and minimum ILS,corresponding to the ILS value in the center of the AAand B M / X / B X / M regions, respectively, as function oftwist angle for the entire set of TMD homo- and hetero-bilayers. At large twist angles, the ILSs in these regionsdiffer significantly from the values in the untwisted AAand B M / X bilayers. The ILS in the B M / X /B X / M regions(bottom panel of Fig. 3(a)) decreases monotonically asthe twist angle is reduced and converges to the ILS ofthe untwisted bilayers. In contrast, the ILS in the AAregions (top panel of Fig. 3(a)) increases with decreasingtwist angle, but does not converge to the value of the un-twisted bilayer in the case of homobilayers and heterobi-layers with same chalcogen atoms. This discontinuity ofthe maximum ILS at θ = 0 ◦ is a consequence of the struc- (a) (b)(c) MoSe layer (d) MoS layer FIG. 2: Atomic relaxations in twisted MoSe /MoS heterobilayer. (a) Inter-layer separation (ILS), defined as thedistance between the two surfaces generated by the inner chalcogen atoms, for θ = 4 . ° . Colored dots refer todifferent stacking regions as described in Fig. 1(a). (b) Out-of-plane displacement ∆ z along the diagonal of a 2 × s = α ( t − t ) with α ranging from 0 to 2, shown as green dashed line in (a), for three different twistangles. (cid:15) denotes the compressive strain in the MoSe layer which is needed to make the lattice constants of bothlayers equal and generate commensurate moir`e cells. (c) and (d) show the in-plane displacements | ∆ u | of the topand bottom layer, respectively, for θ = 4 . ° . Arrows indicate the direction of the in-plane displacements with theirlengths being proportional to the magnitude of the atomic displacements.tural relaxations which result in a growth of the B M / X and B X / M regions and a shrinkage of the AA regions atsmall twist angles. At the center of the large B M / X /B X / M regions the twisted bilayer has a similar structure as theuntwisted B M / X /B X / M bilayer while the small size of theAA restricts the atoms from reaching the same structureas the untwisted AA bilayer.Comparing the ILS of different bilayers, we ob-serve that bilayers where both constituent monolay-ers contain S atoms (MoS /MoS , WS /WS andWS /MoS ) exhibit the smallest interlayer distances(both in AA and B M / X /B X / M regions), whereas bilayerscontaining Se atoms (MoSe /MoSe , WSe /WSe andWSe /MoSe ) in both layers exhibit the largest ILSs.Bilayers with S atoms in one layer and Se atoms inthe other (MoSe /MoS , MoSe /WS , WSe /MoS andWSe /WS ) show intermediate values of the ILS. These trends can be explained by the different van der Waalsradii of S and Se atoms, which are ∼ . ∼ . .The out-of-plane displacements for all homobilay-ers and heterobilayers with same chalcogen species at θ = 5 . ° and for all heterobilayers with different chalco-gens at θ = 5 . ° are shown in Fig. 3(b). As shown in theleft panel of Fig. 3(b), out-of-plane displacements in ho-mobilayers 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, theamplitude of the displacement pattern of the bottom lay-ers (which are unstrained) is similar to that found in thehomobilayers, while the amplitudes of the strained toplayer are somewhat smaller. (a)(b) FIG. 3: (a) Minimum (bottom panel) and maximum (top panel) interlayer separation (ILS), corresponding to B M / X and AA regions, respectively, for all tBL-TMDs as function of twist angle. For heterobilayers with differentchalcogen atoms three twist angles are shown and these are 4 . ° , 5 . ° and 7 . ° . For reference, the ILS of untwistedAA bilayers (top panel) and B M / X bilayers (bottom panel) are also shown by short horizontal lines on the left handside of the plots. (b) Out-of-plane displacements ∆ z for all homobilayers (left) and all heterobilayers (right) alongthe diagonal of a 2 × s = α ( t + t ) for homobilayers, and s = α ( t − t ) for heterobilayers, α ranges from 0 to 2). For all homobilayers and heterobilayers with same chalcogen atoms the twist angle is 5 . ° ;for all heterobilayers with different chalcogens θ = 5 . ° . B. Electronic structure
1. Homobilayers
In this section we study the evolution of the bandstructure of the twisted homobilayers MoS /MoS , MoSe /MoSe , WS /WS and WSe /WSe as functionof twist angle. All calculations were carried out for therelaxed structures and include spin-orbit coupling. Forall twist angles, the homobilayers exhibit a semiconduct- (a) KK ’ K M KM Gq K ’ (b) Γ M K Γ WSe /WSe -70-60-50-40-30-20-10010 E ne r g y ( m e V ) DOS VB1VB2 arbitrary units (c)
VB1 at WSeSeWSeSe
VB1VB2 Å - Å - Å - Å - Å - Å - FIG. 4: Electronic structure of twisted WSe /WSe . (a) Moir´e Brillouin zone (black hexagon) obtained by twistingtwo monolayers (whose Brillouin zones are indicated by blue and red hexagons) by θ = 22 ° . The Γ-M-K-Γ pathused for computing band structures is also shown (yellow line). (b) Left: Band structure for θ = 3 . ° near thevalence band edge. The two highest valence bands (denoted VB1 and VB2) are shown in red. Right: Density ofstates per WSe formula unit. (c) Layer-resolved | ψ Γ ( r ) | of VB1. (d) Layer-averaged squared wavefunctions of VB1(top panels) and VB2 (bottom panels) at Γ, M and K. Colored dots refer to different stacking regions as describedin Fig. 1(a) and the moir´e unit cell is indicated by grey lines.ing band structure with a band gap separating the va-lence and conduction bands. Moreover, the two highestvalence bands (each of which is spin degenerate) are sep-arated from all other “remote” valence bands by energygaps when θ < ◦ . We refer to these two highest valencebands as VB1 and VB2, respectively.As an example, we focus on the valence band structureof WSe /WSe . Fig. 4 (b) shows the band structure at θ = 3 . ° . It can be seen that the valence band maximumoccurs at the Γ-point of the moir´e Brillouin zone and thatVB1 and VB2 touch at the K-point forming a Dirac cone.This is also reflected in the V-shaped density of states.The graphene-like dispersion of VB1 and VB2 can be un-derstood by analyzing the wavefunctions of these states.Fig. 4(d) shows that the wavefunctions at Γ are localizedin the B M / X and B X / M regions which form a honeycomb (a) MoS /MoS FLAT -0.2-0.1 0 Γ M K Γ E ne r g y ( e V ) E ne r g y ( e V ) (b) MoS /MoS RELAXED -0.2-0.1 0 Γ M K Γ E ne r g y ( e V ) E ne r g y ( e V ) (c) (d) FIG. 5: Band structure of twisted MoS /MoS at θ = 5 . ° using (a) the unrelaxed flat atomic structure and (b) therelaxed atomic structure. Layer-averaged | ψ Γ ( r ) | of the highest valence band for (c) flat and (d) relaxed structures.Colored dots refer to different stacking regions as described in Fig. 1(a).lattice. Importantly, the total bandwidth of the two high-est valence bands is less than 30 meV demonstrating theformation of flat bands upon twisting.To understand the chemical origin of the flat bands, weanalyze their projections onto atomic orbitals. Fig. 4(c)shows that these states are localized symmetrically onthe inner layers of chalcogen atoms and also on the twometal layers. These states are mostly composed of innerchalcogen p z orbitals and metal d z orbitals, as show inFig. S5 of the Supplementary Information. This suggeststhat VB1 and VB2 originate from Γ-states of the con-stituent monolayers: in all TMD monolayers, the top va-lence band states at Γ have large projections onto chalco-gen p z orbitals and metal d z orbitals, whereas the topvalence band states at K and K’ have large projectionsonto metal d xy and d x − y orbitals .When atomic relaxations are not taken into account, aqualitatively different valence band structure is obtained.Figs. 5(a) and 5(b) compare the band structures of flatunrelaxed and relaxed MoS /MoS at θ = 5 . ◦ . In con-trast to the relaxed system, the flat system does not ex-hibit a Dirac-like dispersion and exhibits an energy gap between VB1 and VB2. This is a consequence of thedifferent spatial structure structure of the correspondingwavefunction, see Figs. 5(c) and 5(d): in the relaxedsystem VB1 and VB2 localize in the B M / X and B X / M re-gions, while in the flat system the top valence band statesare localized in the AA regions which form a triangularlattice . We find a similar effect of atomic relaxations inall homobilayers, see Supplementary Information (S1).Figure 6 compares the band structures of all homobi-layers at three twist angles ( θ = 5 . ° , 2 . ° and 1 . ° ). At θ = 5 . ° , the top two valence bands are not fully isolatedfrom the remote valence bands and significant differencesin the dispersion of these remote bands can be observedin the different materials.When the twist angle is reduced to θ = 2 . ° , the bandsbecome flatter and the top two valence bands becomeisolated from the remote valence bands. Decreasing thetwist angle further to θ = 1 . ° , the highest four remote va-lence bands also become isolated in energy from all otherremote bands, see Figs. 6(c),(f),(i),(l). The two middlebands of this set also exhibit a Dirac-like dispersion nearthe K-point, while the highest and lowest bands are veryflat. Xian and coworkers have analyzed the origin ofthese bands and proposed that can be described by a setof p x and p y orbitals on a honeycomb lattice.Figure 7(a) shows the bandwidth w , computed as theenergy difference between states at Γ and K, of the topvalence band (denoted as VB1 as in Fig. 4(b)) as func-tion of twist angle for relaxed and unrelaxed (flat) homo-bilayers. As the twist angle decreases, the bandwidthsapproach zero. For relaxed homobilayers, the magnitudeof the bandwidths in the different systems are relativelysimilar with MoSe /MoSe exhibiting the smallest one(reaching ≈ θ = 1 . ° ). When relaxationsare neglected, the bandwidths are smaller. For exam-ple, a bandwidth of only 1 meV is found in unrelaxedMoSe /MoSe at θ = 2 . ° . In addition, when relaxationsare neglected a different ordering of the bandwidths isobtained, with twisted WSe /WSe showing the largestbandwidth over the entire range of twist angles investi-gated here.Figure 7(b) shows the band gap ∆ as function oftwist angle for both relaxed and unrelaxed homobilayers.For all relaxed homobilayers ∆ decreases linearly with aslope of ≈
20 meV/degree as the twist angle is reduced.WS /WS has the largest band gap (1 . − . /MoSe the smallest (1 . − . M / X (or B X / M ) bilayers, shown on theleft panel of Fig. 7(b) at θ = 0 ° . That is expected asthe B M / X and B X / M regions are energetically favorable(compared to AA regions) and their relative size growsas the twist angle is reduced (see Sec. II A).Interestingly, the nature of the band gap ofWSe /WSe changes from direct (K → K) to indirect(Γ → K) at θ = 7 . ° . This is a consequence of the signif-icant mixing of monolayer states from Γ and K/K’ thatoccurs at large twist angles in this system. All othersystems exhibit indirect band gaps. In particular, forMoS /MoS and MoSe /MoSe the valence band maxi-mum is at Γ and the conduction band minimum at K,while for WS /WS the conduction band minimum ishalf-way between the Γ-point and the M-point (referredto as the X-point).Without relaxations, the band gaps are almost con-stant and do not depend sensitively on the twist angle,see Fig. 7(b). Also, the nature of the band gap is differentcompared to the relaxed systems.
2. Heterobilayers
We first consider heterobilayers with the samechalcogen species in each layer, i.e., WS /MoS andWSe /MoSe . As discussed in Sec. II, it is possible togenerate commensurate moir´e structures with very littlestrain for these systems. The band structures of thesesystems at three different twist angles are shown in Fig. 8. The flattening of the top valence bands can be clearlyobserved. Similar to the homobilayers at small twistangles, the top two valence bands are separated fromall other remote valence bands and have a small band-width (Figs. 8(b)-(c) and Figs. 8(e)-(f), for WS /MoS and WSe /MoSe , respectively). These bands also havelarge projections onto inner chalcogen p z and metal d z orbitals. However, the projections are no longer layer-symmetric: the projections onto the WS layer are largerthan those of the MoS layer (see Fig. S6 of Supplemen-tray Information).In contrast to the homobilayers, the top valence bandsVB1 and VB2 are separated by an energy gap at K. Thisgap is larger for WS /MoS : at θ = 1 . ° its value is5 meV, while it is only 1 meV for WSe /MoSe . Simi-larly, additional energy gaps are introduced into the high-est set of remote valence bands. Note that the energygaps in these heterobilayer systems are not caused byspin-orbit coupling. but by the ”chemical asymmetry”of the two constituent layers.Interestingly, WS /MoS exhibits a qualitatively dif-ferent band structure at θ = 5 . ° . At this twist angle, adispersive band crosses the flat band and changes the po-sition of valence band maximum from K to Γ. The sameband can be also be seen in WSe /MoSe at θ = 5 . ° ,but in this system the dispersive band has a slightlylower energy than the flat top valence band. We deferthe discussion about the origin of these dispersive bandsto Sec. III C.Next, we first consider the heterobilayers MoSe /MoS and MoSe /WS , i.e. hetero-bilayers containing differentspecies of chalcogens but no WSe . Fig. 9(a)-(d) showsthe band structures of these systems at two twist angles( θ = 7 . ° and θ = 4 . ° ). For MoSe /MoS the top va-lence bands have large projections onto p z orbitals of theinner chalcogen atoms and metal d z orbitals, and theprojections are layer-asymmetric with significantly moreweight on the MoSe layer, see Fig. S7 of SupplementaryInformation. In contrast to the homobilayers and the het-erobilayers containing a single chalcogen species, the topvalence bands in these systems are not spin-degenerate.Fig. 10(a) shows the spin-resolved dispersion of top va-lence bands in MoSe /MoS at θ = 4 . ° which exhibitsspin splittings with magnitudes up to 13 meV. Also, anenergy gap between VB1 and VB2 of 8 meV is found atthe K-point. Interestingly, these bands are partially spin-polarized: it can be observed that the top valence bandsare only spin-polarized in the vicinity of the K-point eventhough spin splitting occurs along the whole bandstruc-ture path. This scenario was recently discussed by Liuand coworkers who demonstrated that spin-orbit cou-pling can lead to spin splittings without spin polarizationin non-magnetic materials without inversion symmetry.Finally, Figs. 9(e)-(h) show the band structures oftwisted WSe /MoS and WSe /WS , i.e. the heterobi-layers with different chalcogen atoms that contain WSe .The top valence bands in these systems are dispersiveand the flat bands are observed at lower energies, see0 -0.2-0.1 0 G M K G E ne r g y ( e V ) E ne r g y ( e V ) -0.2-0.1 0 G M K G E ne r g y ( e V ) E ne r g y ( e V ) -0.2-0.1 0 G M K G E ne r g y ( e V ) E ne r g y ( e V ) -0.2-0.1 0 G M K G E ne r g y ( e V ) E ne r g y ( e V ) MoS /MoS MoSe /MoSe WS /WS WSe /WSe q = q = q =
MoS MoSe RELAXED θ (°) ∆ ( e V ) FLAT θ (°) ∆ ( e V ) WS WSe RELAXED : G K, : G M,: G XIndirect : K KDirect
FIG. 7: (a) Band width w and (b) band gap ∆ of the top valence band (VB1) as function of twist angle θ forrelaxed and unrelaxed (flat) homobilayers. For relaxed WSe /WSe , isolated bands at the valence band edge appearonly for twist angles θ ≤ . ° .roughly the same for the different homo- and heterobi-layers.The band gap of the twisted heterobilayers does not de-pend sensitively on twist angle, see Fig. 11(b). For mostsystems, a mild reduction of the gap can be observedas the twist angle decreases. The smallest bands gaps( ≈ . /MoS , while MoSe /WS and WS /MoS exhibit the largest gaps ( ≈ . − . /MoS , and from indirect(Γ-X) to direct (Γ-Γ) in MoSe /WS . C. Physical origin of the dispersive valence bands
We now focus on the set of heterobilayers that exhibitdispersive bands at the valence band edge and discussthe origin of these bands. As mentioned in Sec. III B 2,such bands are observed in WS /MoS at θ = 5 . ° andalso in WSe /MoS and WSe /WS , see Figs. 9(e)-(h).Compared to the flat bands, the width of these bands de-creases much less as the twist angle is reduced. Project-ing the corresponding states onto atomic orbitals revealslarge contributions from W d xy and d x − y orbitals, seeFig. S8 of Supplementary Information, suggesting thatthey originate from states at the K-point of WS mono-layer in the case of WS /MoS and similarly from K-point states of the WSe monolayer in the WSe /MoS and WSe /WS systems . As these states have verysmall projections onto the inner chalcogen atoms (seeFig. S9 of the Supplementary Information), they are only weakly affected by interlayer coupling which explains whyintroduction of a twist does not result in a significant re-duction of their band width.The top valence band at K of the WSe and WS monolayers is spin-polarized and the dispersive valencebands of the twisted bilayers inherit this property asshown for WSe /MoS in Fig. 10(b). The valence bandmaximum in bilayers with dispersive bands is located atthe K -point. In the WS /MoS bilayer the band gapis direct, whereas it is indirect in the WSe /MoS andWSe /WS bilayers as the conduction band minimum isat Γ.Closer inspection of the band structure of WSe /MoS and WSe /WS (see, for example, Fig. 10(b)) shows thatthese systems also feature flat bands, but with lower en-ergies than the dispersive bands. We have found that theordering of the flat and dispersive valence bands dependssensitively on the atomic structure and the twist angle.For example, Fig. 12(a) shows that neglecting atomic re-laxations results in a different ordering of dispersive andflat valence bands. Moreover, we have found that it ispossible to switch the order of flat and dispersive bandswhen the interlayer separation of the relaxed structuresis reduced suggesting that the electronic properties ofthese 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 va-lence bands in the different heterobilayers, we proposea simple model in which the bilayer states originatingfrom K / K - or Γ-valleys of the constituent monolayersare obtained from a two-level system. Specifically, the2 -0.2-0.1 0 G M K G E ne r g y ( e V ) E ne r g y ( e V ) -0.2-0.1 0 G M K G E ne r g y ( e V ) E ne r g y ( e V ) q = q = q =
We have studied the electronic band structures of alltwisted transition metal dichalcogenide (TMD) bilayersthat can be obtained by combining MoS , MoSe , WS and WSe monolayers. Specifically, we have carried outtight-binding calculations taking into account the effectof atomic relaxations and also spin-orbit coupling. In alltwisted homobilayers, the top valence bands are derivedfrom monolayer states at Γ and become flat when thetwist angle decreases. For twisted heterobilayers, we find3 -0.2-0.1 0 G M K G E ne r g y ( e V ) E ne r g y ( e V ) -0.1 0 G M K G E ne r g y ( e V ) E ne r g y ( e V ) -0.3-0.2-0.1 0 G M K G E ne r g y ( e V ) E ne r g y ( e V ) -0.1 0 G M K G E ne r g y ( e V ) E ne r g y ( e V ) -0.2-0.1 0 G M K G E ne r g y ( e V ) E ne r g y ( e V ) -0.2-0.1 0 G M K G E ne r g y ( e V ) E ne r g y ( e V ) -0.2-0.1 0 G M K G E ne r g y ( e V ) E ne r g y ( e V ) -0.2-0.1 0 G M K G E ne r g y ( e V ) E ne r g y ( e V )
FIG. 11: (a) Bandwidth w and (b) band gap ∆ of the top valence band for different heterobilayers as a function oftwist angle θ . We only show results for twisted bilayers with flat valence bands (see discussion in the text).bilayers can be harnessed to design flat band propertiesand pave the way to understanding electron-electron in-teraction effects in these materials. ACKNOWLEDGMENTS
Appendix A: Improvements on the Tight-bindingModel
In this Appendix, we describe the modifications thatwere required to generalize the tight-binding modelfor untwisted TMD bilayers developed by Fang andcoworkers to twisted TMD bilayers. In particu-lar, we improved the description of interlayer hoppingsand parametrized these hoppings for all combinationsof homo- and heterobilayers composed of MoS /MoS ,MoSe /MoSe , WS /WS and WSe /WSe monolayers.5 (a) WSe /MoSe FLAT -0.2-0.1 0 Γ M K Γ E ne r g y ( e V ) E ne r g y ( e V ) (b) WSe /MoSe RELAXED -0.2-0.1 0 Γ M K Γ E ne r g y ( e V ) E ne r g y ( e V ) FIG. 12: Band structure of twisted WSe /MoSe at θ = 5 . ° using (a) the unrelaxed flat atomic structure and (b)relaxed atomic structure.TABLE I: Parameters of H Γ and H K , see Eq. 3. ε (1)Γ(K) denotes the energy at Γ (K) of the monolayer with thehigher-lying valence band, and ε (2)Γ(K) denotes the energy at Γ (K) of the monolayer with the lower-lying valenceband. ∆ Γ is the interlayer coupling. ε maxΓ denotes the largest eigenvalue of H Γ . All values are in eV. ε (1)Γ ε (2)Γ ε (1)K ∆ Γ ε maxΓ ε maxΓ − ε (1) K MoS /MoS − . − . − .
97 2 . − .
51 0 . /MoSe − . − . − .
25 1 .
83 0 .
12 0 . /WS − . − . − .
50 2 . − .
34 0 . /WSe − . − . − .
01 1 .
79 0 .
11 0 . /MoS − . − . − .
50 1 . − .
45 0 . /MoSe − . − . − .
01 1 .
08 0 .
10 0 . /MoS − . − . − .
25 0 . − .
21 0 . /WS − . − . − .
25 1 . − .
10 0 . /MoS − . − . − .
01 0 . − . − . /WS − . − . − .
01 1 . − . − . z -orbitals and transition metal d z -orbitals, but did not describe this with a Slater-Koster expression. Instead, they only calculated thevalue of this hopping for the specific geometry of an un-twisted 2H bilayer.To generalize the description of p z to d z hoppings totwisted bilayers, we used the Slater-Koster formula t p z ,d z ( r ) = n (cid:20) n −
12 ( l + m ) (cid:21) V pdσ ( r ) (A1)+ √ n ( l + m ) V pdπ ( r ) , (A2)where the directional cosines are defined as l = r x /r , n = r y /r and n = r z /r .To determine the functions V pdσ ( r ) and V pdπ ( r ), wecalculated t p z ,d z and also t p z ,d xz and t p z ,d yz for a set of untwisted bilayers with different stacking configura-tions and different interlayer separations using a Wan-nier transformation of the DFT Hamiltonian, see Fig. 13.Next, a least square fitting process was used to extract V pdσ and V pdπ at different interatomic distances and theresults were fitted to functions of the type V pd ( r ) = V (cid:16) rd (cid:17) α cos (cid:16) β rd + γ (cid:17) (A3)with V, α, β and γ denoting fitting parameters, and d =3 . t p z ,d z hopping parameters. We have also testedthe influence of other hoppings between chalcogen p-orbitals and transition metal d-orbitals, but found thatthe most important contribution arises from p z to d z hoppings.Figure 14 compares the band structures of untwistedand twisted bilayer MoS at a twist angle of 7 . ° fromtight-binding with and without the p z to d z hopping6FIG. 13: Comparison of interlayer p z - d z hoppingparameters obtained from a Wannierization of the DFTHamiltonian to Slater-Koster (SK) model for untwistedMoS bilayers..to a first-principles density-functional theory result. Wefind that inclusion of p z to d z hoppings improves theagreement with the first principles result significantly. Inparticular, the valence band states near Γ are pushed tohigher energies which reduces the band gap by approx-imately 200 meV. A similar shift of the highest valencebands is found also in the twisted bilayers.Aside from the inclusion of interlayer p z to d z hop-pings, we also discovered that the description of inter-layer hoppings between chalcogen p-orbitals developed byFang and coworkers required improvements to obtain anaccurate description of twisted bilayers. To parametrizethe corresponding Slater-Koster expressions, Fang et al.carried out first-principles calculations of untwisted bi-layers with different stacking configurations, but usinga fixed interlayer separation. As we have demonstrated in the main section of the manuscript, the introductionof a twist results in significant atomic relaxations andconcomitant variations in the interlayer separation. Wehave found that these changes in the interlayer separationare not well captured by the simple Slater-Koster expres-sions used by Fang and coworkers. Fig. 15(a) shows theSlater-Koster parameters V ppσ and V ppπ as function ofthe interatomic distance for different interlayer separa-tions: for small interatomic distances, the Slater-Kosterparameters depend sensitively on the interlayer separa-tion.To account for the dependence of the Slater-Koster pa-rameters on the interlayer separation, the following pro-cedure is used: for a given pair of atoms, we first cal-culate the interlayer separation as the difference of theirz-coordinates as well as the interatomic distance. Then,the Slater-Koster parameters for the specific interlayerseparation are used to obtain the desired hopping matrixelement.The DFT calculations for monolayers and untwistedbilayers were performed with Quantum Espresso withinthe optB88 generalized gradient approximation for theexchange-correlation potential and a plane-wave cutoffvalue of 70 Ry ( ≈
950 eV). Monolayer calculations wereperformed with a 25 × × × × . Onetep ,a linear-scaling DFT code. We use the Perdew-Burke-Ernzerhof exchange-correlation functional withprojector-augmented-wave pseudopotentials , gener-ated from ultra-soft pseudopotentials , and a kineticenergy cutoff of 800 eV. A basis consisting of 9 non-orthogonal generalized Wannier functions for both metaland chalcogen atoms is employed. Yuan Cao, Valla Fatemi, Shiang Fang, Kenji Watanabe,Takashi Taniguchi, Efthimios Kaxiras, and Pablo Jarillo-Herrero. Unconventional superconductivity in magic-anglegraphene superlattices.
Nature , 556(7699):43–50, Apr2018. Yuan Cao, Valla Fatemi, Ahmet Demir, Shiang Fang,Spencer L. Tomarken, Jason Y. Luo, Javier D. Sanchez-Yamagishi, Kenji Watanabe, Takashi Taniguchi, EfthimiosKaxiras, Ray C. Ashoori, and Pablo Jarillo-Herrero. Cor-related insulator behaviour at half-filling in magic-anglegraphene superlattices.
Nature , 556(7699):80–84, Apr2018. Simone Lisi, Xiaobo Lu, Tjerk Benschop, Tobias A.de Jong, Petr Stepanov, Jose R. Duran, Florian Margot,Ir`ene Cucchi, Edoardo Cappelli, Andrew Hunter, AnnaTamai, Viktor Kandyba, Alessio Giampietri, Alexei Bari-nov, Johannes Jobst, Vincent Stalman, Maarten Leeuwen- hoek, Kenji Watanabe, Takashi Taniguchi, Louk Rade-maker, Sense Jan van der Molen, Milan P. Allan, Dmitri K.Efetov, and Felix Baumberger. Observation of flat bandsin twisted bilayer graphene.
Nature Physics , Sep 2020. Matthew Yankowitz, Shaowen Chen, Hryhoriy Polshyn,Yuxuan Zhang, K. Watanabe, T. Taniguchi, David Graf,Andrea F. Young, and Cory R. Dean. Tuning su-perconductivity in twisted bilayer graphene.
Science ,363(6431):1059–1064, 2019. Stephen Carr, Daniel Massatt, Shiang Fang, Paul Cazeaux,Mitchell Luskin, and Efthimios Kaxiras. Twistronics: Ma-nipulating the electronic properties of two-dimensional lay-ered structures through their twist angle.
Phys. Rev. B ,95:075420, Feb 2017. Nicolas Mounet, Marco Gibertini, Philippe Schwaller, Da-vide Campi, Andrius Merkys, Antimo Marrazzo, ThibaultSohier, Ivano Eligio Castelli, Andrea Cepellotti, Giovanni M K E n e r g y ( e V ) (a) DFTwith p z - d z w/out p z - d z (b) with p z - d z w/out p z - d z M K (c)
DFT
M K
FIG. 14: Comparison of the tight-binding band structure with and without p z - d z hopping to the first-principlesDFT result for (a) untwisted AA-stacked bilayer MoS and (b) twisted bilayer MoS at θ = 7 . ° . (c) CorrespondingDFT bandstructure. Pizzi, and Nicola Marzari. Two-dimensional materialsfrom high-throughput computational exfoliation of ex-perimentally known compounds.
Nature Nanotechnology ,13(3):246–252, Mar 2018. Thomas Mueller and Ermin Malic. Exciton physics anddevice application of two-dimensional transition metaldichalcogenide semiconductors. npj 2D Materials and Ap-plications , 2(1):29, Sep 2018. Wonbong Choi, Nitin Choudhary, Gang Hee Han, JuhongPark, Deji Akinwande, and Young Hee Lee. Recent de-velopment of two-dimensional transition metal dichalco-genides and their applications.
Materials Today , 20(3):116– 130, 2017. Lei Wang, En-Min Shih, Augusto Ghiotto, Lede Xian,Daniel A. Rhodes, Cheng Tan, Martin Claassen, Dante M.Kennes, Yusong Bai, Bumho Kim, Kenji Watanabe,Takashi Taniguchi, Xiaoyang Zhu, James Hone, Angel Ru-bio, Abhay N. Pasupathy, and Cory R. Dean. Corre-lated electronic phases in twisted bilayer transition metaldichalcogenides.
Nature Materials , 19(8):861–866, Aug2020. Liheng An, Xiangbin Cai, Ding Pei, Meizhen Huang, Ze-fei Wu, Zishu Zhou, Jiangxiazi Lin, Zhehan Ying, ZiqingYe, Xuemeng Feng, Ruiyan Gao, Cephise Cacho, MatthewWatson, Yulin Chen, and Ning Wang. Interaction effectsand superconductivity signatures in twisted double-bilayerwse2.
Nanoscale Horiz. , 5:1309–1316, 2020. Xiao-Guang Gao, Xiao-Kuan Li, Wei Xin, Xu-Dong Chen,Zhi-Bo Liu, and Jian-Guo Tian. Fabrication, optical prop-erties, and applications of twisted two-dimensional mate-rials.
Nanophotonics , 9(7):1717 – 1742, 01 Jul. 2020. Kha Tran, Galan Moody, Fengcheng Wu, Xiaobo Lu,Junho Choi, Kyounghwan Kim, Amritesh Rai, Daniel A.Sanchez, Jiamin Quan, Akshay Singh, Jacob Embley,Andr´e Zepeda, Marshall Campbell, Travis Autry, TakashiTaniguchi, Kenji Watanabe, Nanshu Lu, Sanjay K. Baner-jee, Kevin L. Silverman, Suenne Kim, Emanuel Tutuc, Li Yang, Allan H. MacDonald, and Xiaoqin Li. Evidencefor moir´e excitons in van der waals heterostructures.
Na-ture , 567(7746):71–75, Mar 2019. Junho Choi, Wei-Ting Hsu, Li-Syuan Lu, Liuyang Sun,Hui-Yu Cheng, Ming-Hao Lee, Jiamin Quan, Kha Tran,Chun-Yuan Wang, Matthew Staab, Kayleigh Jones,Takashi Taniguchi, Kenji Watanabe, Ming-Wen Chu,Shangjr Gwo, Suenne Kim, Chih-Kang Shih, Xiaoqin Li,and Wen-Hao Chang. Moir´e potential impedes interlayerexciton diffusion in van der waals heterostructures.
ScienceAdvances , 6(39), 2020. Malte Kremser, Mauro Brotons-Gisbert, JohannesKn¨orzer, Janine G¨uckelhorn, Moritz Meyer, Matteo Bar-bone, Andreas V. Stier, Brian D. Gerardot, Kai M¨uller,and Jonathan J. Finley. Discrete interactions between afew interlayer excitons trapped at a mose2–wse2 heteroin-terface. npj 2D Materials and Applications , 4(1):8, May2020. Mauro Brotons-Gisbert, Hyeonjun Baek, AlejandroMolina-S´anchez, Aidan Campbell, Eleanor Scerri, DanielWhite, Kenji Watanabe, Takashi Taniguchi, Cristian Bon-ato, and Brian D. Gerardot. Spin–layer locking of inter-layer excitons trapped in moir´e potentials.
Nature Materi-als , 19(6):630–636, Jun 2020. Yanhao Tang, Lizhong Li, Tingxin Li, Yang Xu, Song Liu,Katayun Barmak, Kenji Watanabe, Takashi Taniguchi, Al-lan H. MacDonald, Jie Shan, and Kin Fai Mak. Simulationof hubbard model physics in wse2/ws2 moir´e superlattices.
Nature , 579(7799):353–358, Mar 2020. Chenhao Jin, Zui Tao, Tingxin Li, Yang Xu, YanhaoTang, Jiacheng Zhu, Song Liu, Kenji Watanabe, TakashiTaniguchi, James C. Hone, Liang Fu, Jie Shan, and Kin FaiMak. Stripe phases in wse2/ws2 moir´e superlattices, 2020. Yang Xu, Song Liu, Daniel A. Rhodes, Kenji Watanabe,Takashi Taniguchi, James Hone, Veit Elser, Kin Fai Mak,and Jie Shan. Correlated insulating states at fractionalfillings of moir´e superlattices.
Nature , 587(7833):214–218, Nov 2020. Mit H. Naik and Manish Jain. Ultraflatbands and shearsolitons in moir´e patterns of twisted bilayer transitionmetal dichalcogenides.
Phys. Rev. Lett. , 121:266401, Dec2018. Zhen Zhan, Yipei Zhang, Pengfei Lv, Hongxia Zhong,Guodong Yu, Francisco Guinea, Jose Angel Silva-Guillen,and Shengjun Yuan. Multi-ultraflatbands tunability andeffect of spin-orbit coupling in twisted bilayer transitionmetal dichalcogenides, 2020. Shiang Fang, Rodrick Kuate Defo, Sharmila N. Shirodkar,Simon Lieu, Georgios A. Tritsaris, and Efthimios Kaxiras.Ab initio tight-binding hamiltonian for transition metaldichalcogenides.
Phys. Rev. B , 92:205108, Nov 2015. Somepalli Venkateswarlu, Andreas Honecker, and GuyTrambly de Laissardi`ere. Electronic localization in twistedbilayer mos with small rotation angle. Phys. Rev. B ,102:081103, Aug 2020. Fengcheng Wu, Timothy Lovorn, Emanuel Tutuc, andA. H. MacDonald. Hubbard model physics in transi-tion metal dichalcogenide moir´e bands.
Phys. Rev. Lett. ,121:026402, Jul 2018. Yang Zhang, Noah F. Q. Yuan, and Liang Fu. Moir´equantum chemistry: charge transfer in transition metaldichalcogenide superlattices, 2020. Patrick Zeller and Sebastian G ˜AŒnther. What are thepossible moir´e patterns of graphene on hexagonally packedsurfaces? universal solution for hexagonal coincidence lat-tices, derived by a geometric construction.
New Journal ofPhysics , 16(8):083028, aug 2014. ´Angel Morales Garc´ıa, Elena del Corro, Martin Kalbac,and Otakar Frank. Tuning the electronic properties ofmonolayer and bilayer transition metal dichalcogenidecompounds under direct out-of-plane compression. Phys.Chem. Chem. Phys. , 19:13333–13340, 2017. Steve Plimpton. Fast parallel algorithms for short-rangemolecular dynamics.
Journal of Computational Physics ,117(1):1 – 19, 1995.
Mit H. Naik, Indrajit Maity, Prabal K. Maiti, andManish Jain. Kolmogorov–crespi potential for multi-layer transition-metal dichalcogenides: Capturing struc-tural transformations in moir´e superlattices.
The Journalof Physical Chemistry C , 123(15):9770–9778, 2019. Jin-Wu Jiang and Yu-Ping Zhou. Handbook of stillinger-weber potential parameters for two-dimensional atomiccrystals. 2017. Nicola Marzari, Arash A. Mostofi, Jonathan R. Yates, IvoSouza, and David Vanderbilt. Maximally localized wan-nier functions: Theory and applications.
Rev. Mod. Phys. ,84:1419–1475, Oct 2012. Giovanni Pizzi, Valerio Vitale, Ryotaro Arita, StefanBl¨ugel, Frank Freimuth, Guillaume G´eranton, MarcoGibertini, Dominik Gresch, Charles Johnson, TakashiKoretsune, Julen Iba˜nez-Azpiroz, Hyungjun Lee, Jae-Mo Lihm, Daniel Marchand, Antimo Marrazzo, YuriyMokrousov, Jamal I Mustafa, Yoshiro Nohara, Yusuke Nomura, Lorenzo Paulatto, Samuel Ponc´e, Thomas Pon-weiser, Junfeng Qiao, Florian Th¨ole, Stepan S Tsirkin,Ma lgorzata Wierzbowska, Nicola Marzari, David Van-derbilt, Ivo Souza, Arash A Mostofi, and Jonathan RYates. Wannier90 as a community code: new featuresand applications.
Journal of Physics: Condensed Matter ,32(16):165902, jan 2020. S. S. Batsanov. Van der waals radii of elements.
InorganicMaterials , 37(9):871–885, Sep 2001. Lede Xian, Martin Claassen, Dominik Kiese, Michael M.Scherer, Simon Trebst, Dante M. Kennes, and Angel Ru-bio. Realization of nearly dispersionless bands with strongorbital anisotropy from destructive interference in twistedbilayer mos2, 2020. Kai Liu, Wei Luo, Junyi Ji, Paolo Barone, Silvia Picozzi,and Hongjun Xiang. Band splitting with vanishing spin po-larizations in noncentrosymmetric crystals.
Nature Com-munications , 10(1):5144, Nov 2019. Paolo Giannozzi, Stefano Baroni, Nicola Bonini, MatteoCalandra, Roberto Car, Carlo Cavazzoni, Davide Ceresoli,Guido L Chiarotti, Matteo Cococcioni, Ismaila Dabo, andet al. Quantum espresso: a modular and open-source soft-ware project for quantum simulations of materials.
Journalof Physics: Condensed Matter , 21(39):395502, 2009. Joseph C. A. Prentice, Jolyon Aarons, James C. Womack,Alice E. A. Allen, Lampros Andrinopoulos, Lucian An-ton, Robert A. Bell, Arihant Bhandari, Gabriel A. Bram-ley, Robert J. Charlton, Rebecca J. Clements, Daniel J.Cole, Gabriel Constantinescu, Fabiano Corsetti, Simon M.-M. Dubois, Kevin K. B. Duff, Jos´e Mar´ıa Escart´ın, An-drea Greco, Quintin Hill, Louis P. Lee, Edward Linscott,David D. O’Regan, Maximillian J. S. Phipps, Laura E.Ratcliff, ´Alvaro Ruiz Serrano, Edward W. Tait, GilbertoTeobaldi, Valerio Vitale, Nelson Yeung, Tim J. Zuehls-dorff, Jacek Dziedzic, Peter D. Haynes, Nicholas D. M.Hine, Arash A. Mostofi, Mike C. Payne, and Chris-Kriton Skylaris. The onetep linear-scaling density func-tional theory program.
The Journal of Chemical Physics ,152(17):174111, 2020. L. E. Ratcliff, G. J. Conduit, N. D. M. Hine, and P. D.Haynes. Band structure interpolation using optimized localorbitals from linear-scaling density functional theory.
Phys.Rev. B , 98:125123, Sep 2018. John P. Perdew, Kieron Burke, and Matthias Ernzerhof.Generalized gradient approximation made simple.
Phys.Rev. Lett. , 77:3865–3868, Oct 1996. P. E. Bl¨ochl. Projector augmented-wave method.
Phys.Rev. B , 50:17953–17979, Dec 1994. Fran¸cois Jollet, Marc Torrent, and Natalie Holzwarth.Generation of projector augmented-wave atomic data: A71 element validated table in the xml format.
Comput.Phys. Commun. , 185(4):1246 – 1254, 2014. Kevin F. Garrity, Joseph W. Bennett, Karin M. Rabe, andDavid Vanderbilt. Pseudopotentials for high-throughputdft calculations.
Computational Materials Science , 81:446–452, 2014. r pp ( Å ) V pp ( e V ) (a) z pp =2.62 Åz pp =2.72 Åz pp =2.82 Åz pp =2.92 Åz pp =3.02 Å z pp =3.12 Åz pp =3.22 Åz pp =3.32 Åz pp =3.42 Å r pd ( Å ) V pd ( e V ) (b) z pd =4.21 Åz pd =4.31 Åz pd =4.41 Åz pd =4.51 Åz pd =4.61 Å z pd =4.71 Åz pd =4.81 Åz pd =4.91 Åz pd =5.01 Å FIG. 15: Distance dependence of the Slater-Koster parameters: (a) p-p interlayer hopping with dashed linesdenoting V ppσ and solid lines denoting V ppπ as function of interlayer separation and (b) p z -d z interlayer hoppingparameter V pdπ as function of interlayer separation. upplementary information for “Chemicaltrends in flat band properties of twistedtransition metal dichalcogenide homo- andheterobilayers” Valerio Vitale , Kemal Atalar , Arash A. Mostofi and Johannes Lischner Departments of Materials and Physics and the Thomas Young Centre forTheory and Simulation of Materials, Imperial College London, LondonSW7 2AZ, UK1 a r X i v : . [ c ond - m a t . m t r l - s c i ] F e b Fig. S1 shows the band structures of unrelaxed (flat) and relaxed homobilay-ers at θ = 5 . ° . The effect of relaxation on flat bands is to close the gap atthe K-point. Moreover, in WS /WS and WSe /WSe bilayers the orderingof dispersive vs flat bands is reverseSimilarly, unrelaxed (flat) and relaxed band structures for heterobilayerswith same chalcogen species at θ = 5 . ° are shown in Fig. S2. The orderingof dispersive vs flat bands is reversed by relaxation in WS /MoS but not inWSe /MoSe . The gap between the flat bands at the K-point is closed.Fig. S3 shows the band structures for unrelaxed and relaxed heterobilayerswith different chalcogen species at θ = 4 . ° . In these systems relaxation doesnot change the ordering of flat bands vs dispersive bands. In MoSe /WS and WSe /MoS flat bands are always higher in energy, the effect of relax-ation is to make them separated from all remote bands. In WSe /MoS andWSe /WS bilayers relaxation has a minor effect, as dispersive bands aremuch higher in energy than flat bands.2 G M K G E ne r g y ( e V ) E ne r g y ( e V ) -0.2-0.1 0 G M K G E ne r g y ( e V ) E ne r g y ( e V ) -0.2-0.1 0 G M K G E ne r g y ( e V ) E ne r g y ( e V ) -0.2-0.1 0 G M K G E ne r g y ( e V ) E ne r g y ( e V ) -0.2-0.1 0 G M K G E ne r g y ( e V ) E ne r g y ( e V )
MoS /MoS MoSe /MoSe WS /WS WeS /WSe -0.2-0.1 0 G M K G E ne r g y ( e V ) E ne r g y ( e V ) -0.2-0.1 0 G M K G E ne r g y ( e V ) E ne r g y ( e V ) -0.2-0.1 0 G M K G E ne r g y ( e V ) E ne r g y ( e V ) Figure S1: Comparison of flat(left column) vs relaxed (right column) band-structures of (a)-(b) MoS /MoS , (c)-(d) MoSe /MoSe , (e)-(f)WS /WS and (g)-(h) WSe /WSe at θ = 5 . ◦ . For the flat systems only the interlayerseparation has been relaxed. 3 G M K G E ne r g y ( e V ) E ne r g y ( e V ) -0.2-0.1 0 G M K G E ne r g y ( e V ) E ne r g y ( e V ) -0.2-0.1 0 G M K G E ne r g y ( e V ) E ne r g y ( e V )
MoSe /MoS WSe /WS MoSe /WS WeS /MoS -0.1 0 G M K G E ne r g y ( e V ) E ne r g y ( e V ) -0.1 0 G M K G E ne r g y ( e V ) E ne r g y ( e V ) -0.2-0.1 0 G M K G E ne r g y ( e V ) E ne r g y ( e V ) Figure S3: Comparison of flat (left column) vs relaxed (right column) band-structures of (a)-(b) MoSe /MoS , (c)-(d) WSe /MoS , (e)-(f) MoSe /WS and (g)-(h) WSe /WS at θ = 4 . ◦ . For the flat systems only the interlayerseparation has been relaxed. 5 Fig. S4 shows the comparison between DFT and tight-binding (TB) bandstructures for relaxed homobilayers at two angles, θ = 7 . ° and θ = 5 . ° ,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 cut-off of 800 eV. A basis consisting of 9 non-orthogonal generalized Wannierfunctions for both metal and chalcogen atoms is employed. Both DFT andTB band structures do not include spin-orbit coupling. Agreement betweenDFT and TB band structures improves at smaller angles.6 Γ M K Γ E ne r g y ( m e V ) TB w/o SOCDFT w/o SOC -300-250-200-150-100-50 0 50 100 Γ M K Γ E ne r g y ( m e V ) TB w/o SOCDFT w/o SOC
MoS /MoS WSe /WSe MoSe /MoSe WS /WS q =7.3° q = -200-150-100-50 0 50 100 Γ M K Γ E ne r g y ( m e V ) TB w/o SOCDFT w/o SOC -200-150-100-50 0 50 100 Γ M K Γ E ne r g y ( m e V ) TB w/o SOCDFT w/o SOC -300-250-200-150-100-50 0 50 100 Γ M K Γ E ne r g y ( m e V ) TB w/o SOCDFT w/o SOC -300-250-200-150-100-50 0 50 100 150 200 Γ M K Γ E ne r g y ( m e V ) TB w/o SOCDFT w/o SOC -300-250-200-150-100-50 0 50 100 150 200 Γ M K Γ E ne r g y ( m e V ) TB w/o SOCDFT w/o SOC -200-150-100-50 0 50 100 Γ M K Γ E ne r g y ( m e V ) TB w/o SOCDFT w/o SOC
Figure S4: Comparison of relaxed DFT (solid black) vs TB (blue dots)bandstructures at two different twist angles, θ = 7 . ° (left) and θ = 5 . ° (right) for a)-b) MoS /MoS , c)-d) MoSe /MoSe , e)-f) WS /WS and g)-h)WSe /WSe at θ = 4 . ◦ . 7 Fig. S5 shows the projections of the highest valence bands states onto Wd z orbitals and inner Se p z orbitals in WSe /WSe at θ = 3 . ° . Flat bandshave large projections on these orbitals ( ∼ monolayers, which also havelarge projections onto W d z orbitals and Se p z orbitals. For homobilayersprojections onto the two layers are symmetric. Γ M K Γ E ne r g y ( m e V ) Se (top layer) : p z Γ M K Γ E ne r g y ( m e V ) W (top layer) : d z Γ M K Γ E ne r g y ( m e V ) Se (bottom layer) : p z Γ M K Γ E ne r g y ( m e V ) W (bottom layer) : d z
Fig. S10 shows the effect of interlayer separation (ILS) on the ordering of flatbands over dispersive bands in WSe /WS at θ = 4 . ° . To compute the bandstructures in Fig. S10, we start from a relaxed WSe /WS bilayer and rigidlytranslate the layers in the z -direction without further relaxing them. ForILS greater or equal to the equilibrium ILS we find dispersive bands at thevalence band edge (top panels in Fig. S10). As the ILS is reduced (bottompanels in Fig. S10), flat bands emerge on top of the valence band edge, whichsuggests that pressure can be used to modify flat bands properties in thesesystems. 12 G M K G E ne r g y ( e V ) ILS eq + 1.0 Å −0.8−0.6−0.4−0.2 0 0.2 G M K G E ne r g y ( e V ) ILS eq + 0.5 Å −0.8−0.6−0.4−0.2 0 0.2 G M K G E ne r g y ( e V ) ILS eq − 0.5 Å −0.8−0.6−0.4−0.2 0 0.2 G M K G E ne r g y ( e V ) ILS eq − 1.0 Å Figure S10: Effect of interlayer separation on the ordering of flat vs dispersivebands at the top valence manifold in WSe /WS at θ = 4 . ° . In all panels theILS is increased/decreased by rigidly transposing the layers from the relaxedstructure and no further relaxation is considered. Top left panel the ILS isincreased by 1˚A. In top right panel the ILS is increased by 0.5˚A. In bottomleft panel the ILS is decreased by 0.5˚A. In bottom right panel the ILS isdecreased by 1˚A. In all panels the bandstructure at equilibrium is shownwith empty black circles and that with modified ILS with full red circles.The ordering of flat vs dispersive bands is reversed when the two layers arecloser, compared to the equilibrium case (which exhibit dispersive bands ontop), as shown in bottom panels (which exhibit flat bands on top).13 In this section we present tables of tigh-binding hopping parameters, bothintralayer and interlayer, for all homo- and hetero-bilayers in our work.14able S1: Tight-binding independent parameters in units of eV for MoS ,MoSe , WS and WSe from Wannierised DFT results. MoS MoSe WS WSe ε = ε .
034 1 .
037 0 .
679 1 . ε − . − . − . − . ε = ε − . − . − . − . ε − .
125 0 . − .
695 0 . ε = ε − .
916 0 . − .
470 0 . ε − . − . − . − . ε = ε − . − . − . − . t (1)1 , − . − . − . − . t (1)2 , .
031 0 .
019 0 .
027 0 . t (1)3 , − . − . − . − . t (1)4 , .
867 0 .
953 0 .
877 0 . t (1)5 , − . − . − . − . t (1)6 , − . − . − . − . t (1)7 , .
284 0 .
249 0 .
359 0 . t (1)8 , − . − . − . − . t (1)9 , − . − . − . − . t (1)10 , .
916 0 .
989 0 .
972 1 . t (1)11 , . − .
003 0 .
010 0 . t (1)3 , − . − . − . − . t (1)6 , .
409 0 .
358 0 .
488 0 . t (1)9 , . − . − . − . t (1)1 , − . − . − . − . t (1)3 , − . − . − . − . t (1)4 , − . − . − . − . t (1)6 , − . − . − . − . t (1)7 , − . − . − . − . t (1)9 , .
105 0 .
118 0 .
161 0 . t (1)10 , − . − . − . − . t (5)4 , − . − . − . − . t (5)3 , − . − . − . − . t (5)5 , .
102 1 .
971 2 .
268 2 . t (5)9 , − . − . − . − . t (5)11 , − . − . − . − . t (5)10 , .
364 1 .
280 1 .
522 1 . t (5)9 , − . − . − . − . t (5)11 , .
618 0 .
570 0 .
646 0 . t (6)9 , − . − . − . − . t (6)11 , − . − . − . − . t (6)9 , − . − . − . − . t (6)11 , − . − . − . − . z -p z interlayer hopping parameters for MoS /MoS ∆ z pp [˚A] ν [eV] R [˚A] ησ π σ π σ π .
617 4 . − .
271 2 .
824 2 .
612 3 .
547 5 . .
717 3 . − .
112 2 .
888 2 .
688 3 .
590 5 . .
817 3 . − .
983 2 .
952 2 .
759 3 .
630 5 . .
917 3 . − .
885 3 .
017 2 .
822 3 .
666 5 . .
017 2 . − .
818 3 .
080 2 .
872 3 .
696 5 . .
117 2 . − .
792 3 .
143 2 .
899 3 .
722 5 . .
217 2 . − .
843 3 .
205 2 .
879 3 .
741 4 . .
317 2 . − .
138 3 .
265 2 .
739 3 .
754 4 . .
417 1 . − .
142 3 .
323 1 .
959 3 .
760 2 . Table S3: p z -d z interlayer hopping parameters for MoS /MoS ∆ z pd [˚A] V [eV] α β γσ π σ π σ π σ π . − . − . − . − .
332 2 .
386 1 . − . − . . − . − . − . − .
186 2 .
245 1 . − . − . . − . − . − . − .
050 2 .
171 1 . − . − . . − . − . − . − .
925 2 .
120 1 . − . − . . − . − . − . − .
813 2 .
081 1 . − . − . . − . − . − . − .
713 2 .
049 1 . − . − . . − . − . − . − .
627 2 .
022 1 . − . − . . − . − . − . − .
549 1 .
998 1 . − . − . . − . − . − . − .
482 1 .
976 1 . − . − . z -p z interlayer hopping parameters for MoSe /MoSe ∆ z pp [˚A] ν [eV] R [˚A] ησ π σ π σ π .
485 4 . − .
970 2 .
881 2 .
513 3 .
642 5 . .
585 4 . − .
745 2 .
940 2 .
588 3 .
684 5 . .
685 3 . − .
554 3 .
001 2 .
660 3 .
726 5 . .
785 3 . − .
398 3 .
063 2 .
728 3 .
768 5 . .
885 3 . − .
277 3 .
127 2 .
787 3 .
808 5 . .
985 3 . − .
197 3 .
191 2 .
834 3 .
845 5 . .
085 2 . − .
176 3 .
255 2 .
857 3 .
879 4 . .
185 2 . − .
276 3 .
319 2 .
830 3 .
908 4 . .
285 2 . − .
846 3 .
381 2 .
654 3 .
932 3 . Table S5: p z -d z interlayer hopping parameters for MoSe /MoSe ∆ z pd [˚A] V [eV] α β γσ π σ π σ π σ π . − . − . − . − .
558 2 .
770 1 . − . − . . − . − . − . − .
379 2 .
340 1 . − . − . . − . − . − . − .
195 2 .
141 1 . − . − . . − . − . − . − .
015 2 .
068 1 . − . − . . − . − . − . − .
843 2 .
022 1 . − . − . . − . − . − . − .
682 1 .
987 1 . − . − . . − . − . − . − .
535 1 .
960 1 . − . − . . − . − . − . − .
402 1 .
938 1 . − . − . . − . − . − . − .
283 1 .
918 1 . − . − . z -p z interlayer hopping parameters for WS /WS ∆ z pp [˚A] ν [eV] R [˚A] ησ π σ π σ π .
618 3 . − .
021 2 .
847 2 .
704 3 .
546 6 . .
718 3 . − .
859 2 .
915 2 .
799 3 .
598 6 . .
818 3 . − .
831 2 .
980 2 .
837 3 .
636 6 . .
918 2 . − .
714 3 .
040 2 .
922 3 .
658 6 . .
018 2 . − .
710 3 .
105 2 .
945 3 .
691 5 . .
118 2 . − .
705 3 .
173 2 .
962 3 .
724 5 . .
218 2 . − .
810 3 .
238 2 .
906 3 .
749 4 . .
318 2 . − .
538 3 .
302 2 .
592 3 .
769 3 . Table S7: p z -d z interlayer hopping parameters for WS /WS ∆ z pd [˚A] V [eV] α β γσ π σ π σ π σ π . − . − . − . − .
761 2 .
153 1 . − . − . . − . − . − . − .
369 2 .
115 1 . − . − . . − . − . − . − .
324 2 .
052 1 . − . − . . − . − . − . − .
241 2 .
045 1 . − . − . . − . − . − . − .
990 2 .
029 1 . − . − . . − . − . − . − .
836 2 .
019 1 . − . − . . − . − . − . − .
701 2 .
010 1 . − . − . . − . − . − . − .
595 2 .
003 1 . − . − . z -p z interlayer hopping parameters for WSe /WSe ∆ z pp [˚A] ν [eV] R [˚A] ησ π σ π σ π .
486 3 . − .
042 3 .
020 2 .
757 3 .
949 7 . .
586 4 . − .
012 2 .
925 2 .
799 3 .
557 7 . .
686 3 . − .
834 3 .
010 2 .
910 3 .
668 7 . .
786 3 . − .
779 3 .
024 2 .
970 3 .
590 7 . .
886 3 . − .
019 3 .
096 2 .
891 3 .
646 5 . .
986 2 . − .
838 3 .
225 2 .
998 3 .
840 6 . .
086 2 . − .
737 3 .
294 3 .
076 3 .
881 6 . .
186 2 . − .
709 3 .
369 3 .
110 3 .
938 5 . .
286 2 . − . − .
434 3 .
114 3 .
962 5 . Table S9: p z -d z interlayer hopping parameters for WSe /WSe ∆ z pd [˚A] V [eV] α β γσ π σ π σ π σ π . − . − . − . − .
743 1 .
852 1 . − . − . . − . − . − . − .
161 2 .
654 1 . − . − . . − . − . − . − .
867 1 .
906 1 . − . − . . − . − . − . − .
275 1 .
984 1 . − . − . . − . − . − . − .
101 1 .
983 1 . − . − . . − . − . − . − .
286 1 .
981 1 . − . − . . − . − . − . − .
973 2 .
003 1 . − . − . . − . − . − . − .
782 1 .
996 1 . − . − . . − . − . − . − .
913 1 .
991 1 . − . − . z -p z interlayer hopping parameters for MoSe /MoS ∆ z pp [˚A] ν [eV] R [˚A] ησ π σ π σ π .
617 3 . − .
393 2 .
931 2 .
641 3 .
838 5 . .
717 3 . − .
224 2 .
996 2 .
718 3 .
886 5 . .
817 3 . − .
080 3 .
062 2 .
794 3 .
933 5 . .
917 3 . − .
968 3 .
130 2 .
862 3 .
979 5 . .
017 2 . − .
883 3 .
198 2 .
922 4 .
020 5 . .
117 2 . − .
829 3 .
265 2 .
969 4 .
056 5 . .
217 2 . − .
826 3 .
331 2 .
988 4 .
086 5 . .
317 2 . − .
929 3 .
397 2 .
946 4 .
112 4 . .
417 1 . − .
576 3 .
462 2 .
702 4 .
132 3 . Table S11: p z -d z and d z - p z interlayer hopping parameters for MoSe /MoS ∆ z pd [˚A] V [eV] α β γσ π σ π σ π σ π . − . − . − . − .
118 2 .
397 1 . − . − . . − . − . − . − .
916 2 .
257 1 . − . − . . − . − . − . − .
725 2 .
190 1 . − . − . . − . − . − . − .
578 2 .
144 1 . − . − . . − . − . − . − .
458 2 .
110 1 . − . − . . − . − . − . − .
362 2 .
081 1 . − . − . . − . − . − . − .
282 2 .
058 1 . − . − . . − . − . − . − .
219 2 .
037 1 . − . − . . − . − . − . − .
166 2 .
018 1 . − . − . z dp [˚A]4 .
208 0 .
647 2 . − . − .
926 2 .
269 1 . − . − . .
309 0 .
668 3 . − . − .
779 2 .
180 1 . − . − . .
409 0 .
716 3 . − . − .
647 2 .
109 1 . − . − . .
508 0 .
792 3 . − . − .
581 2 .
054 1 . − . − . .
609 0 .
884 3 . − . − .
547 2 .
005 1 . − . − . .
708 0 .
996 3 . − . − .
553 1 .
962 1 . − . − . .
809 1 .
131 4 . − . − .
599 1 .
924 1 . − . − . .
909 1 .
289 4 . − . − .
674 1 .
891 1 . − . − . .
008 1 .
475 5 . − . − .
774 1 .
860 1 . − . − . z -p z interlayer hopping parameters for WS /MoS ∆ z pp [˚A] ν [eV] R [˚A] ησ π σ π σ π .
617 4 . − .
182 2 .
824 2 .
644 3 .
522 5 . .
717 3 . − .
038 2 .
889 2 .
719 3 .
566 5 . .
817 3 . − .
926 2 .
953 2 .
787 3 .
604 5 . .
917 3 . − .
845 3 .
017 2 .
845 3 .
638 5 . .
017 2 . − .
794 3 .
082 2 .
888 3 .
670 5 . .
117 2 . − .
797 3 .
146 2 .
898 3 .
696 5 . .
217 2 . − .
915 3 .
209 2 .
840 3 .
717 4 . .
317 2 . − .
581 3 .
270 2 .
576 3 .
731 3 . Table S13: p z -d z and d z - p z interlayer hopping parameters for WS /MoS ∆ z pd [˚A] V [eV] α β γσ π σ π σ π σ π . − . − . − . − .
109 2 .
204 1 . − . − . . − . − . − . − .
944 2 .
148 1 . − . − . . − . − . − . − .
778 2 .
111 1 . − . − . . − . − . − . − .
645 2 .
082 1 . − . − . . − . − . − . − .
534 2 .
058 1 . − . − . . − . − . − . − .
454 2 .
039 1 . − . − . . − . − . − . − .
383 2 .
021 1 . − . − . . − . − . − . − .
291 2 .
003 1 . − . − . z dp [˚A]4 .
208 0 .
488 2 . − . − .
057 2 .
416 1 . − . − . .
309 0 .
450 2 . − . − .
929 2 .
290 1 . − . − . .
409 0 .
468 2 . − . − .
691 2 .
209 1 . − . − . .
508 0 .
510 3 . − . − .
598 2 .
147 1 . − . − . .
609 0 .
564 3 . − . − .
492 2 .
095 1 . − . − . .
708 0 .
634 3 . − . − .
477 2 .
049 1 . − . − . .
809 0 .
718 3 . − . − .
497 2 .
009 1 . − . − . .
909 0 .
812 4 . − . − .
561 1 .
978 1 . − . − . z -p z interlayer hopping parameters for WSe /MoS ∆ z pp [˚A] ν [eV] R [˚A] ησ π σ π σ π .
618 3 . − .
357 2 .
936 2 .
646 3 .
793 5 . .
718 3 . − .
191 3 .
005 2 .
724 3 .
850 5 . .
818 3 . − .
062 3 .
080 2 .
795 3 .
920 5 . .
918 2 . − .
979 3 .
144 2 .
852 3 .
953 5 . .
018 2 . − .
937 3 .
212 2 .
891 4 .
000 5 . .
118 2 . − .
922 3 .
278 2 .
915 4 .
032 5 . .
218 2 . − .
020 3 .
348 2 .
881 4 .
073 4 . .
318 2 . − .
598 3 .
419 2 .
672 4 .
116 3 . Table S15: p z -d z and d z - p z interlayer hopping parameters for WSe /MoS ∆ z pd [˚A] V [eV] α β γσ π σ π σ π σ π . − . − . − . − .
414 2 .
282 1 . − . − . . − . − . − . − .
375 2 .
222 1 . − . − . . − . − . − . − .
321 2 .
184 1 . − . − . . − . − . − . − .
165 2 .
152 1 . − . − . . − . − . − . − .
087 2 .
123 1 . − . − . . − . − . − . − .
053 2 .
101 1 . − . − . . − . − . − . − .
997 2 .
082 1 . − . − . . − . − . − . − .
978 2 .
067 1 . − . − . z dp [˚A]4 .
209 0 .
326 2 . − . − .
266 2 .
044 1 . − . − . .
309 0 .
371 1 . − . − .
472 2 .
015 1 . − . − . .
409 0 .
456 2 . − . − .
497 1 .
999 1 . − . − . .
509 0 .
535 2 . − . − .
022 1 .
973 1 . − . − . .
609 0 .
620 3 . − . − .
946 1 .
930 1 . − . − . .
709 0 .
701 3 . − . − .
953 1 .
910 1 . − . − . .
809 0 .
817 3 . − . − .
104 1 .
892 1 . − . − . .
909 0 .
953 3 . − . − .
815 1 .
880 1 . − . − . z -p z interlayer hopping parameters for WS /MoS ∆ z pp [˚A] ν [eV] R [˚A] ησ π σ π σ π .
618 3 . − .
039 2 .
938 2 .
753 3 .
830 6 . .
718 3 . − .
961 3 .
024 2 .
812 3 .
939 6 . .
818 3 . − .
822 3 .
090 2 .
902 3 .
986 6 . .
918 2 . − .
726 3 .
160 2 .
979 4 .
037 6 . .
018 2 . − .
653 3 .
230 3 .
048 4 .
082 6 . .
118 2 . − .
616 3 .
299 3 .
095 4 .
122 6 . .
218 2 . − .
603 3 .
370 3 .
124 4 .
164 6 . .
318 2 . − .
640 3 .
441 3 .
110 4 .
201 5 . .
418 1 . − .
951 3 .
511 2 .
924 4 .
235 4 . Table S17: p z -d z and d z - p z interlayer hopping parameters for WS /MoS ∆ z pd [˚A] V [eV] α β γσ π σ π σ π σ π . − . − . − . − .
546 2 .
375 1 . − . − . . − . − . − . − .
305 2 .
224 1 . − . − . . − . − . − . − .
187 2 .
145 1 . − . − . . − . − . − . − .
066 2 .
117 1 . − . − . . − . − . − . − .
986 2 .
096 1 . − . − . . − . − . − . − .
941 2 .
077 1 . − . − . . − . − . − . − .
901 2 .
059 1 . − . − . . − . − . − . − .
834 2 .
044 1 . − . − . . − . − . − . − .
760 2 .
030 1 . − . − . z dp [˚A]4 .
209 0 .
475 2 . − . − .
165 2 .
158 1 . − . − . .
309 0 .
513 2 . − . − .
958 2 .
089 1 . − . − . .
409 0 .
586 2 . − . − .
860 2 .
068 1 . − . − . .
509 0 .
652 2 . − . − .
817 2 .
025 1 . − . − . .
609 0 .
735 2 . − . − .
871 1 .
989 1 . − . − . .
709 0 .
830 3 . − . − .
945 1 .
959 1 . − . − . .
809 0 .
941 3 . − . − .
024 1 .
932 1 . − . − . .
909 1 .
073 3 . − . − .
137 1 .
909 1 . − . − . .
009 1 .
207 5 . − . − .
592 1 .
857 1 . − . − . z -p z interlayer hopping parameters for WSe /WS ∆ z pp [˚A] ν [eV] R [˚A] ησ π σ π σ π .
618 3 . − .
576 2 .
933 2 .
580 3 .
745 5 . .
718 3 . − .
230 3 .
024 2 .
706 3 .
867 5 . .
818 3 . − .
331 3 .
101 2 .
699 3 .
943 5 . .
918 2 . − .
394 3 .
189 2 .
702 4 .
046 4 . .
018 2 . − .
432 3 .
273 2 .
708 4 .
136 4 . .
118 2 . − .
527 3 .
355 2 .
696 4 .
216 4 . .
218 2 . − .
672 3 .
444 2 .
666 4 .
310 4 . .
318 1 . − .
497 3 .
528 2 .
482 4 .
391 3 . Table S19: p z -d z and d z - p z interlayer hopping parameters for WSe /WS ∆ z pd [˚A] V [eV] α β γσ π σ π σ π σ π . − . − . − . − .
433 2 .
511 1 . − . − . . − . − . − . − .
151 2 .
352 1 . − . − . . − . − . − . − .
739 2 .
300 1 . − . − . . − . − . − . − .
534 2 .
263 1 . − . − . . − . − . − . − .
507 2 .
238 1 . − . − . . − . − . − . − .
272 2 .
215 1 . − . − . . − . − . − . − .
253 2 .
198 1 . − . − . . − . − . − . − .
115 2 .
182 1 . − . − . z dp [˚A]4 .
209 0 .
412 1 . − . − .
459 2 .
119 1 . − . − . .
309 0 .
386 1 . − . − .
068 1 .
985 1 . − .
979 0 . .
409 0 .
419 1 . − . − .
216 1 .
949 1 . − .
941 0 . .
509 0 .
488 1 . − . − .
333 1 .
954 1 . − .
124 0 . .
609 0 .
534 1 . − . − .
443 1 .
933 1 . − .
156 0 . .
709 0 .
583 1 . − . − .
744 1 .
916 1 . − . − . .
809 0 .
677 2 . − . − .
769 1 .
907 1 . − .
326 0 . .
909 0 .
774 2 . − . − .
957 1 .
888 1 . − . − . z -p z interlayer hopping parameters for WSe /MoSe ∆ z pp [˚A] ν [eV] R [˚A] ησ π σ π σ π .
485 4 . − .
791 2 .
873 2 .
552 3 .
574 5 . .
585 4 . − .
467 2 .
930 2 .
663 3 .
614 5 . .
685 3 . − .
264 2 .
998 2 .
749 3 .
676 5 . .
785 3 . − .
111 3 .
066 2 .
826 3 .
732 5 . .
885 3 . − .
978 3 .
130 2 .
903 3 .
775 5 . .
985 3 . − .
880 3 .
193 2 .
969 3 .
806 5 . .
085 2 . − .
824 3 .
260 3 .
014 3 .
847 5 . .
185 2 . − .
841 3 .
326 3 .
015 3 .
885 5 . .
285 2 . − .
509 3 .
388 2 .
719 3 .
905 3 . Table S21: p z -d z and d z - p z interlayer hopping parameters for WSe /WS ∆ z pd [˚A] V [eV] α β γσ π σ π σ π σ π . − . − . − . − .
639 2 .
731 1 . − . − . . − . − . − . − .
151 2 .
061 1 . − . − . . − . − . − . − .
970 2 .
011 1 . − . − . . − . − . − . − .
733 1 .
980 1 . − . − . . − . − . − . − .
544 1 .
963 1 . − . − . . − . − . − . − .
343 1 .
949 1 . − . − . . − . − . − . − .
258 1 .
938 1 . − . − . . − . − . − . − .
120 1 .
931 1 . − . − . . − . − . − . − .
003 1 .
912 1 . − . − . z dp [˚A]4 .
143 0 .
943 3 . − . − .
183 2 .
460 1 . − . − . .
242 1 .
084 2 . − . − .
422 2 .
502 1 . − . − . .
343 0 .
585 3 . − . − .
534 2 .
282 1 . − . − . .
442 0 .
677 3 . − . − .
464 2 .
252 1 . − . − . .
543 0 .
712 3 . − . − .
233 2 .
197 1 . − . − . .
643 0 .
827 3 . − . − .
776 2 .
134 1 . − . − . .
742 0 .
897 2 . − . − .
488 2 .
077 1 . − . − . .
843 0 .
912 3 . − . − .
478 2 .
012 1 . − . − . .
942 1 .
047 3 . − . − .
676 1 .
962 1 . − . − . eferences [1] Joseph C. A. Prentice, Jolyon Aarons, James C. Womack, Alice E. A.Allen, Lampros Andrinopoulos, Lucian Anton, Robert A. Bell, ArihantBhandari, Gabriel A. Bramley, Robert J. Charlton, Rebecca J. Clements,Daniel J. Cole, Gabriel Constantinescu, Fabiano Corsetti, Simon M.-M.Dubois, Kevin K. B. Duff, Jos´e Mar´ıa Escart´ın, Andrea Greco, QuintinHill, Louis P. Lee, Edward Linscott, David D. O’Regan, MaximillianJ. S. Phipps, Laura E. Ratcliff, ´Alvaro Ruiz Serrano, Edward W. Tait,Gilberto Teobaldi, Valerio Vitale, Nelson Yeung, Tim J. Zuehlsdorff,Jacek Dziedzic, Peter D. Haynes, Nicholas D. M. Hine, Arash A. Mostofi,Mike C. Payne, and Chris-Kriton Skylaris. The onetep linear-scalingdensity functional theory program. The Journal of Chemical Physics ,152(17):174111, 2020.[2] L. E. Ratcliff, G. J. Conduit, N. D. M. Hine, and P. D. Haynes. Bandstructure interpolation using optimized local orbitals from linear-scalingdensity functional theory.
Phys. Rev. B , 98:125123, Sep 2018.[3] John P. Perdew, Kieron Burke, and Matthias Ernzerhof. Generalizedgradient approximation made simple.
Phys. Rev. Lett. , 77:3865–3868,Oct 1996.[4] P. E. Bl¨ochl. Projector augmented-wave method.
Phys. Rev. B , 50:17953–17979, Dec 1994.[5] Fran¸cois Jollet, Marc Torrent, and Natalie Holzwarth. Generation ofprojector augmented-wave atomic data: A 71 element validated table inthe xml format.
Comput. Phys. Commun. , 185(4):1246 – 1254, 2014.[6] Kevin F. Garrity, Joseph W. Bennett, Karin M. Rabe, and David Van-derbilt. Pseudopotentials for high-throughput dft calculations.