# Interacting phase diagram of twisted bilayer MoTe2 in magnetic field

Minxuan Wang  National High Magnetic Field Lab, Tallahassee, FL 32310
Department of Physics, University of California, Berkeley, CA 94720, USA
Xiaoyu Wang  National High Magnetic Field Lab, Tallahassee, FL 32310  Oskar
Vafek  National High Magnetic Field Lab, Tallahassee, FL 32310  Department of
Physics, Florida State University, Tallahassee, FL 32306

###### Abstract

We study electron-electron interaction induced states of twisted bilayer MoTe2
in an out-of-plane magnetic field $B\hat{\bf z}$ near one hole per moiré unit
cell filling. The 3D phase diagram showing the evolution of competing phases
with $B$, interaction strength and an out-of-plane electric field is presented
at electron fillings that follow the Diophantine equation along Chern number
$-\text{sign}\left(B\right)$ line, that is pointing away from the charge
neutral filling, where we find prominent Chern insulators consistent with the
experiments. We also explain the experimental absence of prominent Chern
insulators along the Chern number $+\text{sign}\left(B\right)$ line.

_Introduction.—_ Moiré materials are well known for their tunability and
hosting various correlated phases due to their narrow energy electronic
bands[1, 2]. Among these materials, twisted transition metal dichalcogenides
(TMD) homobilayers, as well as heterobilayers, are predicted[3, 4], and
experimentally shown, to exhibit integer quantum anomalous Hall (QAH) effect
[5, 6]. Intriguingly, fractional quantum anomalous Hall effect has recently
also been observed in twisted MoTe2 homobilayer (tMoTe2) [7, 8, 9, 10], with a
number of theoretical papers devoted to these integer and fractional Chern
insulators[11, 12, 13, 14, 15, 16, 17, 18, 19, 20, 21, 22, 23, 24, 25, 26,
27]. Due to the large unit cell area $A_{uc}$ of moiré materials, out-of-plane
magnetic field $B\hat{\bf z}$ has proven to be a powerful probe of their
correlated states [28, 29, 7, 8]. Particularly interesting is the evolution of
such states in the $\nu$-$B$ plane, where the electron filling
$\nu=\rho_{e}A_{uc}$ and $\rho_{e}$ is the areal density of electrons. The key
tool is the thermodynamic formula [30, 31] which states that the quantized
Hall conductivity of a gapped state of electrons (with charge $-e$) at $T=0$
follows

$\sigma_{xy}=-\sigma_{yx}=-e\frac{\partial\rho_{e}}{\partial
B}\bigg{|}_{\mu,T}=-t\frac{e^{2}}{h}.$ (1)

Here $t$ is the Chern number[32] and we use the SI units throughout. The
Diophantine equation for $\nu$ is[33, 34, 35, 36, 37]

$\nu=s+t\frac{\phi}{\phi_{0}},$ (2)

where the flux quantum $\phi_{0}=h/e=4136\mathrm{T\cdot nm}^{2}$ and the flux
through the unit cell is $\phi=BA_{uc}$.

A ubiquitous experimental observation in tMoTe2 at twist angles in the
$3.4^{\circ}\sim 3.9^{\circ}$ range is the QAH state at $\nu=-1=s$ evolving
into a Chern state at $B\neq 0$ with a prominent gap at $t=-\text{sign}(B)$
and a noted absence of a prominent gap at $t=+\text{sign}(B)$ [7, 8, 9, 10].
In other words, the Streda lines emanating from $\nu=-1$ point away from the
charge neutrality point (CNP) for either sign of $B$. The two QAH states with
opposite Chern numbers at $B=0$ are partners under the time reversal symmetry,
which is spontaneously broken, and therefore they must be equally stable at
$B=0$. The relative instability of the $t=+\text{sign}(B)$ state at $B\neq 0$
has thus far not been explained. Note that in $1.23^{\circ}$ twisted WSe2
bilayer the prominent Chern insulator line at $s=-1$ points towards the
CNP[29], and that both behaviors have been observed in the magic angle twisted
bilayer graphene device of Ref.[38], as well as in calculations [39].

Here we show that, within our $B\neq 0$ self-consistent Hartree-Fock
calculation at the twist angle $3.89^{\circ}$ and the dielectric constant
$\epsilon$ entering the interaction (4) chosen to reproduce the $B=0$ QAH
state[23], the ground state at $s=-1$ and $t=-\text{sign}(B)$ is valley
polarized and that the holes preferentially populate the valley whose spin is
aligned to the direction of the magnetic field. This state can be thought of
as a Landau quantized QAH state with a non-vanishing (large) gap in the
$B\rightarrow 0$ limit and we refer to it as a Chern paraelectric because it
preserves $C_{2y}\mathcal{T}$, the combination of an in-plane 2-fold rotation
symmetry and the time reversal symmetry. Its energy spectrum as a function of
$\phi/\phi_{0}$ is shown in the Fig.2(a). Since the spin orientation of the
holes in the two valleys is opposite due to the spin-valley locking[40], the
valley polarization reverses upon reversing the sign of $B$ (still at
$t=-\text{sign}(B)$). Interestingly, the nature of the lowest energy state at
$s=-1$ and $t=+\text{sign}(B)$ depends sensitively on $g$, the strength of the
spin and orbital Zeeman coupling[41]. At $g=0$, it is valley polarized but now
the holes populate the valley whose spin is anti-aligned to the direction of
the magnetic field. For a given $B$, it can be thought of as a Landau
quantized time reversed partner of the QAH discussed above. Its (large) gap
also does not vanish in the $B\rightarrow 0$ limit, see Fig.2(d). If it were
stable, the valley polarization would reverse at a fixed $B$ upon changing the
filling from $-1-\frac{|\phi|}{\phi_{0}}$ to $-1+\frac{|\phi|}{\phi_{0}}$,
similar to the observations reported in Ref.[42] in graphene heterostructures.
Such a state would appear as a prominent Chern insulator, inconsistent with
the experiments[7, 8, 9, 10]. At $g\gtrsim 2$, however, the ground state
changes and the holes populate the valley with spin aligned to the direction
of the magnetic field for a range of $B$ that expands to lower $|B|$ as $g$
increases (Fig.2(h)). It can be thought of as an integer quantum Hall (IQH)
state, populating the lowest pair of Landau levels (LLs) of the Chern
paraelectric excitation spectrum at $t=-\text{sign}(B)$ with electrons, as
seen in the Fig.2(c). Because of the sizable band mass, its gap is small and
vanishing as $B\rightarrow 0$, consistent with the experimental absence of a
prominent Chern insulator at $t=+\text{sign}(B)$ and presence at
$t=-\text{sign}(B)$. As detailed below, we estimate the orbital Zeeman
contribution to $g$ in tMoTe2 to be $3.23$ and of the same sign as the spin
contribution $2$, placing us safely in the $g\gtrsim 2$ regime.

Perpendicular electric field, $D$, has also proved to be a powerful tool to
probe and manipulate the nature of correlated states in tMoTe2[9, 8, 7, 10].
In the experiments at $B=0$, the spin polarized QAH state transitions into a
trivial insulator above a critical $D_{*}$[9, 8, 7, 10]. The proximate trivial
state is argued to be spin, and therefore valley, unpolarized in the Ref. [8]
and spin, and therefore valley, polarized in the Ref. [10]. Such topological
transition is also seen in the Hartree-Fock calculation at $B=0$[23, 27], and
the spin structure of the trivial state depends on the value of $\epsilon$.
For $\epsilon>\epsilon^{\prime}\approx 19$ the transition is directly into a
spin unpolarized intervalley coherent state (IVC) and for $\epsilon_{*}\approx
15.4<\epsilon<\epsilon^{\prime}$ the proximate trivial state is valley
polarized (VP). This can be seen in the $B=0$ plane of Fig.4 where we vary the
$D$-induced potential difference between two layers $u_{D}$. Intriguingly, at
$B\neq 0$ and $t=-\text{sign}(B)$ a sharp change of compressibility can be
seen in the experiment at the same $D_{*}$ as the transition at $B=0$ (see
Extended Data Fig.10 of Ref. [8]). Our calculation at $B\neq 0$ captures this
effect. As shown in the Fig.3(a), at $\epsilon=20$ chosen to recover the
phenomenology of Ref.[8] at $B=0$, we observe a crossover at $D_{*}$ between a
Chern paraelectric with a large gap and an ostensibly IQH state of a Landau
quantized trivial insulator with one LL depopulated by electrons and a small
gap. We refer to the latter as a QH ferroelectric because at $u_{D}=0$ it
breaks $C_{2y}\mathcal{T}$. Both states are spin polarized. The results of a
separate Hartree-Fock calculation at $t=0$ shown in the Fig.3(b) confirm the
prominent trivial insulator with a large gap. Thus, the incompressible states
can appear simultaneously at $\nu=-1$ and $\nu=-1-\frac{|\phi|}{\phi_{0}}$ for
a range of $D$ near $D_{*}$.

