This is the submitted version of the article: Cysne T.P., Costa M., Canonico L.M., Nardelli M.B., Muniz R.B., Rappoport T.G.. Disentangling Orbital and Valley Hall Effects in Bilayers of Transition Metal Dichalcogenides. Physical Review Letters, (2021). 126. 056601: - . 10.1103/PhysRevLett.126.056601. Available at: https://dx.doi.org/10.1103/PhysRevLett.126.056601 Disentangling orbital and valley Hall eects in bilayers of transition metal dichalcogenides Tarik P. Cysne,1 Marcio Costa,1 Luis M. Canonico,2 M. Buongiorno Nardelli,3 R. B. Muniz,1 and Tatiana G. Rappoport 4, 5 1 Instituto de Física, Universidade Federal Fluminense, 24210-346 Niterói RJ, Brazil ∗ 2 Catalan Institute of Nanoscience and Nanotechnology (ICN2), CSIC and BIST, Campus UAB, Bellaterra, 08193 Barcelona, Spain 3 Department of Physics and Department of Chemistry, University of North Texas, Denton TX, USA 4 Instituto de Telecomunicações, Instituto Superior Tecnico, University of Lisbon, Avenida Rovisco Pais 1, Lisboa, 1049001 Portugal 5 Instituto de Física, Universidade Federal do Rio de Janeiro, C.P. 68528, 21941-972 Rio de Janeiro RJ, Brazil (Dated:) It has been recently shown that monolayers of transition metal dichalcogenides (TMDs) in the 2H structural phase exhibit relatively large orbital Hall conductivity plateaus within their energy band gaps, where their spin Hall conductivities vanish [1, 2]. However, since the valley Hall eect (VHE) in these systems also generates a transverse ow of orbital angular momentum it becomes experimentally challenging to distinguish between the two eects in these materials. The VHE requires inversion symmetry breaking to occur, which takes place in the TMD monolayers, but not in the bilayers. We show that a bilayer of 2H-MoS 2 is an orbital Hall insulator that exhibits a sizeable OHE in the absence of both spin and valley Hall eects. This phase can be characterised by an orbital Chern number that assumes the value CL = 2 for the 2H-MoS 2 bilayer and CL = 1 for the monolayer, conrming the topological nature of these orbital-Hall insulator systems. Our results are based on density functional theory (DFT) and low-energy eective model calculations and strongly suggest that bilayers of TMDs are highly suitable platforms for direct observation of the orbital Hall insulating phase in two-dimensional materials. Implications of our ndings for attempts to observe the VHE in TMD bilayers are also discussed. Introduction: The orbital Hall eect (OHE) is the or- It is noteworthy though that the VHE manifests only bital analog of the spin Hall eect and consists in the in the absence of inversion symmetry, which naturally appearance of a transverse current of orbital angular happens for the monolayers, but for bilayers, comprising momentum that is induced by a longitudinally applied two monolayers rotated by π with respect to each other, electric eld [3]. Recently, a renewed interest in orbital the inversion symmetry is restored. This substantially magnetism and other orbital eects [47] gave origin to aects valley related phenomena [2325]. For instance, various theoretical studies on the OHE and related phe- the valley Hall conductivity has opposite signs in each nomena [1, 2, 815]. The possibility of using the OHE layer, cancelling the VHE for the bilayer [2628], as we to generate orbital torque in magnetic materials [16, 17] shall subsequently discuss. Nevertheless, it is also possi- motivated new experimental works on orbital dynamics ble to break inversion symmetry in the bilayers by apply- in magnetic multilayers [18, 19], raising expectations that ing an electric eld perpendicular to the layers, by means orbital angular degrees of freedom may be eventually em- of which one can control the valley polarisation [29] and ployed to process information in logic and memory de- the VHE intensity [25] with a gate voltage. vices. Here we perform calculations of the orbital Hall con- The interrelation between the OHE and the pres- ductivities for ultra-thin lms (single layer and bilayer) of ence of orbital textures in reciprocal space [10] has 2H-MoS2 which is representative of this class of systems. been established and characterised both theoretically We combine Density Functional Theory (DFT) and an and experimentally in several low-dimensional materi- eective low-energy model to disentangle the valley and als [1, 7, 12, 13, 20, 21], widening the class of systems orbital physics of TMD bilayers and explore some of their that may be utilised for orbitronic applications. More topologic orbital features. specically, the occurrence of relatively large OHE has Implications of our ndings regarding interpretations been predicted in the 2H structural phase of TMD mono- of recent experiments on the electric control of the VHE layers [1, 2], where it is associated with the presence of in MoS2 bilayers [24, 25] are also briey discussed. Our a Dresselhaus-like orbital texture around the valleys [1]. results strongly indicate that bilayers of TMDs consti- However, it is experimentally challenging to observe just tute a fertile play-ground for exploring orbital angular the OHE in 2H-TMD monolayers due to the concurrent momentum current generation in 2D-like systems. presence of the VHE that also contribute to the transport DFT results: Our DFT calculations [30, 31] were of orbital angular momentum in these systems [22]. performed with the plane-wave-based code Quantum 2 Espresso [32]. The exchange and correlation poten- Figures 1 (a) and (b) illustrate results of our calcula- tial is treated within the generalised gradient approx- tions for a monolayer and for a bilayer of 2H-MoS 2 , re- imation (GGA) [33]. The ionic cores were described spectively. In the left panels of Figure 1 we compare the with fully relativistic projected augmented wave (PAW) band structures obtained from DFT (purple solid lines) potentials [34]. We used a cuto energy of 63 Ry for and from Paoflow (yellow dashed lines). The agree- the wavefunctions and a value 10 times larger for the ment between the two approaches is excellent. For the charge density. In order to reproduce the interlayer monolayer we obtain a direct energy-band gap of 1.60 eV, distance of the MoS2 bilayer we have used the DFT- whereas for the bilayer we found an indirect gap of 1.28 D3 [35] method, which describe reasonably well the van eV, in agreement with previous calculations [45]. The re- der Walls forces in these systems. We have chosen a sults for the SHC (red solid lines) and for the OHC (blue 10×10×1 reciprocal space sampling, and to avoid spuri- lines) are shown in right-hand side panels of Figure 1. In ous interaction due to periodic boundary conditions we accordance with our previous results for the monolayer insert a vacuum spacing of 15Å. We constructed an ef- [1] we note in right panel of Figure 1 (a) that the SHC fective tight-binding Hamiltonian from our DFT calcula- vanishes in the main energy gap, but the OHC is nite tions using the pseudo atomic orbital projection (PAO) and exhibits relatively high plateau of ≈ 2.6 (e/2π ) in method [36, 37]. The PAO method consists of projecting height within this energy range. For the bilayer, how- the DFT Kohn-Sham orbitals into the compact subspace ever, the right-hand side panel of panel Fig. 1(b) show spanned by the pseudo atomic orbitals which are natu- that the height of OHC plateau is essentially twice the rally built-in into the PAW potentials. The PAW poten- monolayer value, while the SHC remain null in the main tials used for the Mo and S were constructed with a sspd energy gap because it is topologically trivial. Although and sp basis, respectively. other regions in the B. Z. contribute to the OHC [2], Once the PAO Hamiltonian is constructed we can cal- the main contribution originates from the orbitally pro- culate the spin Hall (SH) and orbital Hall (OH) responses jected Berry curvature in the vicinity of K and K', as il- to an applied electric eld [1, 3, 8, 13, 3840]. Up to linear lustrated in the supplementary material [46]. To explore order on the external eld they are given by: the physics behind these results, it is instructive to make use of a low-energy approximation around the K -points (valleys) of the BZ, to build a simple model that is able Z η e X X σOH(SH) = d2 kfn~k Ω η~ , (1) (2π)2 n BZ n,k to reasonably describe the main transport characteristics of these systems. η where σOH(SH) is the orbital Hall (spin Hall) DC con- Low energy calculations: Similarly to the monolayers, ductivity with polarisation along the η -direction, and the low-energy physics of TMD bilayers is dominated by the dz2 , dx2 −y2 and dxy atomic orbitals of the transition Xη metal atoms [26, 47, 48]. We follow references 27 26 to ψm,~k vx (~k) ψn,~k " # ψn,~k j ψm,~k X y,~ k build a simplied tight-binding (TB) model Hamiltonian X Ω η~ = 2~ Im (2) n,k m6=n (En,~k − Em,~k + i0+ )2 in reciprocal space, which is expanded up to rst order in the electronic momentum around the valleys located represents the angular-momentum-weighted Berry cur- at K~ = (4π/3a)x̂ and K ~ 0 = −K ~ . This procedure leads vature [10, 40]. Here, En,~k denotes the eigenvalue of to the following Hamiltonian: the Hamiltonian H(~k) in reciprocal space, and |ψ ~ is n,k   the corresponding eigenvector; n is the band index, ~k is ∆ γ+ 0 0 γ− −τ sz λ 0 t⊥  the wave vector. The velocity operators are dened as H̃(~qτ ) =  , (3) 0 0 ∆ γ−  vx(y) (~k) = ∂H(~k)/∂~kx(y) , where x and y specify the Cartesian axes, and we assume that the electric eld 0 t⊥ γ+ τ s z λ is applied along the x̂ direction. The current density where γ± = at(τ qx ± iqy ), τ = ±1 is the valley quantum operator component along ŷ with polarisation η is de- number associated with valleys K and K 0 , respectively. X ned as j ~η = Xη vy (~k) + vy (~k)Xη /2, where for the  y,k Here, ~k = ~q + τ K ~ where ~q represents the wavevector rela- SH conductivity Xη = ŝη and for the OH conductivity tive to valleys and sz denotes the usual Pauli matrix. For (OHC) Xη = `ˆη ; ŝη and `η represent the η -components a 2H-MoS2 bilayer, an archetypal TMD, ∆ = 1.766eV of the spin and of the atomic angular momentum opera- is the monolayer band-gap, a = 3.160Å is the lattice tors, respectively. This is implemented in the Paoflow constant, t = 1.137eV is the intra-layer nearest-neighbor code [41] that has been successfully used to study topo- hopping, λ = 0.073eV is the spin-orbit coupling, and logical materials [42, 43], and time dependent spin dy- t⊥ = 0.043eV is the interlayer hopping [26]. namics [44] among other topics. For our conductivity cal- The TB basis for this minimal √ model com- culations we have increased the sampling to 200 ×200×1 prises { d1z2 , d1x2 −y2 − iτ d1xy / 2, d2z2 , d2x2 −y2 +  √ k-points in the 2D B.Z. iτ d2xy / 2}, where the superscripts 1 and 2 specify the 3 (a) 2 20 1, k + 2, k (a) (b) 2, k 10 1, k k (Å 2 ) 0 0 10 Energy (eV) 20 L + L −2 30 (c) (d) 1, k L 2, k 2, k L L (Å 2 ) 20 1, k k 10 −4 0 SHC 1, k + 2, k OHC 30 (e) (f) 1, k + 2, k L L −6 k (Å 2 ) Γ M K Γ −2 0 2 4 6 20 σxy z (e/2π) 10 K K0 0 1.0 0.5 0.0 0.5 1.0 1.0 0.5 0.0 0.5 1.0 (b) 2 q(Å 1) q(Å 1) Figure 2: Berry curvature Ωn,k at points K (a) and K 0 (b), 0 orbital weighted Berry curvature ΩL n,k at points K (c) and K 0 Energy (eV) (d) for the two valence bands E1 (q) and E2 (q), associated to SHC the two layers. (e) Total Berry curvature Ω1,k + Ω2,k and (f) −2 OHC orbital weighted Berry curvature ΩL 1,k + Ω2,k for the bilayer L TMD. −4 both E1 and E2 calculated around the K and K 0 points, −6 respectively. The Berry curvature for E2 has a positive Γ M K Γ −5 0 5 10 peak at K and a negative peak at K 0 , which gives rise σxy z (e/2π) to a VHE. The opposite occurs for the Berry curvature of E1 , which has a negative peak around K and a posi- Figure 1: Energy band structures (left panel) together with the spin-Hall and orbital-Hall conductivities (right panel) cal- tive peak at K 0 , giving origin to a VHE with an inverted culated for a MoS 2 monolayer (a) and for a MoS 2 bilayer (b). sign. By adding the contributions of both layers, the net The purple solid and yellow dashed lines depict the DFT and Berry curvature is zero in both valleys, and the VHE Paoflow band structure calculations, respectively. The hor- vanishes. This is a consequence of time-reversal symme- izontal blue dashed line shows the Fermi level. try and the presence of spatial inversion symmetry in the bilayer [27, 47, 49]. A similar situation occurs for TMDs with the T and T 0 structural phases, such as WTe 2 [50]. two layers of the bilayer, respectively. It is noteworthy Figures 2 (c) and (d) show the orbital weighted Berry cur- that the orbital angular momentum (OAM) operator in vatures for both bands around the K and K 0 points, re- this representation is given by Lz = diag(0, −2~τ, 0, 2~τ ), spectively. In contrast with the previous case, the peaks which clearly does not commute with the Hamiltonian of the orbital-weighted Berry curvatures for both bands dened in Eq.(3). have the same sign around both valleys. Hence, the to- Eq. (2) can be used with the four-band low-energy tal orbital-weighted Berry curvature has a nite value, Hamiltonian given by Eq. (3) to dene the Berry and the which leads to an OH insulating phase [1] with no VHE, orbital-weighted Berry curvatures that encode informa- as Figures 2 (e) and (f) illustrate. We note that in order tion of the VHE and OHE, respectively. For simplicity, to assess just the OHE it is crucial from the experimen- we shall initially neglect the eect of spin-orbit coupling tal point of view to have OHE without VHE, because the (λ), thereby restricting Eq. (3) to a spinless Hamilto- VHE also leads a transverse angular momentum current nian, and including a degenerescence factor gs = 2. Eq. [22, 51, 52] that is hard to be distinguished from the one (2) for the orbital weighted Berry curvature may also be generated by the OHE, as it happens for TMD monolay- employed to calculate the usual Berry curvature Ωn,k , ers [1]. Thus, our results show that bilayers of 2H-TMDs provided that Xη is replaced by ~1. The spinless Hamil- are very promising candidates for observing the orbital tonian generates two valence bands ( E1 (q) and E2 (q)) Hall insulating phase with no interferences from VHE or that can be regarded as arising from each of the TMD SHE. layers because of the relatively small interlayer hopping. Let us now address the topological characterisation of Figures 2 (a) and (b) present the Berry curvatures for the OH insulating phase in TMD-bilayers. Our Berry 4 (3). Figure 3 (a) show the eigenvalues of Lv (~k) calculated as functions of q . We clearly see that the spectrum splits in two separated sectors, allowing us to use the eigen- states of Lv (~k) in each valley Φ± q ) to calculate the n,τ (~ Chern numbers: Z 1 X CL± = d2 q ± Fn,τ (q), (4) 2π n,τ where Fn,τ ± (q) = −2Im[ ∂qx Φ± q ) ∂qy Φ± n,τ (~ q ) ]. Fig. 3 n,τ (~ (b) shows the integrands of Eq. (4). Since they have azimuthal symmetry, the calculations of CL± involve nu- Figure 3: (a) Positive (orange) and negative (blue) eigen- merical integrations of one-dimensional radial functions values  of the matrix Lv (~k) calculated as functions of the only. Our results for the insulating phases of the 2H- wavevector amplitudeP q relative to valleys. (b) Integrands of MoS2 bilayer and single layer are CL = 2 and CL = 1, Eq. (4) I ± (q) = n,τ Fn,τ (q) calculated as functions of q . ± respectively, supporting the idea that the relatively weak Band structures of zigzag nanoribbons with 14.8Å in breadth interlayer hopping in the bilayer makes it behave approx- calculated with the Paoflow Hamiltonian for a monolayer imately as a mere superposition of its two constituent (c) and for the spin-up sector of a bilayer (d) of 2H-MoS 2 . monolayers, which are rotated by π with respect to each The color code indicates the orbital angular momentum ex- pectation value hLz i. For better visualization, in the bilayer other. case only the spin-up bands are showed. The spin-down band- The existence of a nontrivial orbital Chern number structure are showed in the supplementary material [46]. should lead to the appearance of edge states when the bulk material is cut to form a ribbon. It is well known that zigzag TMD ribbons present crossing edge-states curvature analysis suggests that it is possible to asso- with interesting orbital properties, even though Z2 = 0. ciate an orbital Chern number to describe the distinctive [58]. Figures 3 (c) and (d) show the energy band spec- nature of these states in analogy with the well known tra of 2H-MoS2 zigzag nanoribbons, calculated with the spin Chern number [53, 54]. Here the situation is slightly use of Paoflow Hamiltonian for a monolayer and a bi- more subtle because the operator Lz does not commute layer including the orbital angular momentum expecta- with the Hamiltonian of Eq. (3) for nite ~q. This is sim- tion value hLz i(k) for each eigenstate. [59]. The energy ilar to the problem of a quantum spin Hall insulator in band spectrum for a monolayer ribbon depicted in Figure the presence of a Rashba SOC. To address this issue, we 3 (c) clearly shows two pairs of orbitally polarized intra- follow the procedure developed in Refs. 5557 to dene valley edge states [60] - one for each spin sector- which the orbital Chern number CL for the insulating phase of is compatible with the orbital Chern number CL = 1. Hamiltonian given by Eq. (3). In this formalism, CL = Results for the bilayer ribbon are displayed in Figure 3 (CL+ − CL− )/2 where CL± are the Chern numbers calculated (d), where we see two pairs of intra-valley edge states per with the eigenstates of an spin-sector - which is also compatible with CL = 2. For  OAM operator projected on the bilayer, the presence of inversion symmetry is trans-  the valence-band states L (k) = P (k)Lz P (k) , where v ~ ~ ~ lated in the existence of positive and negative Lz edge P (~k) is the projector operator. If the bands have or- states in both valleys. For clarity, Figure 3 (d) presents bital polarisation, the spectrum of Lv (~k) consists of two the results for spin-up while the SM presents the two groups of eigenvalues ( ) associated with ml = ±2 that components. are symmetrically separated by a gap. The projectors Experimental signatures: Let us now briey discuss on the eigenstates associated with the positive and nega- the experimental signatures of the OHE in TMD bilay- tive eigenvalues can then be used to calculate the Chern ers. Typically, to characterise the OHE in these mate- numbers CL± . rials, one needs the same experimental setups conceived Thus, to calculate CL , it is necessary to decompose the to analyse VHE in TMD bilayers, where inversion sym- valence-band states into two sectors with respect to op- metry breaking is induced by a gate voltage [25, 61]. For erator Lz . For that purpose, we rst obtain the matrix the bilayer in this case, both the OHE and the VHE lead Lv (~k), with matrix elements given by ψn,~k Lz ψm,~k , to magnetic moment accumulation at the sample's edges. where n, m label the valence-band eigenstates of the low- To provide some insights into what should be expected energy Hamiltonian; more details are given in the ac- in such experiments, we include a gate potential in in Eq. companying Supplementary Material [46]. It is worth (3) given by HU = diag(U, U, −U, −U ). mentioning that hereafter we reinstate the spin degree of For a nite U , the inversion symmetry is broken in freedom and the spin-orbit interaction in the Hamiltonian the bilayer and the VHE takes place. The OH and 5 VH conductivities can be calculated using Eqs.(1-3). explain the unexpected signals at zero bias in the bilay- To calculate the VH conductivity (VHC), we substi- ers. Careful experimental analysis of Kerr rotation and tute the integrand of Eq.(1) by Ωn,~q /~, and rewrite non-local resistance measurements as functions of gate σV H = (στ =+1 − στ =−1 ). Fig. 4(a) shows our results voltage may help to distinguish between the orbital and for the OHC (σOH z ) and VHC (σV H ) calculated as func- valley Hall eects in these materials. The results illus- tions of the Fermi energy ( EF ) for positive and negative trated in Fig. 4, in light of the experiments reported in values of U . There are clear diferences between the two Refs. [25, 61], suggest that ultrathin lms of TMDs are quantities. While the OHC is an even function of EF , promising platforms to explore the OHE in 2D materials. the VHC is odd. Also, the OHE is an even function of Final remarks and conclusion: Our DFT calculations U , whereas the VHE is odd - the valley magnetic mo- showed that centrosymmetric two-dimensional materials, ment inverts when U changes sign [24]. Panels (b) and such as a bilayer of 2H-MoS 2 , can host an orbital Hall in- (c) of Fig. 4 show the energy spectrum for U = 0 and sulating phase in the absence of both spin and valley Hall U = 0.2 eV, respectively. It is clear that U produces a eects. Using MoS2 as a prototype of the TMD family, rigid energy-band shift for the two layers, without chang- we have also unveiled the topological nature of OHE in ing their orbital polarisations. The OHE should remain these systems and calculated the orbital Chern numbers unchanged for small variations of EF , but decreases when for 2H-TMDs. Our work claries the interplay between EF crosses any band, as individual bands in each valley orbital and valley Hall conductivity in bilayer TMDs. We contribute to the total Chern number. We have also per- found that, in the absence of a gate voltage between the formed DFT calculations for the OHE in the presence of layers, the magnetic moment accumulation observed in an electric eld applied perpendicularly to the layers [62]. experiments should be dominated by the OHE as VHE is The results support our low-energy analysis and are pre- zero in centrosymmetric materials. For nite bias, OHE sented in the SM [46]. and VHE are still decoupled and can behave as compet- ing eects. We acknowledge CNPq/Brazil, CAPES/Brazil, FAPERJ/Brazil and INCT Nanocarbono for nancial support. TGR acknowledges funding from Fundação para a Ciência e a Tecnologia and Instituto de Tele- comunicações - grant number UID/50008/2020 in the framework of the project Sym-Break. She thankfully acknowledges the computer resources at MareNostrum and the technical support provided by Barcelona Super- computing Center (FI-2020-2-0033). MC acknowledge the National Laboratory for Scientic Computing (LNCC/MCTI, Brazil) for providing HPC resources of the SDumont supercomputer. MBN also acknowl- edges the High Performance Computing Center at the University of North Texas and the Texas Advanced Computing Center at the University of Texas, Austin, Figure 4: (a) Orbital (blue line) and valley Hall (orange line) for computational resources. LMC is supported by conductivities as a function of the Fermi energy EF for U=0.2 Project MECHANIC (PCI2018-093120) funded by the eV (solid line) and U=-0.2 eV (dashed line). Schematic rep- Ministerio de Ciencia, Innovación y Universidades. ICN2 resentation of the low energy spectrum for U=0 eV (b) and U =0.2 eV (c). The orbital polarisation of the top valence is funded by the CERCA Programme/Generalitat de bands is indicated by ml = ±2 and the solid orange and Catalunya and supported by the Severo Ochoa Centres dashed blue lines indicate the two spin orientations. of Excellence program, funded by the Spanish Research Agency (grant number SEV-2017-0706). Kerr rotation microscopy experiments[25] showed that a bilayer of MoS2 exhibits a sizeable Kerr rotation even in the absence of an applied gate voltage. It was argued that this unexpected behaviour could originate from substrate induced inversion symmetry breaking. Recent non-local resistance measurements in hBN encapsulated bilayer of MoS2 also exhibited non-local signal at zero gate volt- age [61]. The interpretation was the same, although one should not expect hBN to cause such a large inversion symmetry breaking eect. On the other hand, the OHE could be the source of this experimental evidence and 6 Supplementary material for Disentangling orbital and valley Hall eects in bilayers of transition metal dichalcogenides MOS 2 ZIGZAG NANORIBBONS BAND STRUCTURE - PAOFLOW Figure 5: Band structure of a zigzag nanoribbon of 2H-MoS 2 bilayer with width of 14.8Åin breadth calculated with the use of the Paoflow Hamiltonian. Panels (a) and (b) depict the results for the ↓ and ↑ spin components, respectively. The color code indicates the orbital angular momentum expectation value hLz i for the ↓ spin bands. MOS 2 ZIGZAG NANORIBBONS BAND STRUCTURE - 3 BANDS APPROXIMATION (a) (b) 2 2 Energy (eV) Energy (eV) 1 0 0 -3.1π 0 3.1π -π a 0 πa k k Figure 6: Zigzag nanoribbon band-structures calculated for a monolayer (a) and for a bilayer (b) of 2H-MoS 2 using a simplied 3-bands approximation. Solid orange and dashed blue lines indicate the two spin orientations. Zigzag TMD ribbons exhibit crossing edge-states that can also be modelled by an eective three-bands model, which captures their essential transport features [48]. Figures 6 (a) and (b) illustrate the energy band spectra of 2H-MoS 2 zigzag nanoribbons, calculated for a monolayer and for a bilayer, respectively [63, 64], using the three-bands model of Ref. 48 and nearest neighbour interlayer hopping integrals only. It is noteworthy that this simplied model does not capture the behaviour of high and low energy bands, because it does not take into account the orbitals dxz and dyz of the transition metal, and treats the eects of the chalcogens (S) perturbatively only. For this reason, the edge-states generated by this model are separated from the bulk valence bands by a non-realistic energy gap. In more realistic descriptions of a TMD nanoribbon [58, 65], this gap is lled by states with orbital character that are not considered in this simplied model. Nevertheless, the simplied model describes very well the nature of edge-states near their crossing. 7 AB INITIO CALCULATIONS IN THE PRESENCE OF AN APPLIED ELECTRIC FIELD                          Figure 7: (a) Lateral view of the MoS 2 bilayer structure. The red arrow reprsents the electric eld E = Ez ẑ . (b) MoS 2 monolayer energy band structure (left panel) together with the orbital-Hall conductivity (OHC) (right panel) calculated for Ez = 0.0 V/Å (red) and Ez = 0.4 V/Å (blue). The horizontal blue dashed line shows the Fermi level.                                             Figure 8: Energy band structures (left panel) together with the spin-Hall and orbital-Hall conductivities (right panel) calculated for a bilayer of MoS 2 for (a) Ez = 0.0 V/Å and (b) Ez = 0.2 V/Å. The purple-solid and yellow-dashed lines depict the DFT and Paoflow band structure calculations, respectively. The horizontal blue-dashed line denotes the Fermi energy. We examine the eects on the electronic structures of the TMD ultrathin lms caused by a homogeneous electric eld E = Ez ẑ [62], applied perpendicularly to the lm layers, using the modern theory of polarization [66, 67]. A non local energy functional is dened as the regular energy functional subtracted by the product Ez Pz , where Pz is is the polarization component along the ẑ -direction. Once the charge density and wave functions are converged we obtain the paoflow Hamiltonian and calculate the electronic properties of interest. In Fig. 7 we show the band structure together with the OHC of a monolayer of MoS 2 calculated for Ez = 0 and Ez =0.4 V/Å. We note that the MoS2 monolayer band gap and OHC are virtually not aected by the applied electrical eld, which is in agreement with previous results [68]. Fig. 8 shows the energy band structures together with the spin-Hall and orbital-Hall conductivities of a bilayer of MoS2 calculated for Ez = 0.0 V/Åand Ez = 0.4 V/Å. We note that the perpendicularly applied electric eld breaks the inversion symmetry between the two layers causing a substantial band gap reduction as Ez = 0.4 increases, whereas the OHC remains practically unchanged. We expect this behavior not to change as long as the band gap remais nite. For applied electric eld intensities 0.0 < Ez < 0.4 V/Å, the bilayer band energy gap varies almost linearly with Ez , as Fig. 9(a) illustrates, in agreement with previous DFT calculations [68]. Figure. 9(b) shows the calculated OHC in the same range of electrical eld. There is small reduction, ≈ 5%, in the OHC for Ez < 0.4 V/Å. 8 (a) 40 (b) 1.25 6 Band Gap (eV) 1.00 4 (e/2π) 0.75 σxy z 2 0.50 0.25 0 0.0 0.0 0.2 0.4 0.2 0.4 Ez (V /Å) Ez (V /Å) Figure 9: Evolution of the (a) band gap and (b) orbital-Hall conductivity of the MoS 2 bilayer with an applied perpendicular electrical eld Ez . ANGULAR MOMENTUM WEIGHTED BERRY CURVATURE - PAOFLOW To show the validity of the low-energy model to calculate the OHC in bilayer TMDs, we can inspect the contribution of the dierent high symmetry points to the integrand Q(~k) of σOHz : Z X z σOH ∝ Q(~k)d2 k, Q(~k) = Θ(En − EF ) ΩLz~ . (5) n,k BZ n In Fig. 10 we show Q(~k) along the high symmetry lines of the Brillouin zone. It has a sharp strong peak in K , showing that although the inuence of other symmetry points is relevant for the OHE in this system, the dominant contribution to the OHE plateau arrises from the Berry curvature of the valence band at the K points. Q(~k) also presents a secondary broad peak in M . For monolayer TMDs the contribution from M is reduced [2]. Figure 10: Q(~k) along the high symmetry lines of the hexagonal Brillouin zone of the bilayer MoS 2. ORBITAL CHERN NUMBER In the main text, we present an orbital Chern number as the topological invariant which indexes the orbital Hall insulating phase of the TMD bilayers ( CL = 2). To compute this orbital Chern number, we followed the method introduced by D. N. Sheng, et. al. for nite systems [69], and later formalized by E. Prodan [55] in the thermodynamic limit. This method has been applied successfully in dierent situations [56, 57] and it is ideal for systems where the operator associated with the Chern number does not commute with the Hamiltonian, as in the case of the orbital angular momentum ([H̃(~qτ ), Lz ] 6= 0) . 1. First, we calculate the energy spectrum and eigenvectors of the Hamiltonian that describes the electronic struc- ture of the system, such as the one given by Eq. (3) of the main text. Since we are interested in the topological 9 properties of insulators, the band structure must be fully gapped. We may then separate the eigenstates of the Hamiltonian that belong to the conduction- and valence-band subspaces. With our eective model for the s,τ bilayer, the valence-band subspace is spanned by two eigenstates ψn,~q that clearly depend upon s and τ . Here, n = 1, 2 label the two valence energy bands, ~q represents the wavevector relative to valleys, τ = ±1 denote the valley quantum number associated with the K and K 0 symmetry points, and s =↑, ↓ designate the two spin sectors. 2. The second step consists in projecting  the orbital angularmomentum operator in the valence band subspace. To this end, we construct the matrix Lv (~k) = P (~k)Lz P (~k) , where P (~k) is the projection operator on the valence band states. The matrix elements of Lv (~k) are given by ψn,~ s,τ s,τ q Lz ψm,~ q (n, m = 1, 2), which only takes into account the valence band states. We then compute its eigenvalues and eigenvectors. The eigenvalue spectrum of Lv (~q) is presented in Fig. 3 (a) of the main text and it may be separated into positive + q ) and negative s,τ (~ − q ) eigenvalues. The application Prodan's method requires that the eigenvalue spectrum must be gapped. s,τ (~ 3. The third step is the construction of the eigenvectors of Lv (~k) with the use of the coecients αs,τ  ± ±  (~q), βs,τ (~q) obtained with the diagonalization of Lv (~k): s,τ s,τ Φ± ± q ) ψ1,~ s,τ = αs,τ (~ ± q ) ψ2,~ q + βs,τ (~ q . (6) 4. The last step is the computation of the orbital Chern number using Eq. (4) of the main text Z ± 1 X CL = d2 q ± Fs,τ (q), (7) 2π s,τ where Fs,τ ± (q) = −2Im[ ∂qx Φ± q ) ∂qy Φ± s,τ (~ q ) ]. To evaluate the integrals in Eq. (7), we have used a cuto in s,τ (~ −1 the momentum space Λ = 15Å , which provides a numerical deviation smaller than 1% from the quantized value CL = (1/2)(CL+ − CL− ) = 2. Note that this cut-o value is larger than the limit of validity of eective Hamiltonian given by Eq. (3) of the P main text. However, this is not an issue in the context of the low energy continuous theory since the function s,τ Fs,τ ± (q) is strongly peaked in the valleys and decays relatively fast for increasing values of q . Limitations of the analysis The method discussed above relies on a fully gapped spectrum of the orbital angular momentum operator, when projected into the valence band subspace. This condition is fullled for the low-energy Hamiltonian. However, dierently from the case of the spin operator in the quantum spin Hall insulator, the spectrum does not remain opened in the whole Brillouin zone. Fortunately, one can still dene a topological number, in a similar way of what is done for valley Chern numbers [60]. Because the Lz projected Berry curvature is concentrated around the valley points, the integral of eq. 7 is well dened and we can obtain the orbital Chern number. Rigorously, it is an approximation since it does not consider the whole Brillouin zone. In this sense, our analysis has the same limitations of the quantum valley Hall eect. For instance, the existence of edge-states is limited to to edges that preserve the two valleys, as in the case of zigzag nanoribbons. On the other hand, to calculate the OHE, dierently from the VHE, it is possible to use the orbitally projected Berry curvature in the whole BZ, which indicates that it might be possible to nd a dierent strategy to calculate the orbital Chern number using the whole BZ. 035409 (2020), URL https://link.aps.org/doi/10. 1103/PhysRevB.102.035409 . [3] B. A. Bernevig, T. L. Hughes, and S.-C. Zhang, Phys. ∗ Electronic address: tarik.cysne@gmail.com Rev. Lett. 95 , 066601 (2005), URL https://link.aps. [1] L. M. Canonico, T. P. Cysne, A. Molina-Sanchez, org/doi/10.1103/PhysRevLett.95.066601 . R. B. Muniz, and T. G. Rappoport, Phys. Rev. B [4] J. Liu, Z. Ma, J. Gao, and X. Dai, Phys. Rev. X 101 , 161409 (2020), URL https://link.aps.org/doi/ 9, 031021 (2019), URL https://link.aps.org/doi/10. 10.1103/PhysRevB.101.161409 . 1103/PhysRevX.9.031021 . [2] S. Bhowal and S. Satpathy, Phys. Rev. B 102 , [5] D. Go, J.-P. Hanke, P. M. Buhl, F. Freimuth, 10 G. Bihlmayer, H.-W. Lee, Y. Mokrousov, and S. Blügel, Nature Physics 9, 149153 (2013), ISSN 1745-2481, URL Scientic Reports 7 (2017), ISSN 2045-2322, URL http: http://dx.doi.org/10.1038/nphys2524 . //dx.doi.org/10.1038/srep46742 . [25] J. Lee, K. F. Mak, and J. Shan, Nature Nanotechnology [6] T. Yoda, T. Yokoyama, and S. Murakami, 11 , 421425 (2016), ISSN 1748-3395, URL http://dx. Nano Letters 18 , 916 (2018), pMID: 29373028, doi.org/10.1038/nnano.2015.337 . https://doi.org/10.1021/acs.nanolett.7b04300, URL [26] Z. Gong, G.-B. Liu, H. Yu, D. Xiao, X. Cui, X. Xu, and https://doi.org/10.1021/acs.nanolett.7b04300 . W. Yao, Nature Communications 4 (2013), ISSN 2041- [7] S. Beaulieu, J. Schusser, S. Dong, M. Schüler, T. Pincelli, 1723, URL http://dx.doi.org/10.1038/ncomms3053 . M. Dendzik, J. Maklar, A. Neef, H. Ebert, K. Hricovini, [27] A. Kormányos, V. Zólyomi, V. I. Fal'ko, and G. Burkard, et al., Physical Review Letters 125 (2020), URL https: Phys. Rev. B 98 , 035408 (2018), URL https://link. //doi.org/10.1103/physrevlett.125.216404 . aps.org/doi/10.1103/PhysRevB.98.035408 . [8] D. Jo, D. Go, and H.-W. Lee, Phys. Rev. B 98 , [28] A. Kormányos, G. Burkard, M. Gmitra, J. Fabian, V. Zó- 214405 (2018), URL https://link.aps.org/doi/10. lyomi, N. D. Drummond, and V. Fal'ko, 2D Mater. 2 1103/PhysRevB.98.214405 . (2015), URL https://dx.doi.org/10.1088/2053-1583/ [9] J.-H. Park, C. H. Kim, H.-W. Lee, and J. H. Han, Phys. 2/2/022001 . Rev. B 87 , 041301 (2013), URL https://link.aps.org/ [29] L. Du, M. Liao, G.-B. Liu, Q. Wang, R. Yang, doi/10.1103/PhysRevB.87.041301 . D. Shi, Y. Yao, and G. Zhang, Phys. Rev. B 99 , [10] D. Go, D. Jo, C. Kim, and H.-W. Lee, Phys. Rev. Lett. 195415 (2019), URL https://link.aps.org/doi/10. 121 , 086602 (2018), URL https://link.aps.org/doi/ 1103/PhysRevB.99.195415 . 10.1103/PhysRevLett.121.086602 . [30] P. Hohenberg and W. Kohn, Phys. Rev. 136 , [11] F. Xue, C. Rohmann, J. Li, V. Amin, and P. Haney, B864 (1964), URL https://link.aps.org/doi/10. Phys. Rev. B 102 , 014401 (2020), URL https://link. 1103/PhysRev.136.B864 . aps.org/doi/10.1103/PhysRevB.102.014401 . [31] W. Kohn and L. J. Sham, Phys. Rev. 140 , [12] F. Castro de Lima, G. J. Ferreira, and R. H. Miwa, A1133 (1965), URL https://link.aps.org/doi/10. Nano Lett. 19 , 6564 (2019), URL https://pubs.acs. 1103/PhysRev.140.A1133 . org/doi/abs/10.1021/acs.nanolett.9b02802 . [32] P. Giannozzi, O. Andreussi, T. Brumme, O. Bunau, [13] L. M. Canonico, T. P. Cysne, T. G. Rappoport, and R. B. M. Buongiorno Nardelli, M. Calandra, R. Car, Muniz, Phys. Rev. B 101 , 075429 (2020), URL https: C. Cavazzoni, D. Ceresoli, M. Cococcioni, et al., //link.aps.org/doi/10.1103/PhysRevB.101.075429 . Journal of Physics: Condensed Matter 29 , 465901 [14] S. Park and B.-J. Yang, Nano Letters 20 , 7694 (2020), (2017), URL http://stacks.iop.org/0953-8984/29/i= ISSN 1530-6984, URL https://doi.org/10.1021/acs. 46/a=465901 . nanolett.0c03220 . [33] J. P. Perdew, K. Burke, and M. Ernzerhof, Phys. Rev. [15] J. Xiao, Y. Liu, and B. Yan, arXiv preprint Lett. 77 , 3865 (1996), URL https://link.aps.org/ arXiv:2010.01970 (2020). doi/10.1103/PhysRevLett.77.3865 . [16] D. Go, F. Freimuth, J.-P. Hanke, F. Xue, O. Gomonay, [34] G. Kresse and D. Joubert, Phys. Rev. B 59 , 1758 (1999), K.-J. Lee, S. Blügel, P. M. Haney, H.-W. Lee, and URL https://link.aps.org/doi/10.1103/PhysRevB. Y. Mokrousov, Phys. Rev. Research 2, 033401 59.1758 . (2020), URL https://link.aps.org/doi/10.1103/ [35] S. Grimme, Journal of Computa- PhysRevResearch.2.033401 . tional Chemistry 27 , 1787 (2006), [17] D. Go and H.-W. Lee, Phys. Rev. Research 2, https://onlinelibrary.wiley.com/doi/pdf/10.1002/jcc.20495, 013177 (2020), URL https://link.aps.org/doi/10. URL https://onlinelibrary.wiley.com/doi/abs/10. 1103/PhysRevResearch.2.013177 . 1002/jcc.20495 . [18] Z. C. Zheng, Q. X. Guo, D. Jo, D. Go, L. H. [36] L. A. Agapito, A. Ferretti, A. Calzolari, S. Cur- Wang, H. C. Chen, W. Yin, X. M. Wang, G. H. tarolo, and M. Buongiorno Nardelli, Phys. Rev. B 88 , Yu, W. He, et al., Phys. Rev. Research 2, 013127 165127 (2013), URL https://link.aps.org/doi/10. (2020), URL https://link.aps.org/doi/10.1103/ 1103/PhysRevB.88.165127 . PhysRevResearch.2.013127 . [37] L. A. Agapito, M. Fornari, D. Ceresoli, A. Ferretti, [19] J. Kim, D. Go, H. Tsai, D. Jo, K. Kondou, H.-W. Lee, S. Curtarolo, and M. Buongiorno Nardelli, Phys. Rev. and Y. Otani, arXiv preprint arXiv:2002.00596 (2020). B 93 , 125137 (2016), URL https://link.aps.org/doi/ [20] Y. Chen, W. Ruan, M. Wu, S. Tang, H. Ryu, H.-Z. Tsai, 10.1103/PhysRevB.93.125137 . R. Lee, S. Kahn, F. Liou, C. Jia, et al., Nature Physics [38] H. Kontani, T. Tanaka, D. S. Hirashima, K. Yamada, 16 , 218224 (2020), ISSN 1745-2481, URL http://dx. and J. Inoue, Phys. Rev. Lett. 100 , 096601 (2008), URL doi.org/10.1038/s41567-019-0744-9 . https://link.aps.org/doi/10.1103/PhysRevLett. [21] F. Xue, V. Amin, and P. M. Haney, Phys. Rev. B 100.096601 . 102 , 161103 (2020), URL https://link.aps.org/doi/ [39] I. V. Tokatly, Phys. Rev. B 82 , 161404 (2010), 10.1103/PhysRevB.102.161103 . URL https://link.aps.org/doi/10.1103/PhysRevB. [22] K. F. Mak, D. Xiao, and J. Shan, Nature Photon 82.161404 . 12 , 451460 (2018), URL https://www.nature.com/ [40] V. o. T. Phong, Z. Addison, S. Ahn, H. Min, articles/s41566-018-0204-6 . R. Agarwal, and E. J. Mele, Phys. Rev. Lett. 123 , [23] H. Zeng, J. Dai, W. Yao, D. Xiao, and X. Cui, Na- 236403 (2019), URL https://link.aps.org/doi/10. ture Nanotechnology 7, 490493 (2012), ISSN 1748-3395, 1103/PhysRevLett.123.236403 . URL http://dx.doi.org/10.1038/nnano.2012.95 . [41] M. Buongiorno Nardelli, F. T. Cerasoli, M. Costa, [24] S. Wu, J. S. Ross, G.-B. Liu, G. Aivazian, A. Jones, S. Curtarolo, R. D. Gennaro, M. Fornari, L. Liyan- Z. Fei, W. Zhu, D. Xiao, W. Yao, D. Cobden, et al., age, A. R. Supka, and H. Wang, Computational 11 Materials Science 143 , 462 (2018), ISSN 0927- Xing, and D. N. Sheng, Phys. Rev. Lett. 107 , 0256, URL http://www.sciencedirect.com/science/ 066602 (2011), URL https://link.aps.org/doi/10. article/pii/S0927025617306651 . 1103/PhysRevLett.107.066602 . [42] M. Costa, G. R. Schleder, M. Buongiorno Nardelli, [57] H. Li, L. Sheng, D. N. Sheng, and D. Y. Xing, Phys. Rev. C. Lewenkopf, and A. Fazzio, Nano Letters 19 , 8941 B 82 , 165104 (2010), URL https://link.aps.org/doi/ (2019), URL https://doi.org/10.1021/acs.nanolett. 10.1103/PhysRevB.82.165104 . 9b03881 . [58] H. Rostami, R. Asgari, and F. Guinea, Jour- [43] M. Costa, G. R. Schleder, C. M. Acosta, A. Padilha, nal of Physics: Condensed Matter 28 , 495001 F. Cerasoli, M. B. Nardelli, and A. Fazzio, arXiv preprint (2016), URL https://doi.org/10.1088%2F0953-8984% arXiv:2006.07270 (2020). 2F28%2F49%2F495001 . [44] M. Costa, N. M. R. Peres, J. Fernández-Rossier, [59] D. Go, D. Jo, T. Gao, K. Ando, S. Blügel, H.-W. Lee, and and A. T. Costa, Phys. Rev. B 102 , 014450 (2020), Y. Mokrousov, arXiv preprint arXiv:2011.08601 (2020). URL https://link.aps.org/doi/10.1103/PhysRevB. [60] F. Zhang, A. H. MacDonald, and E. J. Mele, Proceedings 102.014450 . of the National Academy of Sciences 110 , 10546 (2013), [45] A. M. García, E. d. Corro, M. Kalbac, and O. Frank, URL https://doi.org/10.1073/pnas.1308853110 . Phys. Chem. Chem. Phys. 19 , 13333 (2017), URL http: [61] Z. Wu, B. T. Zhou, X. Cai, P. Cheung, G.-B. Liu, //dx.doi.org/10.1039/C7CP00012J . M. Huang, J. Lin, T. Han, L. An, Y. Wang, et al., Nature [46] See Supplemental Material, which includes Refs. [6369] Communications 10 (2019), URL https://doi.org/10. for additional analyses, and details on the numerical 1038/s41467-019-08629-9 . simulations. [62] T. Brumme, M. Calandra, and F. Mauri, Phys. Rev. B 91 , 155436 (2015), URL https://link.aps.org/doi/ [47] D. Xiao, G.-B. Liu, W. Feng, X. Xu, and W. Yao, Phys. 10.1103/PhysRevB.91.155436 . Rev. Lett. 108 , 196802 (2012), URL https://link.aps. [63] D. Moldovan, Mi²a Anželkovi¢, and Francois Peeters, org/doi/10.1103/PhysRevLett.108.196802 . pybinding v0.9.5: a python package for tight-binding [48] G.-B. Liu, W.-Y. Shan, Y. Yao, W. Yao, and D. Xiao, calculations (2020), URL https://zenodo.org/record/ Phys. Rev. B 88 , 085433 (2013), URL https://link. 4010216 . aps.org/doi/10.1103/PhysRevB.88.085433 . [64] S. M. João, M. Anželkovi¢, L. Covaci, T. G. Rappoport, [49] N. Zibouche, P. Philipsen, A. Kuc, and T. Heine, Phys. J. M. V. P. Lopes, and A. Ferreira, Royal Society Open Rev. B 90 , 125440 (2014), URL https://link.aps.org/ Science 7, 191809 (2020), URL https://doi.org/10. doi/10.1103/PhysRevB.90.125440 . 1098/rsos.191809 . [50] J.-S. You, S. Fang, S.-Y. Xu, E. Kaxiras, and T. Low, [65] E. Ridol, L. R. F. Lima, E. R. Mucciolo, and C. H. Phys. Rev. B 98 , 121109 (2018), URL https://link. Lewenkopf, Phys. Rev. B 95 , 035430 (2017), URL https: aps.org/doi/10.1103/PhysRevB.98.121109 . //link.aps.org/doi/10.1103/PhysRevB.95.035430 . [51] D. Xiao, M.-C. Chang, and Q. Niu, Rev. Mod. Phys. [66] I. Souza, J. Íñiguez, and D. Vanderbilt, Phys. Rev. Lett. 82 , 1959 (2010), URL https://link.aps.org/doi/10. 89 , 117602 (2002), URL https://link.aps.org/doi/ 1103/RevModPhys.82.1959 . 10.1103/PhysRevLett.89.117602 . [52] X. Xu, W. Yao, D. Xiao, and T. F. Heinz, Nature Physics [67] P. Umari and A. Pasquarello, Phys. Rev. Lett. 89 , 10 , 343 (2014). 157602 (2002), URL https://link.aps.org/doi/10. [53] C. L. Kane and E. J. Mele, Phys. Rev. Lett. 95 , 1103/PhysRevLett.89.157602 . 226801 (2005), URL https://link.aps.org/doi/10. [68] A. Ramasubramaniam, D. Naveh, and E. Towe, Phys. 1103/PhysRevLett.95.226801 . Rev. B 84 , 205325 (2011), URL https://link.aps.org/ [54] M. Ezawa, Phys. Rev. B 88 , 161406 (2013), URL https: doi/10.1103/PhysRevB.84.205325 . //link.aps.org/doi/10.1103/PhysRevB.88.161406 . [69] D. N. Sheng, Z. Y. Weng, L. Sheng, and F. D. M. Hal- [55] E. Prodan, Phys. Rev. B 80 , 125327 (2009), URL https: dane, Phys. Rev. Lett. 97 , 036808 (2006), URL https:// //link.aps.org/doi/10.1103/PhysRevB.80.125327 . link.aps.org/doi/10.1103/PhysRevLett.97.036808 . [56] Y. Yang, Z. Xu, L. Sheng, B. Wang, D. Y.