Polarization-driven band topology evolution in twisted MoTe2 and WSe2

Motivated by recent experimental observations of opposite Chern numbers in R-type twisted MoTe2 and WSe2 homobilayers, we perform large-scale density-functional-theory calculations with machine learning force fields to investigate moiré band topology across a range of twist angles in both materials. We find that the Chern numbers of the moiré frontier bands change sign as a function of twist angle, and this change is driven by the competition between moiré ferroelectricity and piezoelectricity. Our large-scale calculations, enabled by machine learning methods, reveal crucial insights into interactions across different scales in twisted bilayer systems. The interplay between atomic-level relaxation effects and moiré-scale electrostatic potential variation opens new avenues for the design of intertwined topological and correlated states, including the possibility of mimicking higher Landau level physics in the absence of magnetic field.

On the theory side, discrepancies in the Chern numbers were also found by two distinct approaches used to study the moiré electronic structures.The first approach involves deriving electronic structures from small unit cells containing local stacking arrangements [17][18][19][20][21][22].The second approach relies on density functional theory (DFT) calculations performed on reasonably sized moiré superlattices [23][24][25][26].Curiously, for tMoTe 2 , the Chern number of the topmost spin-up (spin-down) moiré valence band is found to be −1 (+1) within the local stacking approximation [20], whereas the DFT calculation conducted on a fully relaxed structure with a 3.89 • twist yield opposite Chern numbers [25].The latter is consistent with experimental observations [27].However, at smaller twist angles, the system size poses a substantial challenge to DFT calculations, and a direct comparison with experiments is currently unavailable.
In this Letter, we perform large-scale DFT calculations for tMoTe 2 and tWSe 2 down to 1.25 • twist angle.This is made possible by using a machine learning force field to obtain the relaxed structures, which enables a comprehensive exploration of the twist angle dependence of the moiré lattice reconstruction.We show that the observed difference in Chern numbers is due to the twist angle depen-dence of the moiré potential.Specifically, we find that as the twist angle varies, the location of the moiré potential maximum shifts from the MX stacking region to the XM stacking region (see Fig. 1 for the definition of MX and XM), causing a sign change of the Chern number.The shift of the moiré potential maximum is attributed to the competition between the in-plane piezoelectricity and the outof-plane ferroelectricity, a mechanism associated with the broken inversion symmetry in TMDs and absent in the local stacking approximation.The large-scale calculations, enabled by machine learning methods, also reveal multiple flat bands with Chern number all equal to 1 in tMoTe 2 at around 2 • twist, indicating the possibility of mimicking higher Landau level physics in the absence of magnetic field.The interplay between atomic-level relaxation effects and moiré-scale electrostatic potential variation opens new avenues for the design of intertwined topological and correlated states.
Following recent experiments [6,[13][14][15][16], we will focus on the valence bands of -type twisted TMD homobilayers.In these systems, the emergence of nontrivial band topology can be understood as a consequence of the real-space FIG. 1.The layer pseudospin skyrmion lattice.a, the -valley layer pseudospin () skyrmion lattice, using parameters from the local stacking approximation [20].The color denotes ∆  () and the arrow denotes ∆ , ().The black, green, and magenta dots denote the three high-symmetry local stacking sites, MM, XM, and MX, respectively.b, similar to a but using parameters from the DFT calculation [25].
where ∆ ∕ () and ∆  () are the intra-and inter-layer moiré potential, respectively.Due to the spin-valley coupling, the band edge at each valley is spin split, and opposite valleys carry opposite spins as required by time-reversal symmetry [28,29].The continuum Hamiltonian for the  ′ valley spin-down electrons can be obtained by applying timereversal symmetry to ℋ ↑  , resulting in moiré bands with opposite Chern numbers.
The moiré potential can be represented as an effective layer pseudospin magnetic field () = (Re ∆  , −Im ∆  , ).There are three high-symmetry local stackings in a moiré supercell, labeled as MM, XM, and MX (Fig. 1).It has been shown that () forms a skyrmion lattice with its north/south poles located at the MX and XM points [20].Curiously, for 3.89 • tMoTe 2 , using parameters from the local stacking approximation [20] and the DFT calculation [25], we find a reversal in the positions of the north/south poles between the two cases as shown in Fig. 1, which results in opposite skyrmion numbers [30,31].This contrast in skyrmion numbers, in turn, manifests as opposite Chern numbers for the topmost valence band, with only the DFT calculation matching the experiment.
Armed with the insight that the moiré potential landscape can affect the band topology, we now perform DFT calculations at even smaller twist angles.Because the system size at these angles (∼ 13,000 atoms at 1.25 • ) is beyond the typical scale of DFT relaxations, we first trained a neural network (NN) inter-atomic potential to capture the moiré lattice reconstruction.The NN potentials are parameterized by using the deep potential molecular dynamics (DPMD) method [32,33], where the training data is obtained from 5,000-6,000 ab initio molecular dynamics (AIMD) steps at 500 K for a 6 • twisted homobilayer calculated using the vasp package [34].We test the NN potential for a moiré bilayer at 5 • , and obtain a root mean square error of force < 0.04 eV/Å.More details can be found in Ref. [35].Figures 2(a The difference between reconstruction patterns at the large and small twist angles also affects the strain tensor distributions, shown in Fig. S5 of Ref. [35].We find that the shear strain (  ) and   −   are much larger than the normal strain (  +   ).As the twist angle decreases, the strains are mostly distributed near the domain boundaries due to the domain wall formation.These findings are consistent with previous calculations based on continuum model and parameterized interatomic potential [26,[36][37][38][39], as well as available experiments [40][41][42][43][44][45][46].
We then calculate the moiré band structure for the relaxed atomic structures.To reduce the computational cost, we adopt the siesta package [47] for band structure calculations.We first benchmark the accuracy of this local basis approach with the plane-wave basis approach by comparing the band structures at 6 • obtained from siesta and vasp, and reach a qualitative agreement between the two (see details in Ref. [35]).Then we perform small twist-angle band calculations by using the siesta package.Figures 2(c) and (d) show the twist-angle dependence of band structures for tWSe 2 and tMoTe 2 , respectively.The top valence bands consist of folded -valley and  ′ -valley minibands with opposite spins.A small band splitting can be seen, mostly between  and .Multiple factors, including trigonal warp-ing and intervalley coupling, may contribute to the splitting.Nevertheless, we assume approximate spin  conservation, and separate moiré bands originating from the two valleys by adding a small Zeeman field in the calculation.In the following, we shall focus on the moiré bands from the spin-up  valley.
To determine the Chern numbers for the moiré bands, we first calculate the eigenvalues of the DFT wave functions under three-fold rotational symmetry ( 3 ).The Chern number is then determined by the product of  3 eigenvalues at rotationally invariant momenta [48]: exp( 2 3 ) = −      ′ , where the 's are the  3 eigenvalue at the high symmetry point of the moiré Brillouin zone (mBZ).These eigenvalues are labeled in Fig. 2. Our assignment of the Chern numbers for tWSe 2 at 3.15 • and 1.70 • agree with a recent calculation at 2.28 • in which the Chern numbers were calculated by directly integrating the Berry phase over the entire mBZ [26].The Chern numbers for tMoTe 2 are further confirmed through the integration of the Berry curvature within the mBZ with the help of Wannier interpolations.The change of the Chern numbers with varying twist angles can be understood by tracking the evolution of the  3 eigenvalues, which signals band inversion.For example, in tWSe 2 [Fig.2(c [6].Further increasing the twist angle results in trivial insulators up to 1.6 • [6].Remarkably, all these observations are consistent with the trend in the twist-angle dependence of our calculations, confirming the validity of our machine-learning based approach.Our calculation also predicts that in tMoTe 2 , as the twist angle decreases, the Chern numbers of the two topmost bands change from (+1, −1) to (+1, +1), and finally to (0, 0) at the smallest angle of the calculation.In particular, our calculations reveal multiple flat bands with Chern number all equal to +1 at around 2 • , indicating the possibility of mimicking higher Landau-level physics in the absence of a magnetic field.The presence of multiple bands of Chern number +1 has been confirmed by a recent experimental measurement [49].
Since we are interested in the sign change of the Chern number of the topmost band, in the following we will focus on tWSe 2 .As mentioned earlier, the evolution of band topology in momentum space is closely related to the change in the real space moiré potential.In particular, the location of the north/south poles of (), which directly affects the skyrmion number, is given by the difference between the moiré potentials at the top and bottom layer.
The moiré potential can be inferred from the surface Hartree potential [50], defined as the difference between the Hartree potential above and below the twisted bilayer surface in DFT calculations.Figure 3(a) shows the coarsegrained surface potential drop at 3.15 • in tWSe 2 .More details can be found in Ref. [35].The maximum is located at MX, zero at MM, and minimum at XM. Surprisingly, the surface potential drop shows a sign reversal at XM (and MX) as the twist angles decrease [see Fig. 3(a-d)].Going from 3.15 • , to 1.70 • , 1.47 • , and eventually down to 1.25 • , the potential at the high-symmetry point XM (and MX) changes sign, and the area of the flipped region grows in size.This sign switch suggests that the north pole of () at ∼ 3 • becomes the south pole at ∼ 1 • .Additional features can be identified near the MM site, where the surface potential drop mimics the pattern of a six-petal flower with  3 symmetry.We find that the amplitude of potential inside the petal is comparable with that at XM (and MX) at 1.70 • , suggesting unique quantum confinement effects which reshape the electron wave function.The overall effects of these features can be clearly seen by a line cut along MM-XM-MX-MM, showing rich variations and multiple extremes in Fig. 3(e).The intricate behavior of the surface potential goes beyond the continuum approximation of moiré potential based on the first-star expansion of the reciprocal lattice vectors alone, evident by their Fourier transform as shown in Fig. 3(f).
The evolution of the surface potential implies that the layer polarization of the wave functions should also change with the twist angle.In Fig. 4, we plot the real-space wave function for the two topmost bands at the  point of the mBZ at various twist angles for tWSe 2 .Since the timereversal symmetry  and  2 symmetry are preserved at the  point, only the wave function in the top layer is plotted, while the wave function in the bottom layer can be obtained by performing a  2 operation, under which MX/XM in the top layer is mapped to XM/MX in the bottom layer.At 3.15 • , the wave function of the first band in the top layer is localized at MX.As the twist angle decreases, the localization region switches to MM and eventually to XM.The shift of the wave function location also coincides with the change of the  3 eigenvalue at the  point [Fig.2(c)].Two remarks are in order.First, the switch from MX to XM of the first band indicates the flip of the layer pseudospin, which gives rise to the sign change of the Chern number as shown in Fig. 2(c).Second, aside from orbitals located at MX and XM, those at MM also play a significant role in deciding band topology, evident by the distribution of the wave functions (Fig. 4).Thus for a tightbinding model to properly describe the moiré band topology of tWSe 2 , one also needs orbitals from the MM site [51].This goes beyond the real-space skyrmion picture discussed earlier.
What is the origin behind the change in surface potential drop as a function of twist angle?For a 2D bilayer system in global charge neutrality, the electrostatic surface potential can be directly associated with the interlayer electric polarization.In twisted TMD homobilayers, two microscopic mechanisms contribute to this polarization: ferroelectricity and piezoelectricity.The moiré ferroelectricity arises from the inversion symmetry breaking in -type TMD bilayers [52,53].In a moiré supercell, this leads to alternating out-of-plane ferroelectric polarization depending on the local stacking registry, with opposite dipoles in the XM/MX region [53][54][55][56].On the other hand, since monolayer TMDs lack inversion symmetry, the strain field can produce piezoelectric polarization for each layer [57].Because the two layers have opposite patterns of atomic displacement fields and the same piezoelectric coefficient, the polarization charge distributions are opposite between the two layers, which can produce a vertical potential drop [38,39].As pointed out in Ref. [39], these two types of polarization charges can be opposite in sign, and their competition will determine the potential drop.
Note that additional in-plane polarizations can also arise from the out-of-plane ferroelectricity and local symmetry breaking [58].While these effects are present in our DFT calculations, they do not create additional Hartree potential that changes the Chern numbers.
Our machine-learning based first-principles simulations enable us to quantitatively study the polarization effects from the relaxed structure with atomic resolution.We start by separately calculating the piezoelectric charge and ferroelectric charge.For each layer, the piezoelectric charge density is directly proportional to the gradient of the strain field by  piezo = − ẽ11 [  (  −   ) − 2    ], where ẽ11 is the independent non-zero component of the piezoelectric tensor [38,57].The out-of-plane ferroelectric polarization is obtained directly from integrating the charge density along the  direction for each local-stacking unit cell.More details about the polarization calculations can be found in Ref. [35].In Fig. 5, we plot the piezoelectric (further considering the dielectric screening effect [38]) and ferroelectric charge densities, as well as their sum at 3.15 • and 1.25 • , respectively.The total polarization charge matches qualitatively the pattern of the surface potential drop.
Specifically, at 3.15 • , the piezoelectric charges are mainly distributed in the XM and MX regions because of the larger gradient of shear strain and   −   in these regions (see Fig. S5 of Ref. [35]).The piezoelectric charge is negative at MX and positive at XM for the top layer.Because the bottom layer has the opposite atomic displacements, it has the opposite charge distributions.In contrast, the ferroelectric charge is positive (negative) at MX and negative (positive) at XM for the top (bottom) layer.Adding them together, we find the total charge density is negative (positive) at MX and positive (negative) at XM for the top (bottom) layer.As the twist angle decreases to 1.25 • , the ferroelectric charge density at MX and XM remains virtually unchanged, but the total amount of ferroelectric charge within the MX and XM domains increases following the formation of the domain wall.In contrast, because the shear strain and   −   are mainly distributed along the domain wall and are uniformly small inside the XM and MX domains (see Fig. S5 of Ref. [35]), the piezoelectric charge density peaks near the domain wall but decreases at the interior of the domain.This explains the six-petal flower pattern that we discovered for the surface potential drop.As a consequence, the total charge density is now positive (negative) at MX and negative (positive) at XM for the top (bottom) layer.This trend from 3.15 • to 1.25 • is consistent with the variation of the surface moiré potentials and wave functions.In contrast, within the local stacking approximation, the sign of polarization charge is fixed and a sign reversal of the Chern number is not possible.Probing the predicted reversal of the polarization charges at MX and XM should be a clear experimental evidence of the change of band topology in momentum space.
Up to this point, our discussion has focused on tWSe 2 .While the evolution of the moiré potential in tMoTe 2 follows a similar trend, there are some quantitative differences in the band structures between tWSe 2 and tMoTe 2 .This is mostly due to a couple of factors.WSe 2 has a lighter effective mass (0.35) compared to MoTe 2 (0.62), resulting in a larger bandwidth.In addition, differences in both elastic and piezoelectric coefficients also lead to quantitative changes in the moiré potential between these materials.More details can be found in Ref. [35].
In summary, we have performed large-scale DFT calculations on -type TMD homobilayers.Our results demonstrate machine learning as a powerful tool to study moiré systems, by revealing the importance of lattice relaxation that eventually leads to qualitative changes in band topology.This change is attributed to the competition between in-plane piezoelectricity and out-of-plane ferroelectricity, resulting in electrostatic potential variation that reshapes the potential landscape for any moiré electronic states.Our findings highlight the crucial long-range interactions arising from polarization charges, which changes rapidly as the twist angle decreases.This behavior necessities the explicit calculations of moiré electronic potential even at minimal twist angle.
We acknowledge useful discussions with Jiaqi Cai, Ken Shih and Xiaodong Xu.The machine learning is supported by the discovering AI@UW Initiative and DMR-2308979, the band structure calculation is supported by DOE Award No. DE-SC0012509, and the calculation of the piezoelectric and ferroelectric charge is supported by ARO MURI with award number W911NF-21-1-0327.This work was facilitated through the use of advanced computational, storage, and networking infrastructure provided by the Hyak supercomputer system and funded by the University of Washington Molecular Engineering Materials Center at the Uni-versity of Washington (DMR-2308979).Access to GPUs is provided by Microsoft sponsorship of Azure credits to the Department of Materials Science and Engineering at UW.Note added.-Werecently became aware of two related works in which machine learning force field is also used to calculate the moiré bands of twisted MoTe 2 homobilayers [59,60].
FIG. 2. Lattice relaxations and band structures.a, the in-plane displacement field of the top-layer W atoms at 3.15 • and 1.25 • for tWSe 2 .The color and arrow denote the amplitude and direction of in-plane displacement fields, respectively.b, the interlayer distance (ILD) distribution at 3.15 • and 1.25 • for tWSe 2 .c and d, twist-angle dependence of the valence moiré bands of tWSe 2 and tMoTe 2 , respectively.The labels indicate the spin orientations and the  3 eigenvalues at high symmetry point with  =  ∕3 ,  * =  −∕3 and 1 = −1.The  3 eigenvalue is the same at the  and the  ′ point.

FIG. 3 .
FIG. 3. Evolution of surface moiré potential.a-d, twist angle dependence of the difference of DFT Hartree potential, ∆  , between the top-layer surface and the bottom-layer surface in tWSe 2 .The white dashed line denotes the path MM-XM-MX-MM.e, the distribution of ∆  along the dashed line in a for different twist angles.f, the Fourier components of ∆  for different twist angles. 1 is the length of the moiré reciprocal lattice vector.Inset in f: schematics of the Hartree potential drop between the top-layer surface and bottom-layer surface.
) and (b) show the calculated in-plane displacement field of the top-layer W atoms and interlayer distance in tWSe 2 at 3.15 • and 1.25 • .The displacement field in the bottom layer shows the opposite pattern.It is clear that while the local stacking varies smoothly at 3.15 • , at 1.25 • large domains of the MX and XM regions form, with domain walls connecting the shrunken MM region.
FIG. 4. Evolution of wave functions.a-d, twist-angle dependence of the wave function of the first band (top row) and second band (bottom row) at the  point for the top-layer in tWSe 2 , respectively.Here, we map the weight of the projected wave function onto the W atomic orbitals.The value is normalized to its maximum in each plot.The dashed parallelogram denotes a moiré unit cell.

FIG. 5 .
FIG.5.The competition between the in-plane piezoelectric and out-of-plane ferroelectric polarizations.a, the piezoelectric charge density, the ferroelectric charge density, and the total charge density, respectively, at 3.15 • in tWSe 2 .b, the same quantities at 1.25 • .