In order to explore the robustness of the above results, as well as the
tunability of the correlated phases, we perform the Hartree-Fock calculations
at $s=-1$ and $t=-\text{sign}(B)$ at varying $\epsilon$, $u_{D}$ and $B$. The
resulting tentative phase diagram is shown in the Fig.4. At $D=0$ the first
order topological transition at $\epsilon_{*}$ extends to $B\neq 0$ despite
separating states with the same Chern number $t$ at $B\neq 0$ because the QH
ferroelectric on one side breaks $C_{2y}\mathcal{T}$ and the Chern
paraelectric on the other does not. This first order phase transition moves
towards stronger interaction as $B$ increases, favoring the Chern paraelectric
i.e. the descendant of the QAH state. The transition must also extend to small
non-zero $D$ despite the explicit breaking of $C_{2y}\mathcal{T}$ by $D$
because it starts out first order at $D=0$. Indeed, Fig.3(c) shows the jump in
the layer polarization at a fixed $B$ persisting to non-zero $u_{D}$, with $D$
favoring the ferroelectric and the first order transition terminating at a
critical endpoint. As we show below, we can describe this using a simple
Landau theory for an Ising ferroelectric with a negative quartic coupling.
Upon varying $B$ the critical endpoint extends to the red curve shown in the
Fig.4. Interestingly, the $D$ induced crossover shown in the Fig.3(a) appears
at $B$ and $\epsilon$ notably separated from the critical endpoint. We
understand it as crossing the Widom line, defined as the peak in the layer
polarizability, extending beyond the critical endpoint as shown in the
Figs.3(d) and 4. As IVC is strongly suppressed by the $B$-field, Chern
paraelectric and QH ferroelectric are the main phases in the Fig.4. Below we
discuss our formalism and computational scheme that lead to the above results.
We note in passing that at the twist angle $3.89^{\circ}$, $\phi=\phi_{0}$ at
$B\simeq 176$T, so to get below $9$T as in the experiments we need
$\phi/\phi_{0}\lesssim 1/20$. Reaching such a small flux efficiently in our
interacting calculation is enabled by a recently developed technique[43, 44,
45] utilizing the $B=0$ hybrid Wannier states to construct the basis at $B\neq
0$; some of our calculations go down to $\phi/\phi_{0}=1/47$.

Figure 1: (a) Band structure at $B=0$ in the valley $\mathbf{K}$ for tMoTe2 at
$\theta=3.89^{\circ}$.(b,c) Hofstader spectra at $\mathbf{K}$ and
$\mathbf{K}^{\prime}$, omitting Zeeman effect. Figure 2: (a-d)Hartree Fock
energy spectrum with varying magnetic field at
$\nu=-1-\frac{|\phi|}{\phi_{0}}$ and $\nu=-1+\frac{|\phi|}{\phi_{0}}$.
Electrons occupy all eigenstates below black dashed line (guide to eye). Blue
(red) points represent eigenstates in $\mathbf{K}$($\mathbf{K}^{\prime}$)
valley. (e,f)Total energy per unit cell versus magnetic field. Dashed lines
are obtained from quadratic interpolation using calculated energy per unit
cell at $B=0$. (g,h)Energy per unit cell with non-zero $g$ factor.

_Continuum model—_ The moiré bands of experimental relevance originate from
the spin-valley locked valence bands of monolayer MoTe2[46]. The relations in
Eq. (1) and Eq. (2), as well as the stated experimental results, hold in any
right handed coordinate system (even though $\sigma_{xy}$, and therefore $t$,
changes sign between two coordinate system choices related by a $180^{\circ}$
rotation about $x$ or $y$ axis). We are therefore free to choose any of the
two out-of-plane directions as $\hat{\bf z}$ and we made the common choice
[46] of aligning it with the spin angular momentum direction of the valley
$\mathbf{K}$. Then, in valley $\mathbf{K}$ and for spin up, the non-
interacting moiré band structure can be described by a continuum electronic
model [3]:

$\begin{split}&{H}_{\mathbf{K}}=\int\mathrm{d}^{2}\mathbf{r}\begin{pmatrix}\psi^{\dagger}_{b\mathbf{K}}(\mathbf{r})&\psi^{\dagger}_{t\mathbf{K}}(\mathbf{r})\end{pmatrix}\times\\\
&\begin{pmatrix}{h}_{b}({\mathbf{p}})+\Delta_{b}(\mathbf{r})+\frac{u_{D}}{2}&\Delta_{T}(\mathbf{r})\\\
\Delta_{T}^{*}(\mathbf{r})&{h}_{t}({\mathbf{p}})+\Delta_{t}(\mathbf{r})-\frac{u_{D}}{2}\end{pmatrix}\begin{pmatrix}\psi_{b\mathbf{K}}(\mathbf{r})\\\
\psi_{t\mathbf{K}}(\mathbf{r})\end{pmatrix}.\end{split}$ (3)

Here $\psi^{\dagger}_{l\mathbf{K}}(\mathbf{r})$ creates an electron in layer
$l=b,t$ , in valley $\mathbf{K}$, and at position $\mathbf{r}$.
${h}_{l}({\mathbf{p}})=-({\mathbf{p}}-\hbar\kappa_{l})^{2}/2m^{*}$ is kinetic
energy, with ${\mathbf{p}}=-i\hbar\nabla_{\mathbf{r}}$ and effective mass
$m^{*}$. $\kappa_{l}$ points sit at the corners of the hexagonal moiré
Brillouin zone, we describe the effect of $D$ by adding a potential difference
$u_{D}$ between two layers,
$\Delta_{b(t)}(\mathbf{r})=2V\sum_{j=1,3,5}\cos(\mathbf{G}_{j}\cdot\mathbf{r}\pm\phi)$
is the intralayer potential for bottom(top) layer, and
$\Delta_{T}(\mathbf{r})=w(1+e^{-i\mathbf{G}_{2}\cdot\mathbf{r}}+e^{-i\mathbf{G}_{3}\cdot\mathbf{r}})$
is interlayer tunneling potential. The reciprocal lattice vector
$\mathbf{G}_{j}$ is obtained by counterclock-wise rotation of
$\mathbf{G}_{1}=4\pi/(\sqrt{3}a_{M})\hat{\mathbf{y}}$ by angle $(j-1)\pi/3$.
$a_{M}=a_{0}/\theta[\text{rad}]$ is moiré superlattice constant, where $a_{0}$
is the lattice constant of monolayer MoTe2. We choose parameters
$(a_{0},m^{*},V,\phi,\omega)=(0.355{\rm nm},0.62m_{e},17.0{\rm
meV},+107.7^{\circ},-16.0{\rm meV})$ at twist angle $\theta=3.89^{\circ}$
[23]. The non-interacting model in the $\mathbf{K}^{\prime}$ valley and for
spin down is related by time reversal symmetry $\mathcal{T}$, which takes the
complex conjugation of the $2\times 2$ Hamiltonian in Eq. (3).

In an external magnetic field, we replace
${\mathbf{p}}\rightarrow{\mathbf{p}}+e\mathbf{A}$, where
$\nabla\times\mathbf{A}=B\hat{\bf z}$, and add Zeeman coupling
$g\mu_{B}B\tau_{0}s_{z}$, where $\tau_{0}$ is the identity Pauli matrix acting
on layer degrees of freedom, and $s_{z}$ is the spin Pauli matrix. At $B=0$,
the Hamiltonian is invariant under three-fold out-of-plane rotation $C_{3z}$,
in-plane 2-fold rotation symmetry about the $y$-axis $C_{2y}$ and the time
reversal symmetry $\mathcal{T}$. At this order the continuum Hamiltonian also
has a pseudo-inversion symmetry[19, 47] $\mathcal{I}$ given by
$\mathbf{r}\rightarrow-\mathbf{r}$ in each valley and flipping the layers with
$\tau_{1}$. While both $C_{2y}$ and $\mathcal{T}$ are broken at $B\neq 0$,
their product $C_{2y}\mathcal{T}$ is preserved, as are $C_{3z}$ and
$\mathcal{I}$.

Figure 3: (a)Hartree Fock spectrum at $\epsilon=20$ with varying $u_{D}$ at
$\nu=-1-\frac{|\phi|}{\phi_{0}}$ and (b) $\nu=-1$. Topological phase
transition at $B=0$ is at $u_{D^{*}}$. (c,d) Layer polarization $P_{l}$ and
polarizability with varying $u_{D}$ and $\epsilon$ at $B=8.78$T.

