Quantum Spin Hall Effect in Two-dimensional Crystals of Transition Metal Dichalcogenides
QQuantum Spin Hall Effect in Two-dimensional Crystalsof Transition Metal Dichalcogenides
M. A. Cazalilla,
1, 2
H. Ochoa, and F. Guinea Department of Physics, National Tsing Hua University,and National Center for Theoretical Sciences (NCTS), Hsinchu City, Taiwan. Donostia International Physics Center (DIPC),Manuel de Lardizabal 4, 20018, San Sebastian, Spain. Instituto de Ciencia de Materiales de Madrid (ICMM). CSIC, Cantoblanco. E-28049 Madrid. Spain. (Dated: September 4, 2018)We propose to engineer time-reversal-invariant topological insulators in two-dimensional (2D)crystals of transition metal dichalcogenides (TMDCs). We note that, at low doping, semiconductingTMDCs under shear strain will develop spin-polarized Landau levels residing in different valleys. Weargue that gaps between Landau levels in the range of 10 −
100 Kelvin are within experimental reach.In addition, we point out that a superlattice arising from a Moir´e pattern can lead to topologicallynon-trivial subbands. As a result, the edge transport becomes quantized, which can be probedin multi-terminal devices made using strained 2D crystals and/or heterostructures. The strong d character of valence and conduction bands may also allow for the investigation of the effects ofelectron correlations on the topological phases. PACS numbers: 85.75.-d,73.43.Qt,73.43.-f
There is currently much interest in two-dimensional(2D) crystals of transition metal dichalcogenides(TMDC) such as MoS or WSe [1]. Compared toGraphene, many of these materials are semiconductorswith sizeable gaps ( ≈ non-dissipative nano-electronics [4]. Indeed, provided thatmagnetic impurities (or other time-reversal symmetry-breaking perturbations) can be eliminated, electrontransport through the edge of a QSH insulator becomesballistic. This is because counter propagating channelsof the same spin are spatially separated.Nevertheless, the QSH effect has been observed invery few materials [5, 6] so far, making the search fornew systems exhibiting this effect a major scientific en-deavor. In the context of 2D crystals, Weeks and cowork-ers have recently proposed an approach to enhance thespin-orbit coupling in graphene by heavy ad-atom depo-sition [7], which could lead to a realization of the Kane-Mele model [8] and the QSH effect. The availabilityof 2D materials exhibiting this effect may allow for theconstruction of flexible electronic devices with low-powerconsumption.In this work, we present a proposal for engineeringthe QSHE in strained 2D crystals and heterostructuresmade using TMDCs. Strain is known to induce pseudo-magnetic fields in 2D crystals [9–11]. These fields havebeen predicted [9] and experimentally found [10] to pro-duce Landau levels (LLs) in the electronic spectrum of2D crystals such as graphene. As explained below, itis also possible to exploit these pseudo-magnetic fields to engineer time-reversal invariant topological phases inother 2D crystals such as the TMDCs. This approachhas several attractive features. Using MoS as an exam-ple, we find that gaps between LLs scale as (cid:126) ω c /k B (cid:39) . B [ T ] K, where B [ T ] is the strength of the pseudo-magnetic field in Tesla. Fields B [ T ] ∼ − T havebeen experimentally demonstrated [9, 10], which in thecase of TMDCs could lead to LL gaps of up to ≈
100 K.By comparison, the gaps achievable by straining GaAsare in the range of tens of mili-Kelvin [12]. The muchlarger magnitude of the LLs gaps in the present proposalstems the strain coupling to the electron current, ratherthan the electron spin current as in Ref. [12]. These twocoupling constants are, in general, vastly different in or-der of magnitude. However, in multi-valley systems likegraphene, strain only leads to (spin) unpolarized LLs.What sets semiconducting TMDCs apart is the spin-orbit coupling (SOC) that produces a large spin split-ting of the valence (and to a smaller extent, the conduc-tion) band. For small doping, this leads to spin-polarizedLLs in different valleys, which opens the possibility of re-alizing time-reversal invariant (TRI) topological phases.Furthermore, the valence and conduction band of semi-conducting TMDC have strong d character [13], which,along with the poor screening of the Coulomb interactionin 2D, means that electron interactions can have inter-esting effects on properties of the TRI topological phasesrealized in TMDCs.Bulk TMDC are composed of X - M - X layers stackedon top of each other and coupled by weak Van der Waalsforces. Therefore, like graphite, these materials can beexfoliated down to a single layer. The transition metalatoms ( M ) are arranged in a triangular lattice and eachone is bonded to six chalcogen atoms ( X ) (see Fig. 1 for a r X i v : . [ c ond - m a t . m e s - h a ll ] D ec People
Good approximation G K M G - E H e V L XX M j G KK ' M --++++-- = t z + +--+- = t x + = t y +- -+- FIG. 1: (Color Online) Top view of the 2D lattice of TMDC.The unit cell and Brillouin zone are shown in the right panel. M corresponds to the transition metal atom, whereas X tothe chalcogen atoms. a top view of the lattice). The point group of the 2Dcrystal is D h . The Fermi level lies at the two inequiv-alent corners of the Brillouin zone, K and K (cid:48) points,as in the case of graphene, but the spectrum is gapped.The character of the conduction and valence bands isdominated by d z − r and d x − y ± id xy orbitals from M atom, respectively, where the sign + ( − ) holds for K ( K (cid:48) )point. Furthermore, both bands have a non-negligible hy-bridization with the bonding combination of p x ± ip y (inthe conduction band) and p x ∓ ip y (in the valence band)orbitals from X atoms.In the continuum limit, the system is well described bythe following two-band Hamiltonian [2, 13, 15]: H = v ( τ z σ x p x + σ y p y ) + ∆ σ z + λ SO τ z s z (1 − σ z ) , (1)where p = ( p x , p y ) is the electron crystal momentumreferred to the K ( K (cid:48) ) point, for which we introduce τ z = +1 ( τ z = − σ α are Pauli matrices acting on thespace span by conduction and valence band states; ∆ isthe band gap and v = ta/ (cid:126) , where a is the lattice parame-ter and t is a phenomenological hybridization parameter;finally, s z is the spin projection along the z axis perpen-dicular to the 2D crystal. The Hamiltonian (1) can beconstructed phenomenologically by considering scalar in-variants of D h formed from operators σ α , p and spin s α .The symmetry properties of these operators are listed inTable I. Besides the gap, another important differencewith graphene is the spin splitting of the bands due tothe strong spin-orbit interaction provided by M atoms.Below, we focus on p-doped TMDCs, which allows usto neglect the smaller spin splitting of the conductionband [15]. For the valence band, the latter is in the rangeof 100-400 meV depending on the material [16, 18, 19].The Hamiltonian (1) can be also derived from a micro-scopic theory. From a simplified tight-binding model [17],we obtain ∆ = 1 .
41 eV and t = 1 .
19 eV.In order to turn a TMDC into a TRI topological in-sulator, we need to consider the effect of strain on theband structure. Strain is described by a rank-2 tensor, u αβ = ( ∂ α u β + ∂ β u α + ∂h∂x α ∂h∂x β ), where u = ( u x , u y ) ( h )is the in-plane (out-of-plane) displacement of the unit cell Irrep TR Even TR Odd A (cid:48) , σ z , (cid:80) α u αα - A (cid:48) - τ z , s z E (cid:48) ( σ x , τ z σ y ), ( u xx − u yy , − u xy ) ( τ z σ x , σ y ), p = ( p x , p y ) E (cid:48)(cid:48) - ( s x , s y )TABLE I: Symmetry classification of the electronic operatorsand strain tensor components according to the irreducible rep-resentations (Irrep) of D h and time reversal operation. in the long-wave length limit. The three components of u αβ can be split according to irreducible representationsof D h as shown in Table I. Thus, the coupling with strainreads: H strain = β t (cid:88) α u αα + β t (cid:88) α u αα σ z + β t [( u xx − u yy ) σ x − u xy τ z σ y ] . (2)Microscopically, the origin of these couplings is thechange in the hybridization between d orbitals from M atoms and p orbitals from X atoms due to the dis-tortion of the lattice. In our microscopic tight-bindingmodel [17], the phenomenological constants β , , aregiven by the Gr¨uneisen parameters [11] associated to thehopping amplitudes considered in the calculation [17].The trace of the strain tensor generates scalar poten-tials of different strength in the valence and conductionbands. In addition, strain can be introduced in Eq. (1)as a minimal coupling p → p − e A to a vector potential e A = (cid:126) β a τ z ( u yy − u xx , u xy ). The presence of τ z in-dicates that the pseudo-magnetic field has opposite signon different valleys, which is necessary as strain does notviolate TRI.Henceforth, we assume that the system is doped withholes and therefore the Fermi level crosses the valenceband. Integrating out the conduction band to leadingorder in ∆ − yields: H v = − Π + ( τ z )Π − ( τ z )2 m ∗ + U ( r ) + λ SO s z τ z , (3)where m ∗ = ∆ / v (cid:39) . U ( r ) = g ( u xx + u yy )with g = t ( β − β ), and Π ± ( τ z ) = ( τ z p x ± ip y ) − e ( τ z A x ± iA y ), which obey the commutation relation[Π + ( τ z ) , Π − ( τ z )] = 2 e (cid:126) τ z B ( r ), where B ( r ) = ∂ x A y − ∂ y A x . Corrections of O (cid:0) ∆ − (cid:1) have been also obtainedand can lead to mixing of the Landau levels, but theycan be neglected as for typical parameters (cid:126) ω c (cid:46) − ∆.In the absence of strain (i.e. A = U = 0), Eq. (3)describes the Bloch states at the top of the valence bandnear K, K (cid:48) . Owing to the spin-orbit coupling ( ∝ λ SO ),for small hole doping (i.e. | (cid:15) F | (cid:28) λ SO ), the spin andvalley spin of the holes are locked to each other, i.e. onlyholes with either ( K, ↑ ) or ( K (cid:48) , ↓ ) can exist (cf. Fig. 2).Furthermore, the states at K and K (cid:48) are Kramers pairs. People
K K ' B - B H a L H b L K K 'LLs H c L K k kkk ' k ' k ' K ' kk kk ' k ' k ' FIG. 2: (Color Online) (a) Low-energy spectrum of a semicon-ducting transition metal dicalchogenide around the K and K (cid:48) points of the first Brioullin zone described by the continuum-limit Hamiltonian in Eq. (3). The inset shows the orientationof the stress tensor field discussed through the text with re-spect to the lattice. (b) Schematic representation of the Lan-dau Levels (LLs) induced in the valence band when strainis applied. (c) A superlattice can be also used to create atopologically non-trivial subband structure. A triangular su-perlattice produces replicas of K and K (cid:48) by mixing the va-lence and conduction band states at the mini-Dirac points asindicated by the arrows. This feature is crucial for the realization of a time-reversal invariant topological phase, as we argue below.Applying a pure shear strain, i.e. u xx + u yy = 0 (i.e. U = 0), and u xx = − u yy = − Cy , u xy = − Cx , we have B z ( r ) = ∂ x A y − ∂ y A x = − τ z B with B = C (cid:126) β ea > a K = Π − ( τ z = +1) / √ e (cid:126) B and a K (cid:48) =Π − ( τ z = − / √ eB : H = − (cid:126) ω c (cid:16) a † K a K + a † K (cid:48) a K (cid:48) (cid:17) + λ SO s z τ z , (4)where ω c = eB /m ∗ . As stated above, even for thelargest achievable pseudo-magnetic fields ( B ∼ T) (cid:126) ω c (cid:28) λ SO , and therefore, provided | (cid:15) F | (cid:28) λ SO , LLsat different valleys are occupied by holes with oppositespin. The resulting model has been shown by Bernevigand Zhang [12] to display the QSH effect, meaning thatthe Hall conductivity σ xyH is zero but the spin-Hall con-ductivity σ xysH is quantized in units of 2 e/ π .The strain configuration described above can be cre-ated by the methods described in Ref. [9]. For a TMDC2D crystal flake under trigonal strain [9] (cf. Fig. 2),one concern is the sample size. The latter is limited bythe maximum tensile strength of the TMDCs, T max . ForMoS , the thermodynamic relation σ αβ = 2 µu αβ ( σ αβ being the stress tensor and µ is the shear modulus) al-lows us to estimate the maximum sample size L ≈ T max µC .The values of µ (cid:39) . T max (cid:39) . L in µ m: B [ T ] ≈ /L [ µ m]. Using (cid:126) ω c /k B = 2 . B [ T ] and taking L ≈ µ m, we estimate (cid:126) ω c /k B (cid:39)
20 K for MoS .For small strained 2D crystal flakes, it is necessary totake into account the effect of an inhomogeneous pseudo-magnetic field resulting from a non-uniform strain dis-tribution. In this regard, we note that the lowest LLeigenfunctions are null eigenvectors of Π − ( τ z ):Π − ( τ z ) ψ ( r ) = 0 , (5)Therefore, following [32], we write A ( r ) = τ z ( ˆz × ∇ χ ( r ) + ∇ φ ( r )), which allows to solve (5)for K [ K (cid:48) ] and yields ψ ( r ) = f ( z ∗ ) e π Φ0 ( χ ( r )+ iφ ( r )) [ ψ ( r ) = f ( z ) e π Φ0 ( χ ( r ) − iφ ( r )) ], where f ( z ∗ ) [ f ( z )] is apolynomial of z = ( x + iy ) [ z ∗ = ( x − iy )] of maxi-mum degree N = [Φ / Φ ], Φ = (cid:82) d r B ( r ) > = h/e the flux quantum [33].Hence, the wave-function describing N ↑ = N ↓ = N (non-interacting) electrons in the lowest LL reads [32]:Φ ( { r iα } ) = e − F (cid:89) i 10 meV ∼ 100 K) [16], which also re-quires much lower temperatures. We hope that the feasi-bility of our proposal will stimulate further experimentalwork along these directions. We stress that comparedto other QSH systems [5, 6], strained 2D TMDCs allowfor much larger tunability of material parameters as wellas the strength of the Berry curvature responsible forthe Landau levels. In this regard, the main obstacle forthe observation of the QSH effect appears to be the re-duced carrier mobility in currently available 2D crystalsof TMDCs. However, given the notorious technologicalpotential of these materials, we expect this obstacle willbe overcome in the near future.MAC acknowledges financial support from a start-up fund from NTHU and NCS and NCTS (Taiwan).HO acknowledges financial support through a JAE-Pregrant (CSIC, Spain). FG acknowledges support fromthe Spanish Ministry of Economy (MINECO) throughGrant no. FIS2011-23713, the European Research Coun-cil Advanced Grant (contract 290846) and from Euro-pean Commission under the Graphene Flagship contractCNECT-ICT-604391. [1] Q. H. Wang, K. Kalantar-Zadeh, A. Kiss, J. N. Coleman,and M. S. Strano, Nature Nanotech , 699 (2012) andreferences therein.[2] B. Radisavljevic, A. Radenovic, J. Brivio, V. Giacometti,and A. Kis, Nature Nanotechnology, , 147 (2011);H. Nam, S. Wi, H. Rokni, M. Chen, G. Priessnitz, W.Lu, and X. Liang, ACS Nano, , 5870 (2013).[3] M. Z. Hassan and C. L. Kane, Rev. Mod. Phys. , 3045(2010); X.-L. Qi and S.-C. Zhang, ibid , 1057 (2011);B. A. Bernevig with T. L. Hughes, Topological Insulatorsand Topological Superconductors , Princeton Univ. Press,(Princeton 2013), and references therein.[4] See e.g. talk by S.-C. Zhang, available from [5] M. K¨onig, S. Wiedmann, C. Br¨une, A. Roth, H. Buhmann,L. W. Molenkamp, X.-L. Qi, S.C. Zhang, Science , 766(2007).[6] I. Knez, R.-R. Du, G. Sullivan, Phys. Rev. Lett. 107,136603 (2011)[7] C. Weeks, J. Hu, J. Alicea, M. Franz, and R. Wu, PhysicalReview X , 021001 (2011).[8] C. L. Kane and E. J. Mele, Phys. Rev. , 226801 (2013).[9] F. Guinea, M. I. Katsnelson, and A. K. Geim, NaturePhysics, , 30 (2010).[10] N. Levy et al. Science , 544 (2010); J. Lu, A.H. CastroNeto, and K. P. Loh, Nature Comm. , 823 (2011).[11] M. A. Vozmediano, M. I. Katsnelson, and F. Guinea,Physics Reports , 109 (2010).[12] B. A. Bernevig and S.-C. Zhang, Phys. Rev. Lett. ,106802 (2006).[13] D. Xiao, G.-B. Liu, W. Feng, X. Xu, and W. Yao, Phys. Rev. Lett. , 196802 (2012). , 045416 (2013).[14] E. Cappelluti, R. Rold´an, J.A. Silva-Guilln, P. Ordej´on,F. Guinea, Phys. Rev. B , 075409 (2013).[15] H. Ochoa and R. Rold´an, Phys. Rev. B , 245421(2013).[16] Z. Y. Zhu, Y. C. Cheng, and U. Schwingenschl¨ogl, Phys.Rev. B , 153402 (2011); W. Feng, Y. Yao, W. Zhu, J.Zhou, W. Yao, and D. Xiao, ibid , 165108 (2012); A.Korm´anyos, V. Z´olyomi, N. D. Drummond, P. Rakyta, G.Burkard, V. I. Fal’ko, ibid , 045416 (2013).[17] See supplementary material.[18] W. Jin et al. Phys. Rev. Lett. , 106801 (2013); S. K.Mahatha, K. D. Patel, abd K. S. R. Menon, J. Phys. Cond.Matt. et al. Phys. Rev. Lett. , 136805 (2010);W. S. Yun et al. Phys. Rev. B , 033305 (2012).[20] We take C > β ≈ et al. ACS Nano Lett. , 358 (2013).[21] R. C. Cooper, C. Lee, C. A. Marianetti, X. Wei, J. Hone,and J. W. Kysar, Phys. Rev. B , 035423 (2013). Thevalue of µ is obtained from the calculated Young’s modulus E and Poisson’s ratio ν as µ = E ν ) .[22] This is the effective mass inferred from the tight-binding model. The actual hole effective mass is slightlyanisotropic [16]. We take m ∗ = 0 . et al. Phys. Rev. B , 155330 (2008).[24] T. Low, F. Guinea, and K. I. Katsnelson, Phys. Rev. B , 054303 (2009).[26] F. D. M. Haldane, Phys. Rev: Lett. , 2015 (1988).[27] M. Yankowitz et al , Nature Phys. , 382 (2012).[28] L. A. Ponomarenko et al , Nature , 594 (2013).[29] C. R. Dean et al , Nature , 598 (2013).[30] B. Hunt et al , Science , 1427 (2013).[31] O. P. Sushkov, and A. H. Castro Neto, Phys. Rev. Lett. , 186601 (2013).[32] Y. Aharonov and A. Casher, Phys. Rev. A , 2461(1979); M. Bander, Phys. Rev. B , 9028 (1990).[33] Φ is assumed to be positive as in the uniform case. Notethat χ ( r ) = − (cid:82) d r (cid:48) π ln( | r − r (cid:48) | ) B ( r (cid:48) ) ∼ − Φ2 π ln( | r | ) as | r | → ∞ , and therefore ψ ( r ) ∼ f ( x + iτ z y ) | r | − ΦΦ0 .[34] C. Xu and J. Moore, Phys. Rev. B , 106401 (2006).[35] R. Rold´an, E. Cappelutti, and F. Guinea, Phys. Rev. B , 054515 (2013).[36] A. Shitade et al. , Phys. Rev. Lett. 102, 256403 (2009) SUPPLEMENTARY INFORMATIONMicroscopic derivation of the two-bands Hamiltonian Tight-binding model We consider a simplified tight-binding Hamiltonian acting on the subspace span by the symmetry-adapted Blochwave functions | M d z − r (cid:105) , (cid:12)(cid:12) M d x − y + τ id xy (cid:11) , | X (b) p y + iτ p x (cid:105) , and | X (b) p y − iτ p x (cid:105) . The model is strictly validaround K τ points, where τ = ± K ( τ = +1) and K (cid:48) ( τ = − (cid:104) M d z − r | H T B | M d z − r (cid:105) = ∆ (cid:10) M d x − y + τ id xy (cid:12)(cid:12) H T B (cid:12)(cid:12) M d x − y + τ id xy (cid:11) = ∆ (cid:104) X (b) p y ± τ ip x | H T B | X (b) p y ± τ ip x (cid:105) = ∆ p . (7)The hopping integrals sketched in Fig. 3 can be expressed in terms of the two-center Slater-Koster parameters V dpσ , V dpπ [1] as follows: t x = √ V pdσ cos ϕt y = − V pdπ cos ϕt z = cos ϕ (cid:18) sin ϕ − 12 cos ϕ (cid:19) V pdσ − √ ϕ sin ϕV pdπ , (8)where the angle ϕ is defined in the main text. The matrix elements between Bloch states of orbitals of M and Xatoms at (cid:126)q read: (cid:104) M d z − r | H T B | X (b) p y ± iτ p x (cid:105) = (cid:88) ˆ δ e − ia(cid:126)q · ˆ δ t z (ˆ y ± iτ ˆ x ) · ˆ δ (cid:10) M d x − y + τ id xy (cid:12)(cid:12) H T B | X (b) p y ± iτ p x (cid:105) = (cid:88) ˆ δ e − ia(cid:126)q · ˆ δ t x (ˆ y ± iτ ˆ x ) · ˆ δ × (cid:20) (cid:16) ˆ x · ˆ δ (cid:17) − iτ (cid:16) ˆ xy · ˆ δ (cid:17) + iτ − (cid:21) ++ (cid:88) ˆ δ e − ia(cid:126)q · ˆ δ t y (ˆ y ± iτ ˆ x ) · ˆ δ ⊥ × (cid:104) (cid:16) ˆ δ · ˆ x (cid:17) (cid:16) ˆ δ ⊥ · ˆ x (cid:17) − iτ (cid:16) ˆ δ · ˆ xy (cid:17) (cid:16) ˆ δ ⊥ · ˆ xy (cid:17)(cid:105) (9)where the unit vectors are defined as: ˆ x = (1 , y = (0 , xy = 1 √ , δ ⊥ = ( δ y , − δ x ) (10)and ˆ δ = ( δ x , δ y ) are the three unit vectors connecting nearest-neighbors in the honeycomb lattice:ˆ δ = (cid:40) ˆ y, (cid:32) √ , − (cid:33) , (cid:32) − √ , − (cid:33)(cid:41) (11)This model for the case of MoS leads to the electronic bands shown in Fig. 4. (cid:43) (cid:43)(cid:45)(cid:45)(cid:43)(cid:45) (cid:61) t x (cid:43) (cid:61) t y (cid:43)(cid:45) (cid:45)(cid:43)(cid:45) (cid:45)(cid:45)(cid:43)(cid:43)(cid:43)(cid:43)(cid:45)(cid:45) (cid:61) t z FIG. 3: Hopping integrals considered in the tight-binding model.Tight-binding parameter Value (eV)∆ -1.512∆ -3.025∆ p -1.276 V pdσ -2.619 V pdπ -1.396TABLE II: Tight-binding parameters considered in the calculation. The values for MoS are taken from Ref. 2. k · p theory At K τ points the Hamiltonian in first quantization reads as the matrix: ∆ t z 00 ∆ − t x + t y )3 t z p − t x + t y ) 0 ∆ p (12) Good approximation (cid:71) K M (cid:71) (cid:45) E (cid:72) e V (cid:76) FIG. 4: Bands calculated within the tight-binding model described in the text. The model is only valid around K τ points(highlighted in the figure). We take the values summarized in Table II for MoS . By diagonalizing this Hamiltonian we obtain the eigenvectors of the Bloch states maximally localized at d z − r and d x − y + iτ d xy respectively, which define the conduction and band states respectively: | ψ c (cid:105) = 1 (cid:113) | c | c | ψ v (cid:105) = 1 (cid:113) | v | v (13)where: c = − t z (cid:112) t z + (∆ p − ∆ ) + ∆ p − ∆ v = 6 ( t x + t y ) (cid:113) 36 ( t x + t y ) + (∆ p − ∆ ) + ∆ p − ∆ (14)Then, the k · p Hamiltonian is just: H k · p = (cid:32) (cid:104) ψ c | H T B | ψ c (cid:105) (cid:104) ψ c | H T B | ψ v (cid:105)(cid:104) ψ v | H T B | ψ c (cid:105) (cid:104) ψ v | H T B | ψ v (cid:105) (cid:33) (15)If we write (cid:126)q = K τ + (cid:126)k and expand in powers of (cid:126)k up to first order we get the Hamiltonian of Eq. (1) of the maintext with: ∆ = 12 (cid:18) ∆ − ∆ − (cid:113) tz + (∆ − ∆ p )) − (cid:113) 36 ( t x + t y ) + (∆ − ∆ p )) (cid:19) = 1 . 41 eV t = √ c ( tx − t y ) + vt z ]2 (cid:113) | c | (cid:113) | v | = 1 . 19 eV (16) Coupling with strain We repeat this calculation by considering the change in the hopping integrals t x,y,z due to the displacement of theatoms. For the hopping integral between atoms at sites α , β we have: t αβ → t αβ + ∂t αβ ∂ r ( u α − u β ) (17)where u α ( β ) is the displacement of the atom at site α ( β ). For small displacements we can approximate: ∂t αβ ∂ r ≈ √ a t αβ β αβ ˆ δ αβ u α − u β ≈ a √ δ αβ · (cid:126)∂ u (18)where β αβ ≡ − ∂ ln t αβ /∂ ln a is the Gr¨uneisen parameter associated to t αβ . Therefore, in the the previous matrixelements we have to consider now: t x,y,z → t x,y,z (1 + β x,y,z δ i δ j u ij ) (19)where u ij are the compnents of the strain tensor. By repeating the previous calculation at K τ we obtain the Hamil-tonian of Eq. (2) of the main text with: tβ = 32 (cid:32) cβ z t z | c | − v ( β x t x + β y t y )1 + | v | (cid:33) tβ = 32 (cid:32) cβ z t z | c | + v ( β x t x + β y t y )1 + | v | (cid:33) tβ = 3 ( cβ x t x − cβ y t y − vβ z t z )4 (cid:113) | c | (cid:113) | v | (20) Gaps induced by superlattice potentials and Haldane-Kane-Mele phase We next provide the details of how a superlattice can be used to induce the subbands conduction and valencebands of a doped TMDCs. We approximate the bands of the homogeneous system by the two-band continuum-limitHamiltonian introduced above, whose eigenvalues and eigenfunctions read: (cid:15) (cid:126) k = ± (cid:113) ∆ + ( v | (cid:126) k | ) (cid:12)(cid:12)(cid:12) (cid:126) k (cid:69) = (cid:32) cos (cid:0) θ k (cid:1) sin (cid:0) θ k (cid:1) e iφ k (cid:33) (21)where θ k = arctan[( v F | (cid:126) k | ) / ∆] and φ k = arctan( k y /k x ).A superlattice potential hybridizes states (cid:12)(cid:12)(cid:12) (cid:126) k (cid:69) and (cid:12)(cid:12)(cid:12) (cid:126) k + (cid:126) G (cid:69) where the vectors (cid:126) G define the superlattice. We considerthe six lowest vectors (cid:126) G , with | (cid:126) G | = (4 π ) / ( √ L ), where L = N a is the lattice constant of the N × N superlattice. Weassume that v | (cid:126) G | (cid:28) ∆ and that the superlattice potential, V (cid:126) G , is such that V (cid:126) G (cid:28) ( v | (cid:126) G | ) / (2∆), so that perturbationtheory in V (cid:126) G applies. We also assume that the lattice potential is sufficiently smooth, | (cid:126) G | (cid:28) | (cid:126) K + − (cid:126) K − | , where (cid:126) K ± are the corners of the Brillouin Zone of the TMDC lattice, and neglect intervalley scattering.Using first order perturbation theory, each set of three points at corners of the Brillouin Zone connected by super-lattice reciprocal vectors leads to a 3 × H κ,κ (cid:48) ≡ (cid:15) κ,κ (cid:48) ± ¯ vk x V κ,κ (cid:48) V ∗ κ,κ (cid:48) V ∗ κ,κ (cid:48) (cid:15) κ,κ (cid:48) ± ¯ v (cid:16) − k x + √ k y (cid:17) V κ,κ (cid:48) V κ,κ (cid:48) V ∗ κ,κ (cid:48) (cid:15) κ,κ (cid:48) ± ¯ v (cid:16) − k x − √ k y (cid:17) (22)where (cid:15) κ,κ (cid:48) = (cid:15) = ( v | (cid:126) K | ) / (2∆), ¯ v ≈ ( v | (cid:126) K | ) / ∆, and the two signs correspond to the κ and κ (cid:48) points. For V κ,κ (cid:48) = | V κ,κ (cid:48) | e iφ κ,κ (cid:48) , the the energies and eigenfunctions at the κ and κ (cid:48) points, in the basis used to write eq.(22)are: (cid:15) a = (cid:15) + 2 | V κ,κ (cid:48) | cos ( φ κ,κ (cid:48) ) | a (cid:105) = 1 √ | (cid:105) + | (cid:105) + | (cid:105) ) (cid:15) b = (cid:15) + 2 | V κ,κ (cid:48) | cos (cid:18) π φ κ,κ (cid:48) (cid:19) | b (cid:105) = 1 √ (cid:16) | (cid:105) + e πi/ | (cid:105) + e − πi/ | (cid:105) (cid:17) (cid:15) c = (cid:15) + 2 | V κ,κ (cid:48) | cos (cid:18) π φ κ,κ (cid:48) (cid:19) | c (cid:105) = 1 √ (cid:16) | (cid:105) + e − πi/ | (cid:105) + e πi/ | (cid:105) (cid:17) (23)In the case of φ κ,κ (cid:48) = 0 , π , then V κ,κ (cid:48) real, an expansion in powers of | (cid:126) k | shows that states | b (cid:105) and | c (cid:105) define aneffective 2 × v/ 2. For φ κ,κ (cid:48) = 0 this Dirac point gives the lowest edge of thehighest valence subband, and for φ κ,κ (cid:48) = π the upper edge of the lowest conduction band. The degeneracy of theseDirac points is lifted for complex values of V κ,κ (cid:48) . The problem is equivalent to a gapped Dirac equation with gap∆ κ,κ (cid:48) = 2 √ | V κ,κ (cid:48) | sin ( φ κ,κ (cid:48) ). If φ κ (cid:48) and φ κ (cid:48) have different signs, the two gaps also have opposite signs, leading to alowest subband with a Chern number equal to one. This is a realization of Haldane’s model [3].0 Κ (cid:162) Κ (cid:162) Κ (cid:162) Κ Κ Κ G G G G G G FIG. 5: (Color online). Sketch of the effect of the Brillouin Zone in a superlattice, and states mixed by the superlatticepotential. The superlattice potential is a 2 × f ( r ) = (cid:88) m =0 ... e i G m · r f ( r ) = i (cid:88) m =0 ... ( − m e i G m · r (24)Then, we can construct the inversion-symmetric superlattice potentials as: V s = v | G | ∆ s f ( r ) V m = v | G | ∆ m f ( r ) σ z V g = v ∆ g ( σ x , τ z σ y ) · (ˆ z × ∇ ) f ( r ) (25)The coefficients ∆ s,m,g in these expressions are dimensionless phenomenological constants with the energy scale setby v | G | = πt √ N .The scalar potential has matrix elements: (cid:104) k + G m | V s ( G m ) | k (cid:105) = v | G | ∆ s (cid:20) cos (cid:18) θ k + G m (cid:19) cos (cid:18) θ k (cid:19) + sin (cid:18) θ k + G m (cid:19) sin (cid:18) θ k (cid:19) e i ( φ k − φ k + G m ) (cid:21) (26)Equivalently for the mass potential: (cid:104) k + G m | V m ( G m ) | k (cid:105) = v | G | ∆ m (cid:20) cos (cid:18) θ k + G m (cid:19) cos (cid:18) θ k (cid:19) − sin (cid:18) θ k + G m (cid:19) sin (cid:18) θ k (cid:19) e i ( φ k − φ k + G m ) (cid:21) (27)The edges of the first subband are determined by the shifts in the energies of the corners of the superlattice BrillouinZone. We assume that | k | = | k + G m | = κ = (4 π ) / (3 N a ). Then: (cid:104) k + G m | V s,m ( G m ) | k (cid:105) ≈ v | G | ∆ s,m (cid:20) − v κ (cid:16) ∓ e i ( φ k − φ k + G m ) (cid:17)(cid:21) (28)For the gauge potential we have: (cid:104) k + G m | V g ( G m ) | k (cid:105) = i ( − m v | G | ∆ g (cid:20) cos (cid:18) θ k + G m (cid:19) sin (cid:18) θ k (cid:19) e i ( φ k − φ G m ) − sin (cid:18) θ k + G m (cid:19) cos (cid:18) θ k (cid:19) e i ( φ G m − φ k + G m ) (cid:21) ≈ ( − m vκ ∆ v | G | ∆ g e i φ k − φ k + G m cos (cid:18) mπ − φ k + φ k + G m (cid:19) (29)1The same can be done with the inversion-asymmetric superlattice potentials, defined as:˜ V s = v | G | ˜∆ s f ( r )˜ V m = v | G | ˜∆ m f ( r ) σ z ˜ V g = v ˜∆ g ( σ x , τ z σ y ) · (ˆ z × ∇ ) f ( r ) (30)By repeating the same calculation we obtain: (cid:104) k + G m | ˜ V s,m ( G m ) | k (cid:105) ≈ i ( − m v | G | ˜∆ s,m (cid:20) − v κ (cid:16) ∓ e i ( φ k − φ k + G m ) (cid:17)(cid:21) (cid:104) k + G m | ˜ V g ( G m ) | k (cid:105) ≈ − i vκ ∆ v | G | ˜∆ g e i φ k − φ k + G m cos (cid:18) mπ − φ k + φ k + G m (cid:19) (31)From this analysis it is clear that inversion-asymmetric potentials are needed in order to induce a topologicalsubband structure. To the leading order in vκ/ ∆, considering scalar potentials only, we have: V κ = (cid:104) κ = κ + G | V s ( G ) + ˜ V s ( G ) | κ (cid:105) = v | G | (cid:16) ∆ s + i ˜∆ s (cid:17) V κ (cid:48) = (cid:104) κ (cid:48) = κ (cid:48) + G | V s ( G ) + ˜ V s ( G ) | κ (cid:48) (cid:105) = v | G | (cid:16) ∆ s − i ˜∆ s (cid:17) (32)So V κ,κ (cid:48) = v | G | (cid:113) ∆ s + ˜∆ s × e ± i arctan (cid:16) ˜∆ s ∆ s (cid:17) . Hence, the gap is ∆ κ,κ (cid:48) = ± √ v | G | ˜∆ s . Either the lowest conductionor valence subbands derived from the bands at a given spin polarized valley in the band structure of the TMDC havea Chern number equal to one. This Chern number is compensated by the opposite value from the other valley due totime-reversal symmetry. Therefore, this system is effectively a realization of the Kane-Mele model [5]. [1] J. C. Slater and G. F. Koster, Phys. Rev. , 1498 (1954).[2] E. Cappelluti, R. Rold´an, J.A. Silva-Guill´en, P. Ordej´on, F. Guinea, Phys. Rev. B , 075409 (2013).[3] F. D. M. Haldane, Phys. Rev. Lett. , 2015 (1988).[4] J. R. Wallbank, A. A. Patel, M. Mucha-Kruczyn’ski, A. K. Geim, and V. I. Fal’ko, Phys. Rev. B , 245408 (2013).[5] C. L. Kane and E. J. Mele, Phys. Rev. Lett. , 226801, ibidibid