Band energy landscapes in twisted homobilayers of transition metal dichalcogenides

Twistronic assembly of 2D materials employs the twist angle between adjacent layers as a tuning parameter for designing the electronic and optical properties of van der Waals heterostructures. Here, we study how interlayer hybridization, weak ferroelectric charge transfer between layers, and piezoelectric response to deformations set the valence and conduction band edges across the moir{\'e} supercell in twistronic homobilayers of MoS$_2$, MoSe$_2$, WS$_2$ and WSe$_2$. We show that, due to the lack of inversion symmetry in the monolayer crystals, bilayers with parallel (P) and anti-parallel (AP) unit cell orientations display contrasting behaviors. For P-bilayers at small twist angles we find band edges in the middle of triangular domains of preferential stacking. In AP-bilayers at marginal twist angles ($\theta_{AP}<1^\circ$) the band edges are located in small regions around the intersections of domain walls, giving highly localized quantum dot states.

The assembly of van der Waals heterostructures has potential for tailoring the properties of 2D materials 1 . Recently, it has been shown that twisted bilayer graphene exhibits intrinsic unconventional superconductivity 2 and a Mott insulating phase 3 . These effects have been related to the localization of electrons in particular stacking areas of the moiré superlattice (mSL), potentially enhanced by lattice reconstruction, promoting energetically favorable Bernal stacking [4][5][6] .
Lattice reconstruction also takes place in bilayers of twisted transition metal dichalcogenides [7][8][9][10][11][12][13][14][15][16] (TMDs), where it gives rise to the formation of preferential stacking domains separated by a network of domain walls (DWs). The shape and size of these domains depend on the mutual orientation of the unit cells of the individual crystals, which can be parallel (P) or anti-parallel (AP), whereas the superlattice period depends on the misalignment angle between the crystallographic axes of the layers. These two orientations fundamentally differ in that AP-bilayers possess inversion symmetry for all local stacking configurations, whereas for P-bilayers both inversion and mirror reflection symmetries are generally broken. Below we discuss how this symmetry breaking in P-bilayers leads to: (a) interlayer charge transfer due to hybridization of the conduction bands in one layer with the valence bands of the other; (b) layer-asymmetric piezoelectric charges caused by lattice reconstruction, concentrated around the network of domain walls and corners. Some of these effects were recently a) Electronic mail: fabio.ferreira@postgrad.manchester.ac.uk b) Electronic mail: vladimir.falko@manchester.ac.uk discussed in relation to the structural characterization of TMD homobilayers and heterobilayers [7][8][9][10][11][12][13][14][15][16] , and in twisted hexagonal boron nitride (hBN) [17][18][19] . In both P-and AP-bilayers, the band energies depend on the interplay between the above effects and the interlayer hybridization of band edge states, which is most prominent for the states around the Γ-and Qvalleys and marginally relevant for the K-valley band edges.
To model twisted TMD homobilayers, we use a multiscale approach 10 for describing lattice reconstruction and computing piezoelectric charge and potential distributions, complemented by DFT-parametrized Hamiltonians for interlayer hybridization of band edge states developed to describe an arbitrary local stacking. The first step is to use an earlier parametrized 10 interlayer adhesion potential to compute atomic reconstruction of twisted TMDs for various twistangles using elasticity theory. We take into account both lateral and vertical deformations of the layers, see S1 in Supplementary Material (SM). This gives us the pattern of local interlayer distance d and lateral offset r 0 , which determine the stacking of the layers, as well as the piezoelectric potential distribution across the mSL. The second step is to compute the band structures of bilayers with various offsets (using DFT implemented in Quantum ESPRESSO package 20 ) in order to parametrize r 0 -dependent Hamiltonians describing interlayer hybridization (see S2 and S3 in SM). Diagonalizing the Hamiltonians, which take into account piezopotentials, we build maps for the local band edge energies across the mSL.
We analyze the band edge properties in the valence band (VB), near both the Γ-and K-valleys, and the conduction band (CB), near both the K-and Q-valleys. We consider both op-tions (Γ and K) for the location of the VB edge in the Brillouin zone, as their relative energies may depend on the encapsulating environment. In the presence of hBN, the energies of chalcogen orbitals, which are strongly represented at the Γvalley, may be shifted with respect to the metal d x 2 −y 2 and d xy orbitals that determine the K-valley energy 21 , with the magnitude and direction of the shift determined by the WSe 2 /hBN band alignment, and the properties of their hybridization.
Local stacking configurations, distinguished by relative lateral offset of the two lattices in the bilayer (see Fig. S2 in SM), vary across the moiré supercell, with lattice reconstruction promoting energetically favorable configurations. For Pstacking, these are MX /XM , in which a metal atom M(M ) is below(above) a chalcogen atom X(X ) in the opposite layer (X and X refer to chalcogen atoms in bottom and top layer; M and M refer to transition metal atoms in bottom and top layer). For AP-stacking, the favorable configuration is 2H, as found in bulk TMD crystals, and features two interlayer pairs of vertically aligned metal/chalcogen atoms. Other important high-symmetry stackings include XX (both P and AP), for which chalcogen atoms are vertically aligned, and for AP stacking MM , which has aligned metal atoms (Fig. S2 in SM). P-bilayer superlattices are comprised of two triangular domains with MX and XM stacking separated by domain walls with XX stacking near hexagonal domain wall lattice sites, whereas in AP-bilayers there are large 2H domains separated by a honeycomb network of domain walls with XX and MM stacking at their intersections 10 . Below, we aim to determine which stacking areas in the moiré supercell host the band edges, which is essential for determining and predicting the formation of quantum dots in marginally twisted structures. We discuss separately features of the K-, Q-and Γ-valley band edges for P-and AP-bilayers.
FIG. 1. Left panels: Variation of the VB energy with twist angle θ AP for different stacking configurations in AP-bilayers at the Γ-and K-valleys. The VB edge at the Γ-valley is located at the corners of 2H domains (labeled as 2H c ) for marginal twist angles. These corners can be seen in the Γ-valley VB edge map for θ AP = 0.2 • . The arrows in the left panels mark the difference between the Γ-valley VB energy at the corners and in the center of 2H domains. Right panels: Variation of the CB energy with θ AP for different stacking configurations in AP-bilayers at the K-and Q-valleys. The middle panels show maps for the VB edge at the Γ-and K-valleys, and for the CB edge at the Kand ±Q 1 -valleys in AP-MoS 2 for θ AP = 0.2 • and θ AP = 3 • . The zigzag orientations of DWs make the Q-valley CB edge maps C 3 -symmetric, leading to the same CB edge landscapes for ±Q 1,2,3 . The vacuum level is set to 0 eV. See Figs. S4-S7 with maps for the other TMDs.
For AP-bilayers the position of the VB at Γ is determined by the interplay between the piezopotential and interlayer hybridization of resonant band edges. Due to symmetry, the piezopotential induced by lateral lattice reconstruction is the same on both layers 10 ; interlayer hybridization is also sensitive to lattice deformations, but in the vertical direction. This is because the wavefunctions at Γ carry a substantial weight of p z orbitals of chalcogens, which strongly overlap between the layers with a pronounced dependence on the interlayer distance 21 .
2H domains, which provide energetically more favorable stacking, also feature the closest interlayer distance, promot-FIG. 2. Left panels: Variation of the VB energy with twist angle θ P for XX -and MX -stacking configurations in P-bilayers at the Γ-and K-valleys. Right panels: Variation of the CB energy with θ P for XX -and MX -stacking configurations in P-bilayers at the K-and Q-valleys. The middle panels show maps for the VB edge at the Γ-and K-valleys, and for the CB edge at the K-and ±Q 1 -valleys in P-MoS 2 for θ P = 0.2 • and θ P = 3 • . Note that for all angles the Q-valley maps are anisotropic, with the anisotropy axis rotated by ±120 • when going from Q 1 to Q 2 and Q 3 . For θ P = 3 • , the low contrast of the CB and VB K-valley maps reflect a negligible variation of the band edge energy (< 20 meV). The vacuum level is set to 0 eV. See Figs. S8-S11 for maps for the other TMDs.
ing the higher position of the top VB due to interlayer hybridization of the VB edges. At the same time, piezocharges are largest and negative (hence, the electron piezopotential is highest) at XX regions 10 , attracting holes from the centers of 2H domains toward XX corners. The trigonal symmetry of the honeycomb domain structure suggests that the resulting VB edges appear as three maxima labeled in Fig. 1 as 2H c . This behavior is corroborated by the plots shown on the left hand side of Fig. 1, where band edge energies in areas with local 2H, MM and XX stacking configurations are compared with each other, and with the energy of those local maxima at 2H c . The data shown in Fig. 1 for all TMDs indicate that they all display the same VB behavior at Γ.
For large twist angles, the piezopotential magnitude decreases, whereas the contribution of interlayer hybridization remains the same. This promotes the maximum value for the VB energy at the center of 2H domains, while the minimum value switches from MM -to XX -stacking regions. This is shown in the left-hand-side panels of Fig. 1, where the crossover between the two regimes takes place at θ AP ≈ 0.6 • .
For AP-bilayers the VB and CB edge energies at K are dominated by the piezopotential contribution. This is because interlayer hybridization between K-valley states in APbilayers is weak: the band edge states are dominated by the orbitals of metals, buried inside the layer 21 , and their spinvalley locking due to spin-orbit splitting makes monolayer bands with the same spin off-resonant.
A domain wall structure is fully developed for twist angles θ AP 0.6 • , giving rise to large piezopotentials around domain wall intersections. For such twist angles, the piezopotential forms quantum-dot-like wells of depth ∼150 meV for electrons and holes in the MM and XX intersections of the domain wall structure, respectively, shown in Fig. 1. In a periodic moiré supercell, quantized states localized in each of these quantum dots should give rise to narrow bands located at the CB and VB edges of marginally twisted bilayers. Upon increasing the twist angle, the magnitude of the piezopotential decreases and the location of its minimum in the supercell shifts to 2H areas, followed by the position of the VB band edge. The CB edge remains in the MM -stacked regions due to a residual attractive piezopotential for electrons even at larger angles.
For AP-bilayers the position of the CB edge at the six Q-valleys is determined by the piezopotential and by interlayer coupling variations across the mSL. For θ AP ≤ 1 • , the piezopotential overcomes hybridization, forming quantum dots for electrons at MM corners of 2H domains. For θ AP > 3 • , the Q-valley CB edge shifts to the center of 2H re-gions.
For P-bilayers the VB edge at Γ is dominated by interlayer hybridization of chalcogen orbitals, which is an order of magnitude stronger than both the piezoelectric and ferroelectric potentials (Table S4 in SM). A weak ferroelectric potential arises due to the interlayer charge transfer caused by off-resonant hybridization of conduction bands in one layer with valence bands in the opposite layer. Such a spontaneous out-of-plane charge polarisation is allowed due to the lack of inversion symmetry in P-bilayers 22,23 . For the full range of twist angles, the VB at Γ lies in the MX /XM domains, where hybridization is largest because the interlayer distance is smallest (see Fig. 2).
For P-bilayers the VB and CB edges at K are determined by a competition between interlayer charge transfer, the piezopotential, and resonant interlayer hybridization at the band edges. All three are comparable for the K-valley states, but vary differently across the mSL. Interlayer hybridization vanishes for MX and XM stackings due to the symmetry of the K-valley Bloch states in the honeycomb lattice 24 , whereas the potential step due to ferroelectric charge transfer is largest (Tables S4 and S5 in SM). The piezopotential, which has opposite signs in the two layers, contributes to the splitting of the resonantly coupled bands rather than to the modulation of their average position, as in AP-bilayers 25 . Splitting in XX areas is only determined by interlayer hybridization, as the ferroelectric potential and piezopotential are absent there. The piezopotential distribution strongly depends on twist angle. Also, piezo-and ferroelectric potentials have opposite signs, and fully or partially cancel each other out in areas where both are present. At marginal twist angles θ P 1 • , this leads to flat CB and VB edges across MX and XM domains, see Fig. 2, forming shallow (≈ 30 meV) triangular traps for electrons and holes, respectively.
As the twist angle is increased, the piezopotential expands into the MX /XM domains, with similar magnitude and opposite sign to the potential caused by the charge transfer effect, reducing the splitting between the bands (SM Fig. S13). For the VB, this is sufficient to move the band edge to the XX areas, where the splitting due to weak interlayer hybridization persists, and for the CB the MX /XM and XX regions become very close in energy, see right hand panels of Fig. 2.
For P-bilayers, the Q-valley CB edge is mainly affected by the variation of resonant interlayer hybridization of monolayer Q-valley states across the mSL, with only weak effects from the layer-asymmetric piezo-and ferroelectric potentials. For θ P ≤ 2 • , the CB edge appears at one-dimensional channels zig-zagging across the mSL. Such channels are ∼ 25 meV deep for P-WSe 2 , but shallower for the other TMDS. The asymmetry seen in the band edge maps is due to the underlying asymmetry of the Q-valley wave functions in each layer, interplaying with the lateral interlayer offset, which is differently oriented at different segments of the DW web. In this regime, XX areas play the role of potential barriers, acting as scatterers between the channels for CB electrons. For θ P ≥ 2 • , the channels broaden to form anisotropic landscapes. Despite their small amplitude, "potential" variations in such landscapes would determine anisotropic moiré minibands for Q-valley electrons.
Finally, we have established how the interplay between interlayer hybridization, piezoelectric potential and ferroelectric charge transfer affects the band edges for electrons in the Kand Q-valleys in the CB, and K-and Γ-valleys in the VB of TMD bilayers. We have analyzed all these possibilities because the influence of the encapsulating material or a substrate may affect their mutual alignment. For example, the DFT modelling of isolated bilayers points towards Q-valley CB edge for all four TMDs studied here, whereas recent experimental studies 26 of 2H MoS 2 bi-and tri-layers indicate the persistence of K-valley band edges. Then, the crossovers between different misalignment angle regimes presented in this paper and highlighted in Figs. 1 and 2 represent the variety of scenarios for the formation of quantum dot arrays that determines Hubbard physics for narrow minibands of electrons and holes in twistronic bilayers with different encapsulation environments.