Fig. 1(a) shows the non-interacting band structure at $B=0$ in the valley
$\mathbf{K}$ where the two uppermost bands carry Chern number $+1$ and $-1$.
The band dispersion and Chern numbers in the opposite valley
$\mathbf{K}^{\prime}$ are related by $\mathcal{T}$. For $B>0$ the
corresponding Hofstadter spectra for the valleys $\mathbf{K}$ and
$\mathbf{K}^{\prime}$, together with the $(s,t)$ labels of the gaps, are shown
in Fig. 1(b) and (c) respectively. As can be seen by the evolution of the
magnetic subbands, a Chern band gains or loses states per unit cell according
to the relation $1+t\phi/\phi_{0}$. To study correlated phenomena, we project
the Coulomb interaction onto the upper two Chern $\pm 1$ bands of both valleys
if at $B=0$, i.e. two bands per valley, or the magnetic subbands emanating
from these Chern bands if at $B\neq 0$,

${H}_{V}=\frac{e^{2}}{4\epsilon\epsilon_{0}A}\sum_{\mathbf{q}\neq
0}\frac{1}{|\mathbf{q}|}\tanh\left(|\mathbf{q}|d\right):\mathcal{P}{\rho}_{\mathbf{q}}{\rho}_{-\mathbf{q}}\mathcal{P}:.$
(4)

Here
$\rho_{\mathbf{q}}=\sum_{\eta,l}\int\mathrm{d}^{2}\mathbf{r}e^{-i\mathbf{q}\cdot\mathbf{r}}\psi_{l\eta}^{\dagger}(\mathbf{r})\psi_{l\eta}(\mathbf{r})$
is the Fourier transform of the electron density operator, with
$\eta=\mathbf{K}$ or $\mathbf{K}^{\prime}$, and total area $A$. In this work
we consider a dual gate screened Coulomb interaction, $d$ being the distance
between each gate and tMoTe2, taken to be $30$nm in our calculations, with the
small distance between the two MoTe2 layers neglected. Symbol $::$ denotes
operator normal ordering. $\mathcal{P}$ projects onto the mentioned Hilbert
subspace which can be generated at $B\neq 0$ by solving the non-interacting
problem[44, 43, 29]. At the $B$ values of interest here, expanding
$\mathcal{P}$ in the LL basis proves to be computationally expensive when
dealing with interactions. Instead, we construct the basis states using the
$B=0$ hybrid Wannier states method following Ref. [44, 45].

_Streda line near $\nu=-1$ and Zeeman effect—_ The $B>0$ energy spectrum of
candidate Chern states along the $(-1,\mp 1)$ lines is shown in Fig. 2(a-d).
We present their total energy per unit cell as a function of $B$ at varying
Zeeman $g$ factors in Fig. 2(e-h). The orbital contribution to the Zeeman
effect stems from $d$-orbital nature of the states in each valley[46]. The
two-band $\mathbf{k}\cdot\mathbf{p}$ Hamiltonian describing the large
semiconductor gap in MoTe2 monolayer has the form of a massive 2+1D Dirac
Hamiltonian[46], whose minimal coupling to a low $B$-field results in an
additional $\frac{1}{2}\hbar\omega_{c}$ energy shift compared to a Schrodinger
Hamiltonian. Since the Dirac mass effectively flips sign between different
valleys, this orbital shift also flips sign; it adds to the usual spin Zeeman
factor of $2$ (see Ref. [45] for details). Using the known effective mass of
MoTe2 we find it to be $3.23$ and thus the total $g=5.23$. Using the computed
ground state energy of the QAH state at $B=0$, we smoothly interpolate between
the calculated energies at $B\neq 0$ and $B=0$ at $g=0$ using at most a second
order polynomial in $B$ as shown by the dashed lines in Fig.2(e,f). Because
the Chern paraelectrics are fully valley polarized, their Zeeman energy can be
obtained analytically. We can therefore find the critical $B$ below which the
valley switches for Chern paralectrics. As seen in the Fig.2(g,h) the crossing
only occurs for the $(-1,1)$ line. While at $g=2$ it occurs at $5.01$T (cyan
star), at the more realistic $g=5.23$ it is $\sim 178$mT (red star). Thus any
putative density induced valley switching would occur in the regime where the
effects of quenched disorder dominate in realistic devices[9, 10], likely
modifying the clean limit energetics obtained here

Figure 4: Tentative phase diagram of tMoTe2 at
$\nu=-1-\frac{|\phi|}{\phi_{0}}$ and $3.89^{\circ}$ twist. The blue surface at
$B\neq 0$ denotes the boundary between IVC state and two valley $U(1)$
preserving states, Chern paraelectric and QH ferroelectric.

_Effects due to $D$-field—_ Energy spectrum induced by the perpendicular
electric field $D$ at $\phi/\phi_{0}=1/20$, corresponding to $B=8.78T$, is
shown in Fig. 3(a) along the $(-1,-1)$ line. The lowest unoccupied LL can be
seen to detach from the unoccupied LL group as $D$ increases and to closely
approach the highest occupied LL. This is similar to the magnetic subband
spectrum of an electron confined to move in a plane in a periodic potential
and an out-of-plane uniform magnetic flux slightly smaller than $\phi_{0}$
shown in the Fig. 7 of Ref.[48] as the strength of the periodic potential
increases and drives a transition from a Chern insulator (a broadened LL) and
a trivial insulator at their $\phi=\phi_{0}$. As shown in the Fig. 3(b) for
$(-1,0)$, a prominent trivial insulator is indeed stabilized in our Hartree-
Fock calculation as well.

The evolution of the layer polarization
$P_{l}=\frac{A_{uc}}{A}\sum_{\eta}\int\mathrm{d}^{2}\mathbf{r}\langle\psi_{t\eta}^{\dagger}(\mathbf{r})\psi_{t\eta}(\mathbf{r})-\psi_{b\eta}^{\dagger}(\mathbf{r})\psi_{b\eta}(\mathbf{r})\rangle$
and polarizability $dP_{l}/du_{D}$ at $(-1,-1)$ with varying $\epsilon$ and
$D$, at a fixed $B=8.78T$, are shown in Fig. 3(c,d). The main features can be
captured by a simple Landau theory for an Ising order parameter, odd under
$C_{2y}\mathcal{T}$,

$\begin{gathered}\delta E=-\gamma
u_{D}P_{l}+\frac{a(\epsilon)}{2}P_{l}^{2}-\frac{b}{4}P_{l}^{4}+\frac{1}{6}P_{l}^{6}.\end{gathered}$
(5)

Here $\gamma$ and $b$ are positive constants, and the prefactor of $P_{l}^{6}$
was used to rescale the energy units. Decreasing $a$ with decreasing
$\epsilon$ causes the spontaneous symmetry breaking at $u_{D}\propto D=0$ via
a 1st order phase transition due to the negative quartic term. At the critical
endpoint, the local minima merge (see Ref.[45]).

_Phase diagram—_ The tentative 3D phase diagram along the $(-1,-1)$ line is
shown in the Fig. 4. Markers correspond to parameters at which the transition,
or a crossover in the case of Widom line, was calculated within Hartree Fock.
The solid lines and surfaces are interpolated based on physical arguments. At
$u_{D}\propto D=0$ the Chern paraelectric and the quantum Hall ferroelectric
are distinguished by $P_{l}$, a Landau order parameter, and so must be
separated by a phase transition. At $u_{D}\neq 0$ however, it is possible to
transform one to another without encountering a phase transition by avoiding
the (orange) surface of 1st order phase transitions terminating in the (red)
critical endpoint curve. The IVC is a sharply defined phase at any $B$ and $D$
because it breaks valley $U(1)$ symmetry, preserved by $B$ and $D$ while the
other phase(s) do not break it. It is suppressed by $B$ even at $g=0$, and
further suppressed by the Zeeman effect. Therefore, while all the lines are
computed at $g=5.23$, the (cyan) phase boundary at $u_{D}=16.8$meV is
extrapolated from the $g=2$ markers; it moves closer to the $B=0$ plane at
higher $g$ (see Ref.[45]).

_Discussion—_ Our calculation explains the presence of the $(-1,-1)$ line and
the absence of the $(-1,+1)$ line in the experiments[7, 10, 8, 9]. It also
demonstrates the delicate balance between the energies of competing Chern
states at $(-1,+1)$ polarized to opposite valleys, with the orbital Zeeman
contribution ultimately deciding on the ground state with the small gap.
Recent calculations[49, 50, 47] pointed out a rich evolution of the Chern
bands and a lattice relaxation playing increasingly important role at low
twist angles. While the Zeeman contribution is expected to be the same in the
low angle regime, the energetic differences due to the magnetic subbands shown
in our Fig.2(e,f) will change. It would be interesting to explore the
possibility of density induced switching of the Chern number at different
angles or different homobilayers using the framework developed here.

_Acknowledgement—_ We thank B. Andrei Bernevig, Fengcheng Wu, Di Xiao, Jian
Kang, Eslam Khalaf, Jiabin Yu, Daniel Parker, Kin Fai Mak, Jie Shan, and
Xiaodong Xu for helpful discussions. M.W. sincerely thanks Taige Wang for
fruitful discussions and ongoing collaborations. X.W. acknowledges financial
support from the National High Magnetic Field Laboratory through NSF Grant No.
DMR-2128556 and the State of Florida. O.V. was funded by the Gordon and Betty
Moore Foundation’s EPiQS Initiative Grant GBMF11070. Part of the numerical
computations were carried out on resources provided by the Planck clusters at
Florida State University. M.W. acknowledge the Beijing Paratera Co., Ltd. for
providing HPC resources that have contributed to the research results reported
within this paper.

## References

  * Balents  _et al._ [2020] L. Balents, C. R. Dean, D. K. Efetov, and A. F. Young, Superconductivity and strong correlations in moiré flat bands, Nature Physics 16, 725 (2020).
  * Kennes  _et al._ [2021] D. M. Kennes, M. Claassen, L. Xian, A. Georges, A. J. Millis, J. Hone, C. R. Dean, D. Basov, A. N. Pasupathy, and A. Rubio, Moiré heterostructures as a condensed-matter quantum simulator, Nature Physics 17, 155 (2021).
  * Wu  _et al._ [2019] F. Wu, T. Lovorn, E. Tutuc, I. Martin, and A. MacDonald, Topological insulators in twisted transition metal dichalcogenide homobilayers, Physical review letters 122, 086402 (2019).
  * Devakul  _et al._ [2021] T. Devakul, V. Crépel, Y. Zhang, and L. Fu, Magic in twisted transition metal dichalcogenide bilayers, Nature communications 12, 6730 (2021).
  * Li  _et al._ [2021] T. Li, S. Jiang, B. Shen, Y. Zhang, L. Li, Z. Tao, T. Devakul, K. Watanabe, T. Taniguchi, L. Fu, _et al._ , Quantum anomalous hall effect from intertwined moiré bands, Nature 600, 641 (2021).
  * Xie  _et al._ [2022] Y.-M. Xie, C.-P. Zhang, J.-X. Hu, K. F. Mak, and K. T. Law, Valley-polarized quantum anomalous hall state in moiré mote 2/wse 2 heterobilayers, Physical Review Letters 128, 026402 (2022).
  * Cai  _et al._ [2023] J. Cai, E. Anderson, C. Wang, X. Zhang, X. Liu, W. Holtzmann, Y. Zhang, F. Fan, T. Taniguchi, K. Watanabe, _et al._ , Signatures of fractional quantum anomalous hall states in twisted mote2, Nature 622, 63 (2023).
  * Zeng  _et al._ [2023] Y. Zeng, Z. Xia, K. Kang, J. Zhu, P. Knüppel, C. Vaswani, K. Watanabe, T. Taniguchi, K. F. Mak, and J. Shan, Thermodynamic evidence of fractional chern insulator in moiré mote2, Nature 622, 69 (2023).
  * Park  _et al._ [2023] H. Park, J. Cai, E. Anderson, Y. Zhang, J. Zhu, X. Liu, C. Wang, W. Holtzmann, C. Hu, Z. Liu, _et al._ , Observation of fractionally quantized anomalous hall effect, Nature 622, 74 (2023).
  * Xu  _et al._ [2023] F. Xu, Z. Sun, T. Jia, C. Liu, C. Xu, C. Li, Y. Gu, K. Watanabe, T. Taniguchi, B. Tong, _et al._ , Observation of integer and fractional quantum anomalous hall effects in twisted bilayer mote 2, Physical Review X 13, 031037 (2023).
  * Neupert  _et al._ [2011] T. Neupert, L. Santos, C. Chamon, and C. Mudry, Fractional quantum hall states at zero magnetic field, Phys. Rev. Lett. 106, 236804 (2011).
  * Sheng  _et al._ [2011] D. Sheng, Z.-C. Gu, K. Sun, and L. Sheng, Fractional quantum hall effect in the absence of landau levels, Nature Communications 2, 10.1038/ncomms1380 (2011).
  * Regnault and Bernevig [2011] N. Regnault and B. A. Bernevig, Fractional chern insulator, Phys. Rev. X 1, 021014 (2011).
  * Tang  _et al._ [2011] E. Tang, J.-W. Mei, and X.-G. Wen, High-temperature fractional quantum hall states, Phys. Rev. Lett. 106, 236802 (2011).
  * Sun  _et al._ [2011] K. Sun, Z. Gu, H. Katsura, and S. Das Sarma, Nearly flatbands with nontrivial topology, Phys. Rev. Lett. 106, 236803 (2011).
  * Abouelkomsan  _et al._ [2023] A. Abouelkomsan, E. J. Bergholtz, and S. Chatterjee, Multiferroicity and topology in twisted transition metal dichalcogenides (2023), arXiv:2210.14918 [cond-mat.str-el] .
  * Wang  _et al._ [2023a] C. Wang, X.-W. Zhang, X. Liu, Y. He, X. Xu, Y. Ran, T. Cao, and D. Xiao, Fractional chern insulator in twisted bilayer mote $\\_2$, arXiv preprint arXiv:2304.11864 (2023a).
  * Crépel and Fu [2023] V. Crépel and L. Fu, Anomalous hall metal and fractional chern insulator in twisted transition metal dichalcogenides, Physical Review B 107, L201109 (2023).
  * Yu  _et al._ [2024] J. Yu, J. Herzog-Arbeitman, M. Wang, O. Vafek, B. A. Bernevig, and N. Regnault, Fractional chern insulators versus nonmagnetic states in twisted bilayer mote 2, Physical Review B 109, 045147 (2024).
  * Goldman  _et al._ [2023] H. Goldman, A. P. Reddy, N. Paul, and L. Fu, Zero-field composite fermi liquid in twisted semiconductor bilayers, Physical Review Letters 131, 136501 (2023).
  * Qiu  _et al._ [2023] W.-X. Qiu, B. Li, X.-J. Luo, and F. Wu, Interaction-driven topological phase diagram of twisted bilayer mote 2, Physical Review X 13, 041026 (2023).
  * Wang  _et al._ [2023b] T. Wang, T. Devakul, M. P. Zaletel, and L. Fu, Topological magnets and magnons in twisted bilayer mote $\\_2$ and wse $\\_2$, arXiv preprint arXiv:2306.02501 (2023b).
  * Wang  _et al._ [2023c] T. Wang, M. Wang, W. Kim, S. G. Louie, L. Fu, and M. P. Zaletel, Topology, magnetism and charge order in twisted mote2 at higher integer hole fillings (2023c), arXiv:2312.12531 [cond-mat.str-el] .
  * Fan  _et al._ [2024] F.-R. Fan, C. Xiao, and W. Yao, Orbital chern insulator at $\nu$=- 2 in twisted mote 2, Physical Review B 109, L041403 (2024).
  * Liu  _et al._ [2023] X. Liu, C. Wang, X.-W. Zhang, T. Cao, and D. Xiao, Gate-tunable antiferromagnetic chern insulator in twisted bilayer transition metal dichalcogenides, arXiv preprint arXiv:2308.07488 (2023).
  * Luo  _et al._ [2024] X.-J. Luo, W.-X. Qiu, and F. Wu, Majorana zero modes in twisted transition metal dichalcogenide homobilayers, Physical Review B 109, L041103 (2024).
  * Li  _et al._ [2024] B. Li, W.-X. Qiu, and F. Wu, Electrically tuned topology and magnetism in twisted bilayer mote 2 at $\nu$ h= 1, Physical Review B 109, L041106 (2024).
  * Lu  _et al._ [2021] X. Lu, B. Lian, G. Chaudhary, B. A. Piot, G. Romagnoli, K. Watanabe, T. Taniguchi, M. Poggio, A. H. MacDonald, B. A. Bernevig, _et al._ , Multiple flat bands and topological hofstadter butterfly in twisted bilayer graphene close to the second magic angle, Proceedings of the National Academy of Sciences 118, e2100006118 (2021).
  * Foutty  _et al._ [2023] B. A. Foutty, C. R. Kometter, T. Devakul, A. P. Reddy, K. Watanabe, T. Taniguchi, L. Fu, and B. E. Feldman, Mapping twist-tuned multi-band topology in bilayer wse $\\_2$, arXiv preprint arXiv:2304.09808 (2023).
  * Streda [1982] P. Streda, Theory of quantised hall conductivity in two dimensions, Journal of Physics C: Solid State Physics 15, L717 (1982).
  * Widom [1982] A. Widom, Thermodynamic derivation of the hall effect current, Physics Letters A 90, 474 (1982).
  * Kohmoto [1985] M. Kohmoto, Topological invariant and the quantization of the hall conductance, Annals of Physics 160, 343 (1985).
  * Wannier [1978] G. H. Wannier, A result not dependent on rationality for bloch electrons in a magnetic field, physica status solidi (b) 88, 757 (1978).
  * Thouless [1983] D. J. Thouless, Quantization of particle transport, Phys. Rev. B 27, 6083 (1983).
  * MacDonald [1983] A. H. MacDonald, Landau-level subband structure of electrons on a square lattice, Phys. Rev. B 28, 6713 (1983).
  * Tešanović  _et al._ [1989] Z. Tešanović, F. Axel, and B. I. Halperin, “hall” crystal versus wigner crystal, Phys. Rev. B 39, 8525 (1989).
  * Spanton  _et al._ [2018] E. M. Spanton, A. A. Zibrov, H. Zhou, T. Taniguchi, K. Watanabe, M. P. Zaletel, and A. F. Young, Observation of fractional chern insulators in a van der waals heterostructure, Science 360, 62 (2018).
  * Xie  _et al._ [2021] Y. Xie, A. T. Pierce, J. M. Park, D. E. Parker, E. Khalaf, P. Ledwith, Y. Cao, S. H. Lee, S. Chen, P. R. Forrester, K. Watanabe, T. Taniguchi, A. Vishwanath, P. Jarillo-Herrero, and A. Yacoby, Fractional chern insulators in magic-angle twisted bilayer graphene, Nature 600, 439–443 (2021).
  * Wang and Vafek [2023a] X. Wang and O. Vafek, Theory of correlated chern insulators in twisted bilayer graphene (2023a), arXiv:2310.15982 [cond-mat.mes-hall] .
  * Xiao  _et al._ [2012a] D. Xiao, G.-B. Liu, W. Feng, X. Xu, and W. Yao, Coupled spin and valley physics in monolayers of ${\mathrm{mos}}_{2}$ and other group-vi dichalcogenides, Phys. Rev. Lett. 108, 196802 (2012a).
  * Pan  _et al._ [2020] H. Pan, F. Wu, and S. D. Sarma, Band topology, hubbard model, heisenberg model, and dzyaloshinskii-moriya interaction in twisted bilayer wse 2, Physical Review Research 2, 033087 (2020).
  * Polshyn  _et al._ [2020] H. Polshyn, J. Zhu, M. A. Kumar, Y. Zhang, F. Yang, C. L. Tschirhart, M. Serlin, K. Watanabe, T. Taniguchi, A. H. MacDonald, and A. F. Young, Electrical switching of magnetic order in an orbital chern insulator, Nature 588, 66–70 (2020).
  * Wang and Vafek [2022] X. Wang and O. Vafek, Narrow bands in magnetic field and strong-coupling hofstadter spectra, Physical Review B 106, L121111 (2022).
  * Wang and Vafek [2023b] X. Wang and O. Vafek, Revisiting bloch electrons in magnetic field: Hofstadter physics via hybrid wannier states, arXiv preprint arXiv:2303.16347 (2023b).
  * [45] See supplemental material for orbital contribution to g factor, landau theory of transition between chern paraelectric and ferroelectirc, as well as hybrid wannier basis at finite magnetic field.
  * Xiao  _et al._ [2012b] D. Xiao, G.-B. Liu, W. Feng, X. Xu, and W. Yao, Coupled spin and valley physics in monolayers of mos 2 and other group-vi dichalcogenides, Physical review letters 108, 196802 (2012b).
  * Jia  _et al._ [2023] Y. Jia, J. Yu, J. Liu, J. Herzog-Arbeitman, Z. Qi, N. Regnault, H. Weng, B. A. Bernevig, and Q. Wu, Moir$\backslash$’e fractional chern insulators i: First-principles calculations and continuum models of twisted bilayer mote $\\_2$, arXiv preprint arXiv:2311.04958 (2023).
  * Tešanović  _et al._ [1989] Z. Tešanović, F. Axel, and B. Halperin, “hall crystal”versus wigner crystal, Physical Review B 39, 8525 (1989).
  * Zhang  _et al._ [2023] X.-W. Zhang, C. Wang, X. Liu, Y. Fan, T. Cao, and D. Xiao, Polarization-driven band topology evolution in twisted mote $\\_2$ and wse $\\_2$, arXiv preprint arXiv:2311.12776 (2023).
  * Mao  _et al._ [2023] N. Mao, C. Xu, J. Li, T. Bao, P. Liu, Y. Xu, C. Felser, L. Fu, and Y. Zhang, Lattice relaxation, electronic structure and continuum model for twisted bilayer mote $\\_2$, arXiv preprint arXiv:2311.07533 (2023).
  * Kang and Vafek [2020] J. Kang and O. Vafek, Non-abelian dirac node braiding and near-degeneracy of correlated phases at odd integer filling in magic-angle twisted bilayer graphene, Physical Review B 102, 035161 (2020).