SUPPLEMENTARY MATERIAL
See the Supplementary Material for more details about the multiscale approach for describing lattice reconstruction and computing piezoelectric charge and potential distributions. More details about DFT bandstructure calculations and effective Hamiltonians used to describe band edges at Γ-, K and Q-valleys can also be seen in the Supplementary Material. Supplementary material for: Band energy landscapes in twisted homobilayers of transition metal dichalcogenides

S1. LATTICE RECONSTRUCTION AND PIEZOPOTENTIAL
Modeling of lattice relaxation in twisted TMD bilayers. We describe lattice reconstruction in twisted bilayers in frames of mesoscale elasticity theory supplemented by DFT-parametrized adhesion energy derived in Refs. 10,11,15 . In this approach we introduce in-plane deformation fields u t and u b , which ensure local minimum for sum of elastic and adhesion energies in the moiré supercell: Here, λ t/b , µ t/b are elastic moduli of TMD monolayers (listed in Table S1), u is top/bottom layer in-plane strain tensor i, j = x, y, and W P/AP (r 0 (r)) is the adhesion energy taken at local offset between the twisted layers 10,15 Reference of the local offsets r 0 (0) = (0, 0) in Eq. (S2) corresponds to XX configuration for P/AP-orientation of bilayers shown in Fig. S2. Explicit form for the adhesion energy reads where γ AP = 0, γ P = π/2, G 1,2,3 are the shortest reciprocal vectors of TMD monolayer with magnitude G related by C 3 -rotation (see Fig. S1), the term −κZ 2 (r 0 ) takes into account relaxation of interlayer distances d 0 + Z(r 0 ) with stacking configurations as it was derived in Ref. 10 (d 0 is optimal distance for stacking-averaged adhesion energy). Values of parameters in adhesion energy (S3) for studied homobilayers are listed in Table S1.  Minimization of the functional (S1) results in a system of four Euler-Lagrange equations that are discretized on a rectangular grid using the finite-difference method. The package GEKKO Optimization Suite 31 is used to solve this system of non-linear equations, in which the final solution produces a detailed structure of the strain fields.
Piezopotential. Lacking inversion center, TMD crystals become piezoactive in the limit of a single monolayer 29 . Therefore, deformations in each layer induced by lattice relaxation result in the emergence of piezoelectric charge densities in the top/bottom layers: where e t/b 11 are piezocoefficients listed in Table S1. The latter are the same (e t 11 = e b 11 > 0) for P-and opposite (e b 11 = −e t 11 > 0) for AP-bilayers. We note that opposite signs of the deformation fields in the two layers u t = −u b lead to the same sign of the piezocharges in the constituent layers for AP-bilayers and opposite signs of piezocharges for P-bilayers. To calculate the potential φ t/b created by the piezocharges in the top/bottom layers we expand them in Fourier series over the moiré superlattice reciprocal vectors solving Poisson equation with appropriate boundary conditions on interfaces taking into account in-plane polarisation of orbitals and hBN encapsulation as it is described in Ref. 15 . Lattice reconstruction results from the competition between elastic and adhesion forces. This is expressed in the following functional 10 .
where r 0 is the local shift between chalcogen atoms in each layer, d is the interlayer distance, U is the elastic energy and W AP/P is the local adhesion energy that takes into account in-plane and out-plane reconstruction. The expression for the elastic energy is parametrised as follows where λ t/b is first Lamé coefficients for the top (t) and bottom (b) layer, µ t/b is the shear modulus and u is the strain tensor. The formula for W AP/P is expressed as follows +A 2 e −Gd 0 sin g n r 0 + G n u (t) − u (b) + ϕ P/AP , where g n are the reciprocal vectors of the moiré supercell, G n are the reciprocal vectors of TMD where G is their magnitude, 2 is stacking-averaged term and ϕ P = π/2 and ϕ AP = 0, The coefficients A 1,2 , ρ, ε, C 1,2,3 are obtained by fitting the adhesion energy to ab initio DFT calculations as explained in more detail in Ref. 10. The reconstructed local shift between chalcogen atoms is given by where δ is the lattice mismatch between two layers and θẑ is a rotation that results from the misalignment angle. The displacement vectors can be obtained by minimizing Eq. (S5) with respect to u t/b . This minimization results in a system of Euler-Lagrange equations and it is solved using finite difference method with the help of the package GEKKO Optimization Suite 32 . Lattice reconstruction will also influence the rearrangement of the charges in this TMDs. Thus, piezoelectricity plays an important role on these materials given that they lack inversion symmetry. The expression for bare piezocharges is given by 10 where the piezocoefficient e t 11 = e b 11 for P-bilayers and e b 11 = −e t 11 > 0 for AP-bilayers. Computation of piezopotential φ t/b is determined by taking into account screening-induced charges as well. The procedure for such a computation is the one used in Ref. 10 , where we considered TMDs to be encapsulated by hexagonal boron nitride.