Supplemental Material for “Interacting phase diagram of twisted MoTe2 in
magnetic field”

###### Contents

  1. I Moiré geometry
  2. II Orbital contribution to $g$ factor
  3. III Landau Theory of transition between Chern paraelectric and QH ferroelectric
  4. IV Hybrid Wannier states at $B=0$
  5. V Constructing hWS basis at $B\neq 0$
  6. VI $\hat{t}(\mathbf{a}_{2})$ symmetry
  7. VII Extended interaction diagrams

##  I Moiré geometry

Moiré primitive and reciprocal primitive lattice vectors are shown in Fig. 1,
and given explicitly as follows:

$\mathbf{a}_{1}=a_{M}\left(\frac{\sqrt{3}}{2},-\frac{1}{2}\right),\mathbf{a}_{2}=a_{M}(0,1),\
\textbf{G}_{1}=\frac{4\pi}{\sqrt{3}a_{M}}(1,0),\textbf{G}_{2}=\frac{4\pi}{\sqrt{3}a_{M}}\left(\frac{1}{2},\frac{\sqrt{3}}{2}\right),$
(S1)

where the moiré period $a_{M}={a_{0}}/{\theta}$. At $B=0$, we take the rhombus
Brillouin zone with orange boundary lines in Fig. 1, with wavevector defined
as $\mathbf{k}=k_{1}\mathbf{G}_{1}+k_{2}\mathbf{G}_{2}$ where
$k_{i=1,2}\in[0,1)$. At $B\neq 0$, we use the Landau gauge $\mathbf{A}=(0,Bx)$
and work at rational magnetic flux per unit cell ${\phi}/{\phi_{0}}={p}/{q}$,
where $\\{p,q\\}$ are coprime integers and $\phi_{0}=h/e$ is the flux quantum.
Hence, magnetic translation operators are

$\hat{t}(\mathbf{a}_{1})=e^{-i\mathbf{q}_{\phi}\cdot\mathbf{r}}\hat{T}(\mathbf{a}_{1})=e^{-i\frac{2\pi}{a_{M}}\frac{p}{q}\textbf{y}}\hat{T}(\mathbf{a}_{1}),\hat{t}(\mathbf{a}_{2})=\hat{T}(\mathbf{a}_{2})$
(S2)

with $\hat{T}(\mathbf{a}_{\\{1,2\\}})$ translation operators at $B=0$(see, for
example, Ref.[44]). Since
$[\hat{t}(\mathbf{a}_{1}),\hat{t}^{q}(\mathbf{a}_{2})]=0$, we can construct
magnetic moiré Brillouin zone as
$\\{k_{1}\mathbf{G}_{1}+k_{2}\mathbf{G}_{2}|0<k_{1}<1,0<k_{2}<\frac{1}{q}\\}$,
whose boundary lines are shown as purple at $\phi/\phi_{0}=1/11$. At
$\phi/\phi_{0}\ll 1$, we find it sufficient to perform Hartree-Fock
calculations with a discretized mesh given by the cyan and brown points. Note
if the Hartree-Fock ground state preserves $\hat{t}(\mathbf{a}_{2})$, the
calculations can be further simplified by constructing the Hartree-Fock
Hamiltonian at the cyan point only[44, 43]. When constructing hybrid Wannier
states at $B=0$, we consider real space as $N_{1}\mathbf{a}_{1}\times
N_{2}\mathbf{a}_{2}$, which makes $\frac{A_{uc}}{A}=\frac{1}{N_{1}N_{2}}$ and
discrete $\mathbf{k}$ points
$\mathbf{k}=\frac{n_{1}}{N_{1}}\mathbf{G}_{1}+\frac{n_{2}}{N_{2}}\mathbf{G}_{2}=k_{1}\mathbf{G}_{1}+k_{2}\mathbf{G}_{2}$
with integers $n_{1}$ and $n_{2}$.

Figure 1: Schematics of moiré Brillouin zone at $B=0$ and $B\neq 0$. Sample
points are also shown in the figure. Inset plot shows probability density of
hybrid Wannier states in tMoTe2.

##  II Orbital contribution to $g$ factor

Correct band basis and fitting parameters in MoTe2 are crucial to estimate the
orbital contribution to Zeeman effect. The Ref. [46] is based on monolayer
calculation and Ref. [3] on aligned bilayer; their band basis of MoTe2 near
$\mathbf{K}$ and $\mathbf{K}^{\prime}$ valley are the same. Namely, conduction
bands in both valleys are described by $d_{z^{2}}$ orbital, valence band in
$\mathbf{K}$ is by $d_{x^{2}-y^{2}}+id_{xy}$ and valence band in
$\mathbf{K}^{\prime}$ by $d_{x^{2}-y^{2}}-id_{xy}$. The effective
$\mathbf{k}\cdot\mathbf{p}$ models (for one layer) differ only by $\sigma_{0}$
term in valence-conduction $2\times 2$ subspace, which simply shifts total
energy and will not affect effective mass in each valley. We therefore choose
to start with massive Dirac Hamiltonian in $\mathbf{K}$ valley[3]

$h_{l}^{\mathbf{K}}(\mathbf{k},d_{0})=\begin{pmatrix}\Delta_{g}+\Delta_{l,c}(d_{0})&\hbar
v_{F}(k_{x}-ik_{y})\\\ \hbar
v_{F}(k_{x}+ik_{y})&\Delta_{l,v}(d_{0})\end{pmatrix}$ (S3)

here $c$ denotes conduction bands above semiconductor gap $\Delta_{g}$, and
$v$ denotes valence bands below the gap. Since $\Delta_{g}$($\approx 2eV$) is
much larger than other energy scales, we could write eigen-equations for
valence band wavefuction $\ket{\psi_{v}}$ as

$\begin{split}&\hbar v_{F}(k_{x}+ik_{y})\frac{-\hbar
v_{F}(k_{x}-ik_{y})}{\Delta_{g}+\Delta_{l,c}(d_{0})-E}\ket{\psi_{v}}+\Delta_{l,v}(d_{0})\ket{\psi_{v}}\\\
&\approx\frac{-\hbar^{2}v_{F}^{2}}{\Delta_{g}}(k_{x}+ik_{y})(k_{x}-ik_{y})\ket{\psi_{v}}+\Delta_{l,v}(d_{0})\ket{\psi_{v}}\\\
&=(-\frac{p^{2}}{2m^{*}}+\Delta_{l,v}(d_{0}))\ket{\psi_{v}}=E\ket{\psi_{v}},\\\
\end{split}$ (S4)

where $m^{*}=\frac{\Delta_{g}}{2v_{F}^{2}}\approx 0.62m_{e}$ is the effective
mass for valence band electron [3, 46]. Within this approximation, the Dirac
equation gives an approximately parabolic kinetic energy for valence band
electrons.

Now consider adding magnetic field through minimal coupling. Note the original
form of kinetic energy is $\frac{-1}{2m^{*}}(p_{x}+ip_{y})(p_{x}-ip_{y})$ for
$\mathbf{K}$ valley electron. Minimal coupling within the original form is not
equal to the minimal coupling within $\frac{-p^{2}}{2m^{*}}$, because the
momentum operator does not commute with the vector potential $\mathbf{A}$,

$\begin{split}-\frac{1}{2m^{*}}(p_{x}+ip_{y})(p_{x}-ip_{y})&\rightarrow-\frac{1}{2m^{*}}(p_{x}+eA_{x}+ip_{y}+ieA_{y})(p_{x}+eA_{x}-ip_{y}-ieA_{y})\\\
&=-\frac{1}{2m^{*}}((p_{x}+eA_{x})^{2}+(p_{y}+eA_{y})^{2}+i[p_{y}+eA_{y},p_{x}+eA_{x}])\\\
&=-\frac{1}{2m^{*}}((p_{x}+eA_{x})^{2}+(p_{y}+eA_{y})^{2}+ie([p_{y},A_{x}]+[A_{y},p_{x}]))\\\
&=-\frac{1}{2m^{*}}((p_{x}+eA_{x})^{2}+(p_{y}+eA_{y})^{2})-\frac{1}{2m^{*}}ie(-i\hbar\partial_{y}A_{x}+i\hbar\partial_{x}A_{y})\\\
&=-\frac{1}{2m^{*}}((p_{x}+eA_{x})^{2}+(p_{y}+eA_{y})^{2})-\frac{1}{2m^{*}}ie(+i\hbar)(\nabla\times
A)\\\
&=-\frac{1}{2m^{*}}((p_{x}+eA_{x})^{2}+(p_{y}+eA_{y})^{2})+\frac{e}{2m^{*}}\hbar
B.\end{split}$ (S5)

The first term in the last line is already included by replacing
$\mathbf{p}\rightarrow\mathbf{p}+e\mathbf{A}$ in the continuum model at $B=0$
presented in the main text, and the second term can be understood as the
orbital contribution to the Zeeman energy shift. Indeed, comparing to the
definition of the $g$ factor $H_{Zeeman}=-\mu\cdot
B=\frac{e}{2m_{e}}g\frac{\hbar}{2}B$ we get

$g_{orbit}=2\frac{m_{e}}{m^{*}}=2\times\frac{1}{0.62}\approx 3.23.$ (S6)

The massive Dirac Hamiltonian in the $\mathbf{K}^{\prime}$ valley is different
from Eq. S3, since the valence band basis changes[46]:

$h_{l}^{\mathbf{K}^{\prime}}(\mathbf{k},d_{0})=\begin{pmatrix}\Delta_{g}+\Delta_{l,c}(d_{0})&\hbar
v_{F}(-k_{x}-ik_{y})\\\ \hbar
v_{F}(-k_{x}+ik_{y})&\Delta_{l,v}(d_{0})\end{pmatrix}.$ (S7)

Therefore for $\mathbf{K}^{\prime}$ valley electrons, kinetic energy reads
$\frac{-1}{2m^{*}}(p_{x}-ip_{y})(p_{x}+ip_{y})$, and hence additional energy
shift is $-\frac{e}{2m^{*}}\hbar B$. The opposite sign of this orbital
contribution is again aligned with the opposite spin orientation in the valley
$\mathbf{K}^{\prime}$ due to spin-valley locking, making the total $g$ factor
in both valleys approximately $5.23$.

##  III Landau Theory of transition between Chern paraelectric and QH
ferroelectric

In order to further explore the first order phase transition and the related
crossover physics, we write a phenomenological Landau theory and show that it
can capture the results of Hartree Fock calculation presented in the main
text. At $u_{D}=0$ the Hamiltonian is invariant under $C_{2y}\mathcal{T}$
which is spontaneously broken in the quantum Hall ferroelectric phase and
preserved in the Chern paraelectric phase. The layer polarization $P_{l}$,
which changes sign under $C_{2y}\mathcal{T}$, can therefore serve as a Landau
order parameter. The Landau free energy must be invariant under
$C_{2y}\mathcal{T}$ and therefore can only contain even powers of $P_{l}$ at
$u_{D}=0$. In order to obtain the first order transition at $\epsilon_{*}$,
one needs a negative prefactor of the $P_{l}^{4}$ term (and of course a
positive prefactor of the $P_{l}^{6}$ term in order to stabilize the free
energy). Non-zero $u_{D}$ then couples directly to $P_{l}$ at linear order.
Therefore,

$\begin{gathered}\delta E=-\gamma
u_{D}P_{l}+\frac{1}{2}a(\epsilon)P_{l}^{2}-\frac{1}{4}bP_{l}^{4}+\frac{1}{6}P_{l}^{6}.\end{gathered}$
(S8)

Here the prefactor of $P_{l}^{6}$ term is set to $1/6$ for convenience which
can always be achieved by rescaling $\delta E$. Constant $\gamma$ is added to
set the minimum position of $\delta E$ at $0<P_{l}<1$. Near the phase
boundary, the prefactor of $P_{l}^{2}$ term is treated as a function of
dielectric constant $\epsilon$, while the prefactor $b$ is approximated as a
constant. Note $b<0$ will turn transition into a continuous one, hence is not
relavant to the phase transition studied in the main text.

At $u_{D}=0$ the free energy has a minimum at $P_{l}=0$ and, depending on the
value of $a(\epsilon)$, it may have another pair of equal and opposite minima
at non-zero $P_{l}$. At $a(\epsilon)=3b^{2}/16$ the minima are degenerate and
this corresponds to the first order phase transition at $u_{D}=0$.