S2. DFT COMPUTATIONS
Ab initio density-functional-theory calculations were carried out using the Quantum ESPRESSO package 20 . We used full relativistic pseudopotentials with spin-orbit interaction included. The exchange-correlation functional used was the generalized gradient approximation of Perdew-Burke-Ernzerhof (GGA-PBE) 33 . A plane-wave cutoff energy of 1090 eV was used for all calculations. The integration over the Brillouin-zone (BZ) was performed using the scheme proposed by Monkhorst-Pack with a grid of 12 × 12 × 1 k-points 34 . Monolayer structure parameters of bulk TMDs taken from Refs. 35 and 36 were used for TMD homobilayers. The interlayer hybridization models set out below were parametrized from DFT bands calculated using a range of in-plane relative shifts of the layers, r 0 , and interlayer distances, d (see Fig. S2 and Tables S2 and S3).

S3. INTERLAYER HYBRIDIZATION MODELS
To determine modulation of valence band edges at Γ, K-valleys and conduction band edge at K-and Q-valleys in twisted TMD homobilayers, we, first, establish effective Hamiltonians describing the interlayer coupling of the corresponding states in aligned AP-and P-bilayers for a given lateral offset r 0 and interlayer distance d between the layers (see Fig. S2). In these effective Hamiltonians we use translational symmetry of the aligned bilayers and represent all matrix elements in terms of Fourier series over reciprocal vectors of a TMD monolayer, keeping only the lowest harmonics and accounting for symmetry of restrictions.
For Γ-valley, an effective Hamiltonian can be represented as: H P/AP In-plane relative shifts of the layers used to parametrize the interlayer hybridization models from DFT bands. In-plane relative shifts from 1 to 12 were used to extract energies from DFT bandstructure calculations at Γ-and K-valleys for AP-bilayers. In-plane relative shifts from 1 to 7 and 13 to 14 were used to extract energies from DFT bandstructure calculations at both Γ-and K-valleys for P-bilayers. In-plane relative shifts 1, 10, 15 to 24 were used to extract energies from DFT bandstructure calculations at Q-valley for AP-bilayers. In-plane relative shifts 1, 7, 10, 15, 16 and 18 were used to extract energies from DFT bandstructure calculations at Q-valley for P-bilayer. and Hamiltonians (S10-S12) act on a two-component wave function, where the first component describes the top layer state, the second -states of the bottom layer. The first term (S10), responsible for the interlayer hybridization of the monolayer states, is the same for P-and AP-bilayers. The second one (S11), describing the electron potential energy shift, ∆ P , induced by the interlayer charge transfer, is applied only for P-bilayers. The last term (S12) takes into account r 0 → −r 0 asymmetry of resonant Interval of interlayer distances that were used to compute DFT bandstructure calculations using in-plane relative shifts listed in Table S2. # d Å MoS 2 MoSe 2 WS 2 WSe 2 1 6.0 Γ, K, Q -Γ, K, Q -2 6.1 Γ, K -Γ, K -3 6.1489 Γ, K, Q ---4 6.1725 --Γ, K, Q -5 6.2 Γ, K -Γ, K -6 6.3 Γ, K, Q Γ, K Γ, K, Q Γ, K 7 6.4 Γ, K Γ, K, Q Γ, K Γ, K, Q 8 6.463 -Γ, K, Q -9 6.477 ---Γ, K, Q 10 6.5 Q Γ, K Q Γ, K 11 6.6 Γ, K Γ, K, Q Γ, K Γ, K, Q 12 6.7 Γ, K, Q Γ, K Γ, K, Q Γ, K 13 6.8 Γ, K Γ, K, Q Γ, K Γ, K, Q 14 6.9 - coupling in AP-bilayers. Matrix elements in Eqs. (S10-S12) read as follows: sin(G j r 0 ). (S13) The interlayer distance dependence of the matrix elements (S13) were extracted from fitting of the Γ-valley valence band edge in P-and AP-TMD bilayers calculated in DFT for various offsets and interlayer distances listed in Tables S2 and S3. Interlayer distances are chosen by taking the r 0 dependence found using the methods set out in Ref. 10, then rigidly shifting the resulting values of d such that the value for 2H stacking agrees with the experimentally determined quantity for bulk crystals. Results of the fitting are gathered in Table S4.  Effective Hamiltonians describing coupling of the valence band at the K-valley for P-and AP-bilayers are given by where the matrix elements read as follows: sin (G j · r 0 ). (S16) The quantity ∆ SO VB is the spin-orbit splitting in the correspond band, τ and s = −τ is the valley and spin index respectively. The other terms were already defined above. The effective hybridization Hamiltonian for conduction band at the K-valley for P and AP are given by where the matrix elements read (S19) Fitting parameters for quantities on the left hand side of Eqs. (S16-S19) using the expression Ae −q(d−d 0 ) are shown in Table S5.
For Q 1 -valley (see Fig. S1), the effective Hamiltonian is expressed as: where the matrix elements are given by where ∆ Q SO is the monolayer spin-orbit splitting at Q. In Table S6 are listed the fitting parameters for the left hand side quantities in Eqs. S21. Matrix elements for Q 2 and Q 3 can be obtained by applying the C 1 3 and C −1 3 rotation operations, respectively, to G 1,2,3 in phase factors of terms in Eqs. S21. Energies of −Q 1,2,3 correspond to those of Q 1,2,3 respectively.
To apply the effective Hamiltonians in Eqs. S10, S14-S15, S17-S18 and S20, characterising coupling in aligned AP-and Pbilayers, to twisted structures, we substitute r 0 (r) in matrix elements by local r 0 (r) given by Eq. S2 and local interlayer distance d by d(r) = d 0 + Z(r 0 (r)). Finally, we add potential energy of electron in piezopotential (diag(−eϕ, −eϕ) for AP-bilayers and diag(−eϕ, eϕ) for P-bilayers) to the effective Hamiltonians, where ϕ ≡ ϕ t . Then, we diagonalize resulting Hamiltonians and obtain the following expressions for band edge variation in twisted TMD bilayers: In Figs. S4-S11 are shown band edge maps of homobilayers considered in this work for twist-angles of 0.2 • and 3 • . Maps of piezopotentials across the moiré supercell for AP-and P-MoS 2 are shown in Figs. S12-S13. In Table S7 are summarized the dominant effects that determine the location of the band edge in the valence and conduction bands for TMD homobilayers under consideration for marginal twist angles. Interlayer hybridization effect is strong at Γ-valley for both P-and AP-bilayers. This is due to the contribution from d 2 z and p z atomic orbitals of metals and chalcogens respectively, that makes band edge value sensitive to interlayer distance variation. Piezoelectric potential plays an important role around corners of 2H in AP-bilayers for marginal angles. This promotes the band edge at Γ-valley to be located in corners of 2H domains. On the other hand, the piezoelectric potential is not as important for P-bilayers, making interlayer hybridization as the only important term to determine the band edge location, in this case at MX stacking. In the valence band at the K-valley, the lack of contribution from atomic orbtials of chalcogens makes it not so sensitive to variation of interlayer distance. Thus, the interlayer hybridization contribution does not vary significantly across the moiré supercell. For AP-bilayers, we see that the variation of the band edge at K-valley with twist angle for valence and conduction bands depends on the magnitude of the piezopotential. Therefore, piezoelectricity is the dominant effect for such a valley promoting the band edge to be located at XX and MM for the valence and conduction bands respectively at the K-valley. Ferroelectricity in P-bilayers is another important effect that influences the location of the band edge. The only TMD considered in this study where the interlayer hybridization effect is important at K-valley is P-WSe 2 . The band edge location at Q in AP-bilayers for marginal angles is mainly affected by the piezopotential. Just like in the case for the K-valley, piezoelectricity promotes the band edge to be located at MM for conduction bands. For P-bilayers, the band edge location at Q is mostly determined by the interlayer hybridization where ferroelectricity and piezopotential show a negligible contribution to MX /XM regions. The CB minimum is located along one-dimensional channels zig-zagging across the moiré superlattice. This value does not vary with the twist angle, which shows that interlayer hybridization is the dominant term.
Lastly, in Fig. S3 is shown a comparison of quantum dot potential depth in XX' region for K-valley holes 21 with their maximal kinetic energy in mini Brillouin zone for given twist angle. Regime of weakly coupled quantum dots is realized when the depth is much larger than the kinetic energy, which leads to the following estimate for twist angles: θ AP < 1 • .  S3. Dashed lines represent the depth of quantum dots in XX regions quantum for K-valley holes in AP-bilayers. Solid lines represent maximal kinetic energy of the hole given by E kinetic = (∆K) 2 /2m * K , where (∆K) 2 = 4πθ 2 AP /3a and m * K is the effective mass of the hole at K taken from Ref. 21 . Narrow bands are formed in the regime where quantum dot potential depth is much larger than the kinetic energy. This condition is satisfied in the interval θ AP < 1 • .
FIG. S4. Maps of band edge for AP-MoS 2 homobilayer in the VB at the Γ-and K-valley and CB at the K and Q -valley for θ AP = 0.2 • and θ AP = 3 • . See Fig. S2 for more details about 2H, XX and MM stacking configurations.
FIG. S11. Same as in Fig. S8 for P-WS 2