At small non-zero $u_{D}$ one of the minima shifts away from $P_{l}=0$, the
degeneracy between the other two minima is lifted, and the transition remains
first order. At very large $u_{D}$ there is a single minimum. To find the
critical endpoint, we analyze the extrema of $\delta E$:

$\begin{split}\frac{d\delta E}{dP_{l}}&=-\gamma
u_{D}+a(\epsilon)P_{l}-bP_{l}^{3}+P_{l}^{5}\\\ &=-\gamma
u_{D}+P_{l}(P_{l}^{4}-bP_{l}^{2}+a(\epsilon)).\\\ \end{split}$ (S9)

To find the critical endpoint, we find the smallest value of $a(\epsilon)$ for
which $P_{l}(P_{l}^{4}-bP_{l}^{2}+a(\epsilon))$ is strictly increasing for any
$P_{l}$. We have

$\begin{split}\frac{d}{dP_{l}}(P_{l}(P_{l}^{4}-bP_{l}^{2}+a(\epsilon)))&=5P_{l}^{4}-3bP_{l}^{2}+a(\epsilon)\\\
&=5(P_{l}^{2}-\frac{3b}{10})^{2}+a(\epsilon)-\frac{9b^{2}}{20}\geq 0.\\\
\end{split}$ (S10)

The smallest $a(\epsilon)$ for which the above is always positive gives the
critical $a(\epsilon)=9b^{2}/20$, at which the above derivative vanishes if we
set $P_{l}$ to $\sqrt{3b/10}$. The value of $\gamma u_{D}$ at the critical
endpoint is then calculated by requiring the derivative of $\delta E$ to
vanish at this $a(\epsilon)$ and $P_{l}$:

$\begin{split}\gamma u_{D}&=P_{l}(P_{l}^{4}-bP_{l}^{2}+a(\epsilon))\\\
&=\sqrt{\frac{3}{10}b}(\frac{9}{100}b^{2}-b\frac{3b}{10}+\frac{9}{20}b^{2})=\frac{3\sqrt{30}}{125}b^{\frac{5}{2}}.\\\
\end{split}$ (S11)

Figure 2:  Order parameter $P_{l}$ with varying $a(\epsilon)$ and $\gamma
u_{D}$ in Landau theory. These two plots capture the results obtained by
finite magnetic field Hartree Fock in the main text. The parameter $b$ is
fixed at $b=0.7$ to make $P_{l}\in[0,1)$. Analytic expression of critical end
point(red star) is $(a(\epsilon),\gamma
u_{D})=(9b^{2}/20,3\sqrt{30}b^{5/2}/125)$. At $\gamma u_{D}=0$, the first
order transition happens at $a(\epsilon)=3b^{2}/16$.

##  IV Hybrid Wannier states at $B=0$

Hybrid Wannier states at $B=0$ are eigenstates of operator
$\hat{O}=\hat{P}e^{-i\frac{1}{N_{1}}\mathbf{G}_{1}\cdot\mathbf{r}}\hat{P}$[51,
44]. Here
$\hat{P}=\sum_{n,\mathbf{k}}\ket{\psi_{n,\mathbf{k}}}\bra{\psi_{n,\mathbf{k}}}$
is the projector onto Bloch bands of interest. Since
$[\hat{O},\hat{T}(\mathbf{a}_{2})]=0$, $k_{2}$ remains as a good quantum
number, $\hat{T}(\mathbf{a}_{2})\ket{w(k_{2})}=e^{-i2\pi
k_{2}}\ket{w(k_{2})}$. Suppose $\ket{w(k_{2})}$ is an eigenstate of $\hat{O}$,
then $\hat{T}(\mathbf{a}_{1})\ket{w(k_{2})}$ is also an eigenstate, with
eigenvalue shifted:

$\begin{gathered}\hat{O}\ket{w(k_{2})}=e^{-i\frac{1}{N_{1}}\langle\mathbf{G}_{1}\cdot\mathbf{r}\rangle}\ket{w(k_{2})}\\\
\hat{O}\hat{T}(\mathbf{a}_{1})\ket{w(k_{2})}=\hat{T}(\mathbf{a}_{1})\hat{O}e^{-i\frac{1}{N_{1}}\mathbf{G}_{1}\cdot\mathbf{a}_{1}}\ket{w(k_{2})}=e^{-i\frac{1}{N_{1}}\mathbf{G}_{1}\cdot\mathbf{a}_{1}}e^{-i\frac{1}{N_{1}}\langle\mathbf{G}_{1}\cdot\mathbf{r}\rangle}\hat{T}(\mathbf{a}_{1})\ket{w(k_{2})}.\end{gathered}$
(S12)

Therefore we can denote all eigenstates as
$\ket{w_{\alpha}(n_{0},k_{2})}=\hat{T}^{n_{0}}(\mathbf{a}_{1})\ket{w_{\alpha}(0,k_{2})}$,
where $\ket{w_{\alpha}(0,k_{2})}$ are all eigenstates whose eigenvalues
satisfy
$-\pi\leq\langle\mathbf{G}_{1}\cdot\mathbf{r}\rangle<\pi=\frac{\mathbf{G}_{1}\cdot\mathbf{a}_{1}}{2}$.
We show the evolution of $\langle\mathbf{G}_{1}\cdot\mathbf{r}\rangle$ with
varying $k_{2}$, for
$\hat{P}_{1}=\sum_{\mathbf{k}}\ket{\psi_{n=1,\mathbf{k}}}\bra{\psi_{n=1,\mathbf{k}}}$
and
$\hat{P}_{2}=\sum_{\mathbf{k}}\ket{\psi_{n=2,\mathbf{k}}}\bra{\psi_{n=2,\mathbf{k}}}$
in Fig. 3(a); the top most valance band ($\hat{P}_{1}$, blue) and the band
below it ($\hat{P}_{2}$, orange) at $\theta=3.89^{\circ}$ in $\mathbf{K}$
valley of tMoTe2. The winding correponds to their Chern numbers $+1$ and $-1$
respectively. Then we show the evolution for
$\hat{P}_{3}=\sum_{n=\\{1,2\\},\mathbf{k}}\ket{\psi_{n,\mathbf{k}}}\bra{\psi_{n,\mathbf{k}}}$
in Fig. 3(b). Corresponding plots for real space probability density in each
layer of $\ket{w_{\alpha}(n_{0}=0,k_{2}=0)}$ is shown. All hybrid Wannier
states $\ket{w_{\alpha}(n_{0},k_{2})}$ are localized along $\mathbf{a}_{1}$
but extended along $\mathbf{a}_{2}$. In addition, hybrid Wannier state
$\ket{w_{\alpha}(n_{0},k_{2})}$ can be expressed in terms of Bloch states
$\psi_{n,\mathbf{k}}$ via a unitary transformation $U^{hWS}$ at each
$\mathbf{k}$ point

$\ket{w_{\alpha}(n_{0},k_{2})}=\frac{1}{\sqrt{N_{1}}}\sum_{n={1,2}}\sum_{k_{1}}e^{-i2\pi
k_{1}n_{0}}U_{n,\alpha}^{hWS}(\mathbf{k})\ket{\psi_{n,\mathbf{k}}}.$ (S13)

Figure 3: (a, b)Average of
$\langle\mathbf{\mathbf{G}}_{1}\cdot\mathbf{r}\rangle$ of hybrid Wannier
states. (a) blue(orange) points denote the average value for hWS constructed
from the first(second) band of tMoTe2 in $\mathbf{K}$ valley. (b)average value
for two hWSs constructed from first two bands. (c) up:Probability density of
$\ket{w(n_{0}=0,k_{2}=0)}$ from first band in $\mathbf{K}$ valley.
down:Probability density of $\ket{w(n_{0}=0,k_{2}=0)}$ from second band in
$\mathbf{K}$ valley. (d)Probability density of $\ket{w_{1}(n_{0}=0,k_{2}=0)}$
and $\ket{w_{2}(n_{0}=0,k_{2}=0)}$ from first two bands in $\mathbf{K}$
valley.

##  V Constructing hWS basis at $B\neq 0$

In order to perform the Hartree Fock calculations at $B\neq 0$ we construct
the projector onto the subspace spanned by the eigenstates of the non-
interacting Hamiltonian with minimal coupling to $B$. Specificly, we are
interested in all eigenstates at $B\neq 0$ that evolve from the first two
valence bands in each valley at $B=0$. The exact projector can be obtained
from complete Landau level(LL) wavefunctions[44]. This work adapts hybrid
Wannier states at $B=0$ to approximate this projector at low magnetic field,
with excellent accuracy shown in the Fig. 4(see Ref. [44] and Ref. [43] for
further discussion). This method allows us to perform the Hartree Fock
calculations at low $B$ more efficiently than the direct Landau level
expansion. In the Landau gauge $\mathbf{A}=(0,Bx)$ and low $B$, hybrid Wannier
states $\ket{w_{\alpha}(n_{0}=0,k_{2})}$ (which are constructed at $B=0$) have
an approximately same expectation value of the Hamiltonian (or any polynomial
function of it) because they only extend in $\hat{y}$ and by being localized
in the $\hat{x}$ near origin the region of space where the vector potential
$\mathbf{A}$ is large never contributes. Therefore,
$\ket{w_{\alpha}(n_{0}=0,k_{2})}$ must be almost completely within the $B\neq
0$ Hilbert subspace evolving from the narrow bands of interest. To generate
the rest of the basis at $B\neq 0$, we can employ the magnetic translation
group at rational magnetic flux ratios $\phi/\phi_{0}=p/q$. Eigenstates of the
magnetic translations $\hat{t}(\mathbf{a}_{1})$ and $\hat{t}(\mathbf{a}_{2})$
can therefore be obtained as

$\ket{W_{\alpha,r}(\mathbf{k})}=\frac{1}{\sqrt{N_{1}}}\sum_{s=\lfloor-\frac{N_{1}}{2}\rfloor}^{\lfloor\frac{N_{1}}{2}\rfloor}e^{i2\pi
k_{1}s}\hat{t}^{s}(\mathbf{a}_{1})\ket{w_{\alpha}(0,k_{2}+\frac{r}{q})}.$
(S14)

Here $\ket{w_{\alpha=\\{1,2\\}}(0,k_{2}+\frac{r}{q})}$ are constructed from
$\hat{P}_{3}$ which projects onto the band composite with zero total Chern
number. Here $k_{1}\in[0,1)$. $k_{2}\in[0,\frac{1}{q})$ and $r=0,1,...,q-1$ is
additional index. The total number of states at zero and finite B are the
same, as must be the case because the total Chern number vanishes. We note in
passing that for projectors $\hat{P}_{1}$ and $\hat{P}_{2}$ whose total Chern
numbers are nonzero, we need an additional procedure to construct proper basis
at $B\neq 0$ [44, 43]. Under the action of $\hat{t}(\mathbf{a}_{1})$,
$\hat{t}^{q}(\mathbf{a}_{2})$ and $\hat{t}(\mathbf{a}_{2})$, the basis
$\ket{W_{\alpha,r}(\mathbf{k})}$ satisfies

$\begin{gathered}\hat{t}(\mathbf{a}_{1})\ket{W_{\alpha,r}(\mathbf{k})}=e^{-i2\pi
k_{1}}\ket{W_{\alpha,r}(\mathbf{k})}\\\
\hat{t}^{q}(\mathbf{a}_{2})\ket{W_{\alpha,r}(\mathbf{k})}=e^{-i2\pi
k_{2}q}\ket{W_{\alpha,r}(\mathbf{k})}\\\
\hat{t}(\mathbf{a}_{2})\ket{W_{\alpha,r}(\mathbf{k})}=e^{-i2\pi(k_{2}+\frac{r}{q})}\ket{W_{\alpha,r}(\mathbf{k}+\frac{p}{q}\mathbf{G}_{1})}.\\\
\end{gathered}$ (S15)

The last line follows from
$\hat{t}(\mathbf{a}_{2})\hat{t}(\mathbf{a}_{1})=e^{i2\pi\frac{\phi}{\phi_{0}}}\hat{t}(\mathbf{a}_{1})\hat{t}(\mathbf{a}_{2})$.

For $\frac{\phi}{\phi_{0}}=\frac{p}{q}$, we impose a requirement on the
discrete k-space mesh at $B=0$ so that $N_{1}$ is evenly divided by $2q$ and
$N_{2}$ is evenly divided by $q$. This originates from representation of
$\mathbf{q}_{\phi}$ in hexagonal (triangular) lattice

$\mathbf{q}_{\phi}=\frac{p}{q}\frac{2\pi}{a_{M}}\hat{y}=\frac{p}{q}\frac{a_{1y}}{a_{M}}\mathbf{G}_{1}+\frac{p}{q}\mathbf{G}_{2}=-\frac{p}{2q}\mathbf{G}_{1}+\frac{p}{q}\mathbf{G}_{2}$
(S16)

consequently, the even division is needed to represent the operator
$e^{-i\mathbf{q}_{\phi}\cdot\mathbf{r}}$. This is crucial to construct hybrid
Wannier basis at $B\neq 0$.

Using $\ket{W_{\alpha,r}(\mathbf{k})}$, the matrix elements of the Hamiltonian
can be easily evaluated[44]. However, it is worth mentioning that
$\ket{W_{\alpha,r}(\mathbf{k})}$ are not guaranteed to be orthonormal
(although they nearly are at low $B$). Therefore, a transformation is needed
to generate an orthonormal basis out of them. This will mix indices
$\\{\alpha,r\\}$ to index $a$, with $dim(a)=dim(\alpha)dim(r)$[44, 43]. The
orthonormal basis which we use is therefore

$\ket{V_{a}(\mathbf{k})}=\sum_{\alpha,r}\ket{W_{\alpha,r}(\mathbf{k})}\left(U\frac{1}{\sqrt{D}}\right)_{\alpha
r,a}\\!\\!\\!\\!\\!\\!\\!\\!(\mathbf{k}),$ (S17)

where $U$ and $D$ are respectively unitary matrix and diagonal matrix of the
overlap matrix $\Lambda_{\alpha,\beta,r_{1},r_{2}}\equiv\langle
W_{\beta,r_{2}}(\mathbf{k})|W_{\alpha,r_{1}}(\mathbf{k})\rangle$, namely
$\Lambda=UDU^{\dagger}$.

Figure 4: Comparison of the Hofstadter spectra using exact LL approach and
approximate hybrid Wannier state(hWS) approach. (a) Non-interacting
Hamiltonian in the $\mathbf{K}$ valley (b) Non-interacting Hamiltonian in the
$\mathbf{K}^{\prime}$ valley (c) Hartree Fock mean field Hamiltonian at
$\frac{\phi}{\phi_{0}}=\frac{1}{12}$. The good agreement demonstrates the
accuracy of the hWS approach.

##  VI $\hat{t}(\mathbf{a}_{2})$ symmetry

During our calculation, we choose points as
$k_{1}=0,\frac{1}{q},...,\frac{q-1}{q}$ and $k_{2}=0$(Fig. 1). We observe that
$\hat{t}(\mathbf{a}_{2})$ symmetry is preserved down to flux
$\phi/\phi_{0}=1/18$, which we interpret to mean that ground states in TMD
along (-1, -1), (-1, 0) and (-1, 1) Streda line do not break the
$\hat{t}(\mathbf{a}_{2})$ symmetry. The similar phenomena is also observed
through $B\neq 0$ Hartree-Fock calculations in twisted bilayer graphene
without strain [39]. Therefore, we extend the Hartree Fock calculation down to
$\phi/\phi_{0}=1/47$ by assuming that $\hat{t}(\mathbf{a}_{2})$ symmetry is
preserved, which allows us to significantly speed up the computation at such
low $B$.

##  VII Extended interaction diagrams

In extended interaction diagrams, we present, among other results, valley
polarization
$P_{\eta}\equiv\frac{A_{uc}}{A}\sum_{l}\int\mathrm{d}^{2}\mathbf{r}\langle\psi_{l\mathbf{K}}^{\dagger}(\mathbf{r})\psi_{l\mathbf{K}}(\mathbf{r})-\psi_{l\mathbf{K}^{\prime}}^{\dagger}(\mathbf{r})\psi_{l\mathbf{K}^{\prime}}(\mathbf{r})\rangle$.
$P_{\eta}=1$ for $\mathbf{K}$ valley polarzied states, and $P_{\eta}=-1$ for
$\mathbf{K}^{\prime}$ valley polarized states.

Figure 5: (a)Interacting diagram at $B=0$ at integer filling $\nu=-1$[23].
Orange line at $\epsilon=20$ intersects the $B=0$ topological transition from
a Chern $\pm 1$ insulator to a trivial insulator with a vanishing Chern number
but with intervalley coherence. The behavior along such line, at the same
$\epsilon$, but at $B=8.78T$ is presented in the main text. Panels (b-c) show
the charge gap and layer polarization at $B=0$, respectively.  Figure 6: IVC-
VP transition at $u_{D}=16.8meV$ with different $g$ factors. Charge gap and
$P_{\eta}$ at (a,b)$g=0$ and at (c,d)$g=2$ are shown. (e) Extrapolated IVC-VP
transition lines at various $g$. Stars denote original data points from finite
magnetic field Hartree Fock calculation.  Figure 7: (a-b)Charge gap and
$P_{\eta}$ for various filling factors and magnetic field. (c-d)Charge gap and
layer polarization $P_{l}$ for filling factor $\nu$ near $\nu=-1$ and various
$u_{D}$.