# Designer Meron Lattice on the Surface of a Topological Insulator

Daniele Guerci  Center for Computational Quantum Physics, Flatiron Institute,
New York, New York 10010, USA  Jie Wang  Center for Computational Quantum
Physics, Flatiron Institute, New York, New York 10010, USA  J. H. Pixley
Department of Physics and Astronomy, Center for Materials Theory, Rutgers
University, Piscataway, New Jersey 08854, USA  Center for Computational
Quantum Physics, Flatiron Institute, New York, New York 10010, USA  Jennifer
Cano  Department of Physics and Astronomy, Stony Brook University, Stony
Brook, New York 11974, USA  Center for Computational Quantum Physics, Flatiron
Institute, New York, New York 10010, USA

###### Abstract

We present a promising route to realize spontaneous magnetic order on the
surface of a 3D topological insulator by applying a superlattice potential.
The superlattice potential lowers the symmetry of the surface states and
creates tunable van Hove singularities, which, when combined with strong spin-
orbit coupling and Coulomb repulsion give rise to a topological meron lattice
spin texture. The periodicity of this designer meron lattice can be tuned by
varying the periodicity of the superlattice potential. We employ Ginzburg-
Landau theory to classify the different magnetic orders and show that the
magnetic transition temperature reaches experimentally accessible values. Our
work introduces a new direction to realize exotic quantum order by engineering
interacting Dirac electrons in a superlattice potential, with promising
applications to spintronics.

##  I Introduction

Three-dimensional topological insulators (TIs) host gapless surface Dirac
cones protected by time reversal symmetry [1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11].
In the presence of strong interactions, time-reversal symmetry can be
spontaneously broken, gapping the Dirac cone. The result is an exotic
magnetically ordered surface exhibiting the quantized anomalous Hall effect
[5, 12, 13, 14, 15, 16, 17, 18]. However, for all measured TIs, the Coulomb
interaction is too weak to induce the magnetically ordered phase [5, 12, 14,
18].

In this manuscript, we show that a superlattice potential enhances correlation
effects and provides a new and experimentally accessible route to realize
spontaneous magnetic order on the surface of a TI as depicted in Fig. 1(a).
The superlattice potential downfolds and strongly renormalizes the low-energy
band structure, creating satellite Dirac cones without opening a gap [19] and
inducing strong van Hove singularities (VHSs) [20]. We show that the
superlattice induced VHSs drive a spin density wave instability that results
in an exotic meron lattice taking place at arbitrarily small values of the
electron-electron interaction. A meron, as shown in Fig. 1(b), is
topologically equivalent to half a skyrmion: magnetic moments in its core
point up or down, but magnetic moments along its boundary are in-plane [21].
While topological spin textures such as skyrmions have been observed in non-
centrosymmetric magnets [22, 23, 24, 25], ultrathin magnetic films [26, 27,
28] and multiferroic insulators [29, 30], a meron lattice has only been
observed recently [31, 32].

The meron lattice we describe on the surface of the TI is stabilized by the
interplay between the superlattice potential and the strong spin-orbit
coupling (SOC) on the TI surface that “locks” the spin to the momentum,
forcing the magnetic moments to rotate in space about an in-plane axis. It
features several novel aspects: (1) it leaves the surface Dirac cone gapless;
(2) each unit cell exhibits two merons with opposite topological charge; and
(3) the meron lattice periodicity is determined by the applied potential,
i.e., it can be chosen by design. These features make it distinct from other
magnetic textures induced on [5, 12, 13, 14, 15, 16, 17, 18] or proximity-
coupled [33, 34, 35, 36, 37] to the surface of a TI.

Figure 1: Schematic of experiment and meron spin texture. (a) The experiment
consists, from top to bottom, of a metallic gate, a patterned dielectric and a
topological insulator. Applying a bias between the metallic gate and the
topological insulator imposes a modulated potential on the surface of the TI.
(b) Magnetization exhibiting half-integer winding characteristic of a meron
spin texture. Merons pinned by the superlattice potential form a meron
lattice.

In this work we focus on the experimental setup shown in Fig. 1(a), which
depicts a superlattice potential imposed on the surface of a TI by gating a
patterned dielectric stacked above the TI surface. This approach was
introduced to realize a superlattice potential on graphene with periodicity
down to 35nm and strength $\sim 50$meV [38, 39]. It offers great tunability:
the periodicity, strength, and symmetry of the potential can be engineered. A
superlattice potential has been studied theoretically to band engineer
topological materials [40, 19, 20, 41, 42, 43, 19, 20, 44, 45, 46, 47, 48,
49]. Strain has also been employed to modify the dispersion of surface states
[50, 51].

##  II Model of an interacting TI surface

We consider interacting electrons on the surface of a TI described by the
Hamiltonian $\hat{H}=\hat{H}_{0}+\hat{H}_{\rm int}$. The non-interacting
Hamiltonian $\hat{H}_{0}=\int
d^{2}\mathbf{r}\,\hat{\Psi}^{\dagger}(\mathbf{r})H_{0}(\mathbf{r})\hat{\Psi}(\mathbf{r})$
describes a spin-momentum locked Dirac cone subject to a superlattice
potential,

$H_{0}(\mathbf{r})=v_{F}\left(-i\bm{\nabla}_{\mathbf{r}}\times\bm{\sigma}\right)_{z}+\sigma_{0}\,w(\mathbf{r}),$
(1)

where $\bm{\nabla}_{\mathbf{r}}=(\partial_{x},\partial_{y})$,
$\bm{\sigma}=(\sigma_{x},\sigma_{y},\sigma_{z})$ are the Pauli matrices,
$\sigma_{0}$ the identity, and
$w(\mathbf{r})=2w\sum_{j=1}^{3}\cos\left(\mathbf{q}_{j}\cdot\mathbf{r}\right)$
is the hexagonal superlattice potential, with amplitude $2w$ and wave vectors
$\mathbf{q}_{1,2,3}$ illustrated as red vectors in Fig. 2(a). The wave vectors
satisfy $|\mathbf{q}_{1,2,3}|=4\pi/\sqrt{3}L$, where $L$ is the periodicity of
the potential. We set the Fermi velocity to $v_{F}=2.55\text{eV}\rm\AA$, the
experimentally measured value in $\text{Bi}_{2}\text{Te}_{3}$ [8, 7]. For the
moment we neglect higher-order corrections [5] to the Dirac cone dispersion
(1). Then Eq. (1) depends only on a single dimensionless parameter
$w/(v_{F}/L)$.

Figure 2: Superlattice Brillouin zone, electronic properties of surface states
and magnetic ordering basis. (a) Black and cyan lines show the superlattice
and magnetic Brillouin zones, respectively. The red and green arrows represent
the wave vectors $\mathbf{q}_{j}$ and the nesting vectors, $\bm{\delta}_{j}$,
respectively. (b) Density of states and electronic dispersion for
$w/(v_{F}/L)\simeq 1.54$ ($L=20$nm, $w=20$meV). Dashed lines show the density
of states and the band structure of the TI Dirac cone folded into the
superlattice Brillouin zone. (c) Momentum resolved spectral function
displaying the Fermi surface at the energy of the VHS with $w/(v_{F}/L)\simeq
1.54$. (d) Decomposition of the order parameter into $\bm{v}_{\mu j}$.

The superlattice potential has a profound effect on the Dirac cone. The
dispersion of the first band above charge neutrality [red line in Fig. 2(b)]
becomes nearly flat at the $\mathbf{K}$ and $\mathbf{K}^{\prime}$ points of
the moiré Brillouin zone. Specifically, near $\mathbf{K}$ and
$\mathbf{K}^{\prime}$, symmetry constrains the Taylor expansion of the
dispersion to the form $\epsilon_{\pm}(\mathbf{k})=\alpha
k^{2}\pm\eta(k^{3}_{x}-3k_{x}k^{2}_{y})$, where $\pm$ indicates valley,
expressions for the coefficients $\alpha$ and $\eta$ are given in Appendix A.
This yields three VHSs near $\mathbf{K}$ and $\mathbf{K}^{\prime}$ that map
into each other under $C_{3z}$ and provide the large diverging density of
states (DOS) shown in Fig. 2(c). In addition, there is a local maximum
(minimum) for $\alpha<0$ $(>0)$ at $\mathbf{k}=0$. At a critical value of
$w/(v_{F}/L)$ where $\alpha=0$, the three nearby VHSs merge to form higher-
order VHSs at $\mathbf{K}$ and $\mathbf{K}^{\prime}$ with a power-law
diverging DOS [52, 53, 54, 55], as found in Ref. [20]. Our results do not
require fine tuning to the critical value of $w/(v_{F}/L)$. Instead, we focus
on the vicinity of the VHSs that produce a sharp peak in the DOS [highlighted
by the horizontal line in Fig. 2(b)].

To describe the magnetic instability with momentum
$\mathbf{K}^{\prime}-\mathbf{K}$ we consider the electron-electron interaction
coupling the hot-spot regions around $\mathbf{K}$ and $\mathbf{K}^{\prime}$.
Note that screening by low-energy electrons at momentum transfer
$\mathbf{q}=\mathbf{K}^{\prime}-\mathbf{K}$ are weak because electrons at
$\mathbf{K}$ and $\mathbf{K}^{\prime}$ occupy orthogonal Bloch states (see
Appendix C).

For simplicity we focus on the Hubbard interaction:

$\hat{H}_{\text{int}}=U\int
d^{2}\mathbf{r}\,\hat{\Psi}^{\dagger}_{\uparrow}(\mathbf{r})\hat{\Psi}^{\dagger}_{\downarrow}(\mathbf{r})\hat{\Psi}_{\downarrow}(\mathbf{r})\hat{\Psi}_{\uparrow}(\mathbf{r});\quad
U>0,$ (2)

where $\hat{\Psi}_{\uparrow,\downarrow}(\mathbf{r})$ is the electron
annihilation operator,
$\hat{\Psi}_{\sigma}(\mathbf{r})=\sum_{\mathbf{k},\mathbf{G}}e^{i(\mathbf{k}-\mathbf{G})\mathbf{r}}\,\hat{c}_{\mathbf{k},\mathbf{G},\sigma}/\sqrt{A}$
with $A$ area of the sample, $\mathbf{k}$ and $\mathbf{G}$ label the Bloch
momentum and reciprocal lattice vector of the moiré Brillouin zone,
respectively, and $\hat{c}^{\dagger}_{\mathbf{k},\mathbf{G},\sigma}$ creates
an electron at $\mathbf{k}-\mathbf{G}$ with spin $\sigma$. Notice that Ref.
[20] considered an attractive interaction $(U<0)$ and found that the VHS
enhances the superconducting critical temperature.

###  II.1 Spin density wave operators

The interaction term $\hat{H}_{\text{int}}$ drives an instability towards
density wave ordering, modulated by the nesting vectors, $\bm{\delta}_{j}$,
which connect the $\mathbf{K}$ and $\mathbf{K}^{\prime}$ regions where the DOS
diverges at the VHS. The nesting vectors are indicated by green vectors in
Fig. 2(a),(c) and define a magnetic Brillouin zone [solid cyan line in Fig.
2(a)] three times smaller than the original one. Spin-momentum locking is
expected to drive the formation of an exotic spin-density wave (SDW) [41]. To
study the SDW, we decompose the order parameters for each modulation
$\bm{\delta}_{j}$ [5] into an in-plane direction parallel to the nesting
vector, $\bm{v}_{\parallel j}=\bm{e}_{j}\equiv\bm{\delta}_{j}/|\delta_{j}|$;
an in-plane direction perpendicular to the nesting vector, $\bm{v}_{\perp
j}=-i(\mathbf{z}\times\bm{e}_{j})$; and an out-of plane part,
$\bm{v}_{zj}=-i\mathbf{z}$. Correspondingly, for each wave vector
$\bm{\delta}_{j}$, the SDW operator can be decomposed as

$\hat{S}_{\mu
j}=\sum_{\mathbf{k}}\sum_{\mathbf{G}}\hat{c}^{\dagger}_{\mathbf{k}+\bm{\delta}_{j},\mathbf{G},\sigma}\,\bm{v}_{\mu
j}\cdot\bm{\sigma}_{\sigma\sigma^{\prime}}\,\hat{c}_{\mathbf{k},\mathbf{G},\sigma^{\prime}},$
(3)

which represents the three types of magnetic order illustrated in Fig. 2(d).
We observe that we have $9$ independent spin-density waves corresponding to
different combinations of $\bm{\delta}_{j}$ and $\bm{v}_{\mu j}$.

SDW | $C_{3z}$ | $m_{x}$ | $T$ | $m_{y}$
---|---|---|---|---
$\hat{S}_{\mu 1}$ | $\hat{S}_{\mu 2}$ | $\hat{S}^{\dagger}_{\mu 1}$ | $-\hat{S}^{\dagger}_{\mu 1}$ | if $z,\parallel:-\hat{S}_{\mu 1}$; else $\hat{S}_{\perp 1}$
$\hat{S}_{\mu 2}$ | $\hat{S}_{\mu 3}$ | $\hat{S}^{\dagger}_{\mu 3}$ | $-\hat{S}^{\dagger}_{\mu 2}$ | if $z,\parallel:-\hat{S}_{\mu 3}$; else $\hat{S}_{\perp 3}$
$\hat{S}_{\mu 3}$ | $\hat{S}_{\mu 1}$ | $\hat{S}^{\dagger}_{\mu 2}$ | $-\hat{S}^{\dagger}_{\mu 3}$ | if $z,\parallel:-\hat{S}_{\mu 2}$; else $\hat{S}_{\perp 2}$
Table 1: Symmetries of the SDW operator $\hat{S}_{\mu j}$ playing a key role
in determining the Ginzburg-Landau free energy. The transformation $m_{y}$
acts differently on $\mu=\parallel,z$ and $\perp$.

###  II.2 Symmetries

The interacting Hamiltonian $\hat{H}$ composed by $\hat{H}_{0}$ (1) and
$\hat{H}_{\rm int}$ (2) is invariant under time-reversal $T=\sigma_{y}K$,
where $K$ indicates complex conjugation, as well as a three-fold rotational
symmetry, $C_{3z}=e^{i\pi\sigma_{z}/3}$ and the mirror symmetries,
$m_{y}=i\sigma_{y}$ and $m_{x}=i\sigma_{x}$, which act on both spin and
spatial coordinates. Combinations of these imply invariance under two- and
six-fold rotation symmetries, $C_{2z}=i\sigma_{z}$,
$C_{6z}=e^{i\pi\sigma_{z}/6}$, respectively. In order to determine the leading
spin-density wave instability it is crucial to notice that the Bloch states
$\ket{u_{\mathbf{K}}}$ and $\ket{u_{\mathbf{K}^{\prime}}}$ of the Hamiltonian
(1) at the hot spots are singly degenerate and have opposite mirror $m_{y}$
eigenvalues. This can be readily understood observing that the time-reversal
symmetry $T$ sends $\mathbf{K}\to\mathbf{K}^{\prime}$ and
$Tm_{y}T^{-1}=m_{y}$. Thus, given $\ket{u_{\mathbf{K}}}$ with
$m_{y}\ket{u_{\mathbf{K}}}=i\ket{u_{\mathbf{K}}}$, it follows that
$m_{y}\ket{u_{\mathbf{K}^{\prime}}}=m_{y}T\ket{u_{\mathbf{K}}}=-iT\ket{u_{\mathbf{K}}}=-i\ket{u_{\mathbf{K}^{\prime}}}$.
The potential $w(\mathbf{r})$ breaks particle-hole symmetry. The action of the
symmetries of the model on the spin density wave operators $\hat{S}_{\mu j}$
is illustrated in Table 1.

##  III Ginzburg-Landau mean field theory

We study the magnetic order by deriving a Ginzburg-Landau theory. We decouple
the local interaction (2) introducing the Hubbard-Stratonovich fields
$n(\mathbf{r},\tau)$ and $\bm{m}(\mathbf{r},\tau)$ [56, 57, 58] for the charge
and the magnetization densities to obtain the Lagrangian:

$\begin{split}&\mathcal{L}=\int
d^{2}\mathbf{r}\Psi^{\dagger}\left(\partial_{\tau}+H_{0}(\mathbf{r})-\mu\right)\Psi\\\
&+\frac{U}{2}\int
d^{2}\mathbf{r}\left[n\,\rho-\bm{m}\cdot\bm{S}\right]+\frac{U}{4}\int
d^{2}\mathbf{r}\left[\bm{m}^{2}-n^{2}\right],\end{split}$ (4)

where for the sake of space we omit the dependencies of the fields on the
imaginary time $\tau$ and position $\mathbf{r}$, and we have introduced the
fermionic operators
$\hat{\rho}(\mathbf{r})=\sum_{\sigma}\hat{\Psi}^{\dagger}_{\sigma}(\mathbf{r})\hat{\Psi}_{\sigma}(\mathbf{r})$
and
$\hat{S}_{a}(\mathbf{r})=\bm{\hat{}}\Psi^{\dagger}(\mathbf{r})\sigma_{a}\hat{\Psi}(\mathbf{r})$.
Integrating out the electronic degrees of freedom of Eq. (LABEL:HS_action) we
obtain the effective action for $n(\tau,\mathbf{r})$ and
$\bm{m}(\mathbf{r},\tau)$. Then, we take the semiclassical limit of a static
order parameter. More specifically, we assume a homogenous electron density
$n$ and a spatially-modulated magnetization
$\bm{m}(\mathbf{r})=\sum^{3}_{j=1}\left(\bm{m}_{\bm{\delta}_{j}}e^{i\bm{\delta}_{j}\mathbf{r}}+\bm{m}_{-\bm{\delta}_{j}}e^{-i\bm{\delta}_{j}\mathbf{r}}\right)$
where $\bm{m}_{\bm{\delta}_{j}}$ is the Fourier amplitude associated to the
modulation $\bm{\delta}_{j}$. Since $\bm{m}(\mathbf{r})$ is real,
$\bm{m}_{\bm{\delta}_{j}}=\bm{m}^{*}_{-\bm{\delta}_{j}}\equiv\bm{m}_{j}$. We
expand $\bm{m}_{j}=\sum_{\mu}s_{\mu j}\bm{v}_{\mu j}$ where $\bm{v}_{\mu j}$
are the normal modes in Fig. 2(d) and the order parameter $s_{\mu
j}=\langle\hat{S}_{\mu j}\rangle/A$ is the average value of $\hat{S}_{\mu j}$
with $A$ area of the sample. As a result we find the free energy:

$F=\frac{A}{2}\left(\sum_{\mu j}|s_{\mu
j}|^{2}-\frac{n^{2}}{2}\right)+\sum_{m=1}^{\infty}\Tr\left[\frac{\left(-G_{0}\,X\right)^{m}}{m\,U}\right],$
(5)

where $G_{0}=(-\partial_{\tau}-H_{0}+\mu)^{-1}$ is the non-interacting single-
particle Green’s function and $X=U\bm{m}(\mathbf{r})\cdot\bm{\sigma}/2-Un/2$
describes the interaction between the order parameter and the electrons.

The self-consistency equations are obtained by minimizing the free energy (5)
with respect to the variational parameters $n$ and $s_{\mu j}$, i.e., $\delta
F/\delta n=0$ and $\delta F/\delta s^{*}_{\mu j}=0$. The magnetic order is
determined by finding the roots of these equations numerically, as detailed in
Appendix B. We find that in the range of filling between $n\sim 0.61$ and
$n\sim 0.8$ the minimum of the free energy (5) develops a SDW magnetic
ordering. The blue data in Fig. 3(a) shows the magnitude of the total
magnetization density, $|M|=\sqrt{\sum_{\mu j}|s_{\mu j}|^{2}}$, as a function
of the filling per unit cell of the superlattice.

In real space, the magnetic order $\bm{m}(\mathbf{r})$ obtained from the
saddle-point solution forms a Néel-type meron lattice [25] with two merons in
each unit cell, as shown in Fig. 3(b). The merons are well-defined because the
magnetic moment is forced to be in-plane along mirror-invariant lines. Each
meron can be characterized by the Pontryagin density
$\Phi(\mathbf{r})=\,\hat{\bm{m}}(\mathbf{r})\cdot\left[\partial_{x}\hat{\bm{m}}(\mathbf{r})\times\partial_{y}\hat{\bm{m}}(\mathbf{r})\right]/4\pi$
shown in Fig. 3(c) where $\hat{\bm{m}}$ is a unit magnetization vector. We
find $\int_{S}d^{2}\mathbf{r}\Phi(\mathbf{r})=\pm 1/2$, where $S$ is the
triangular domain illustrated by the orange dashed line in Fig. 3(c). In the
center of each unit cell is an intermediate region of destructive interference
with vanishing magnetization. The magnetic order spontaneously breaks the
translation symmetry of the potential, giving rise to a magnetic unit cell
[cyan line in Fig. 3(b)] three times larger ($\sqrt{3}\times\sqrt{3}$) than
the original one [black solid line in Fig. 3(b)].

Figure 3: Saddle-point solution, real space magnetization and electronic
structure in the meron lattice phase. (a) Blue and red data show the total
magnetization $|M|$ and the lowest eigenvalue of $\mathcal{M}^{\mu\nu}_{jl}$,
respectively. (b) Contour plot of the magnetization in real space. Color
indicates the out-of-plane component. (c) Winding number
$\bm{m}\cdot(\partial_{x}\bm{m}\times\partial_{y}\bm{m})/4\pi$ in real space.
Black and cyan lines show the superlattice and magnetic unit cells,
respectively. (d) Density of states and electronic dispersion of the meron
lattice state. The dashed line shows the density of states of the normal
phase. The horizontal line shows the Fermi energy $E_{F}$; $\Delta E$ denotes
the induced gap at $\bm{\Gamma}_{m}$. Inset shows the Fermi surface of the
magnetic state. The bands are obtained from the Hartree-Fock solution at
$w/(v_{F}/L)\simeq 1.54$, $2w=40\text{meV}$, $U=30\,\text{meV}$ and
$T=1.5\text{K}$.

###  III.1 Hartree-Fock Hamiltonian

In this section we discuss the properties of the electronic band structure in
the magnetic phase. We treat the effect of the magnetization on the electronic
spectrum at the Hartree-Fock level by replacing $H_{0}(\mathbf{r})$ in Eq. (1)
with:

$H_{\text{HF}}(\mathbf{r})=H_{0}(\mathbf{r})-\frac{U}{2}\bm{m}(\mathbf{r})\cdot\bm{\sigma}.$
(6)

Despite breaking time-reversal symmetry, the resulting electronic spectrum
remains gapless because the magnetic order preserves $C_{2z}T$ symmetry. The
dispersion is plotted in the right panel of Fig. 3(d) in the magnetic
Brillouin zone: the original red band in Fig. 2(b) is decomposed into three
different bands given by the red, blue and green lines in Fig. 3(d). The Dirac
cones at $\mathbf{K}_{m}$ between the red and blue bands and at
$\mathbf{K}^{\prime}_{m}$ between the blue and the green ones are protected by
$C_{2z}T$. Thus, they give rise to a gapless electronic spectrum whose Fermi
surface is shown in the inset of Fig. 3(d).

Although the order parameter $\bm{m}(\mathbf{r})$ does not open a full gap, it
significantly reduces the DOS at the Fermi level (horizontal black dashed line
in Fig. 3(d)) by opening a gap between the blue and green bands at
$\Gamma_{m}$ of order $\Delta E\propto U|M|$. The gap splits the peak in the
DOS resulting from the VHS into two peaks above and below, as shown in the
left panel of Fig. 3(d). The significant decrease in the kinetic energy from
splitting the large peak in the DOS makes the magnetic state energetically
favorable with respect to the normal one.

In the following we discuss the symmetries of the magnetic phase and we
determine the region of stability of the spin density wave order.

##  IV Phase diagram

The instability of the Dirac cone surface state is determined by expanding the
free energy for small values of the order parameter $s_{\mu j}$. As time-
reversal symmetry acts on the order parameter by $T:s_{\mu j}\to-s^{*}_{\mu
j}$ [see Table 1], only even powers of $s_{\mu j}$ are allowed in the free
energy (5), so that to second order in $s_{\mu j}$ at fixed density $n$,
$F_{(2)}=\sum_{\mu\nu}\sum^{3}_{jl=1}s^{*}_{\mu
j}\mathcal{M}^{\mu\nu}_{jl}s_{\nu l}/2$, where

$\mathcal{M}^{\mu\nu}_{jl}=\delta_{\mu\nu}\delta_{jl}-\frac{U}{2}\chi^{\mu\nu}_{jl},$
(7)

and
$\chi^{\mu\nu}_{jl}\equiv\int^{\beta}_{0}d\tau\langle\hat{S}^{\dagger}_{\mu
j}(\tau)\hat{S}_{\nu l}(0)\rangle/A$ is the thermodynamic spin susceptibility
with $\tau$ imaginary time and $\beta$ inverse temperature. The matrix
$\mathcal{M}^{\mu\nu}_{jl}$ contains both diagonal terms ($j=l$), which come
from momentum-conserving scattering processes, and off-diagonal terms ($j\neq
l$) from Umklapp processes with momentum
$\bm{\delta}_{j}-\bm{\delta}_{l}\in\mathbf{G}$ with $\mathbf{G}$ a reciprocal
lattice vector. The matrix $\mathcal{M}^{\mu\nu}_{jl}$ is constrained by the
symmetries of the Hamiltonian as discussed in Appendix C.

An instability exists when one of the eigenvalues of the matrix
$\mathcal{M}^{\mu\nu}_{jl}$ becomes negative, which comprises a generalization
of the Stoner criterion [59]. The corresponding eigenvector indicates the
magnetic configuration of the instability and is classified by how it
transforms under symmetry. Since the saddle-point solution is
$C_{3z}$-symmetric, we consider a $C_{3z}$-invariant eigenvector of
$\mathcal{M}^{\mu\nu}_{jl}$, which implies the magnetization takes the form
$s_{\mu j}=MU_{\mu}/\sqrt{3}$. It follows that

$F_{(2)}=\frac{|M|^{2}}{2}\left(\begin{matrix}U_{\parallel}\\\ U_{\perp}\\\
U_{z}\end{matrix}\right)^{\dagger}\left(\begin{matrix}L_{\parallel\parallel}&0&L_{\parallel
z}\\\ 0&L_{\perp\perp}&0\\\ L_{\parallel
z}&0&L_{zz}\end{matrix}\right)\left(\begin{matrix}U_{\parallel}\\\
U_{\perp}\\\ U_{z}\end{matrix}\right),$ (8)

where the elements $L_{\mu\nu}=\sum_{jl}\mathcal{M}^{\mu\nu}_{jl}/3$ are shown
in Fig. 4(a) as a function of the density at $T=1.5$K. The perpendicular
($U_{\perp}$) component of the magnetic order is decoupled from the other
components ($U_{\parallel},U_{z}$) because it is even under the mirror
$m_{y}$, while the parallel and $z$-components are odd. This property can be
readily understood by looking at the last column of Table 1, where we show the
action of the symmetry $m_{y}$ on the different components of the spin-density
wave order. From Eq. (8), the surface Dirac cone is unstable when either
$\lambda_{1}=L_{\perp\perp}<0$, which corresponds to in-plane order purely in
the “spiral-xy” channel and even under $m_{y}$, or when one of the eigenvalues
$\lambda_{\pm,2}=(L_{\parallel\parallel}+L_{zz})/2\pm\sqrt{(L_{\parallel\parallel}-L_{zz})^{2}/4+L_{\parallel
z}^{2}}<0$, which corresponds to a generically non-coplanar SDW with “120-xy”
and $z$ components, and breaking $m_{y}$ [see Fig. 2(d)]. In order to
characterize the magnetic instability we compute the spin-susceptibility
introduced in Eq. (7):

$\begin{split}\chi^{\mu\nu}_{jl}=&-\frac{T}{A}\sum_{\mathbf{k},\epsilon_{n}}\Tr[G_{0}(\mathbf{k},i\epsilon_{n})\,O^{\dagger}_{\mu
j}\,V^{\bm{\delta}_{j}-\bm{\delta}_{l}}\\\
\hfil&G_{0}(\mathbf{k}+\bm{\delta}_{l},i\epsilon_{n})\,O_{\nu
l}\Big{]},\end{split}$ (9)

where $G_{0}(\mathbf{k},i\epsilon_{n})$ is the single particle Green’s
function

$G_{0}(\mathbf{k},i\epsilon)=\sum_{n}\frac{\ket{u_{n\mathbf{k}}}\bra{u_{n\mathbf{k}}}}{i\epsilon-\xi_{n\mathbf{k}}},$
(10)

where $\ket{u_{n\mathbf{k}}}$ are the Bloch states of $H_{0}$, for simplicity
we have introduced the operator $O_{\mu j}=\bm{v}_{\mu j}\cdot\bm{\sigma}$ and
$V^{\mathbf{G}}$ is the sewing matrix
$V^{\mathbf{G}}_{\mathbf{Q},\mathbf{Q}^{\prime}}=\delta_{\mathbf{Q},\mathbf{Q}^{\prime}+\mathbf{G}}$.
We find that for a range of filling near $2/3$ electron per unit cell, the
dominant instability occurs in the channel corresponding to $\lambda_{-,2}$,
as shown by the region where the red solid line in Fig. 3(a) goes negative,
and in agreement with our saddle-point solution. This solution, which breaks
$m_{y}$, is energetically favored over the in-plane spin solution that is even
in $m_{y}$ because it allows a gap to open at $\bm{\Gamma}_{m}$ in the
magnetic BZ, indicated by $\Delta E$ in Fig. 3(d). Indeed, the latter solution
being odd under $m_{y}$ allows mixing between the Bloch states
$\ket{u_{\mathbf{K}}}$ and $\ket{u_{\mathbf{K}^{\prime}}}$ at the hot spot
regions which are characterized by opposite mirror $m_{y}$ eigenvalues as
detailed in Appendix C. The resulting magnetic state is a pattern of half-
integer topological vortices, with the same winding and opposite polarities
corresponding to a meron lattice. The other magnetic state, with eigenvalue
$\lambda_{+,2}$, also breaks $m_{y}$, but is energetically unfavourable
because it exhibits a clockwise magnetization winding inconsistent with the
intrinsic SOC of the surface state in Eq. (1). We have also studied the
$C_{3z}$-breaking magnetic orders. Since these orders do not hybridize the
states $\ket{u_{\mathbf{K}}}$ and $\ket{u_{\mathbf{K}^{\prime}}}$, they do not
lead to a magnetic instability. The symmetry of the magnetic orders are
further discussed in Appendix C.

Figure 4: Magnetic instability, critical temperature and interaction. (a)
Solid lines show the components of the second order free energy $L_{\mu\nu}$.
Red and cyan dashed lines illustrate the eigenvalues $\lambda_{-,2}$ and
$\lambda_{+,2}$, respectively, while $\lambda_{1}=L_{\perp\perp}$ is the solid
green line. The calculations were performed at $T=1.5\text{K}$,
$U=30\text{meV}$ and $w/(v_{F}/L)\simeq 1.54$. (b) Critical temperature for
the magnetic transition as a function of $w/(v_{F}/L)$ for $U=30\text{meV}$.
For a given $w/(v_{F}/L)$, the filling $n$ is set so that the Fermi level is
at the energy of the VHS. (c) Evolution of the energy of the VHS singularity
as a function of $w/(v_{F}/L)$. Inset shows the electron filling per
superlattice unit cell at the VHS versus $w/(v_{F}/L)$. (d) Critical
interaction $U_{c}$ as a function of the filling $n$ at several temperatures,
fixing $2w=40\text{meV}$ and $w/(v_{F}/L)\simeq 1.54$.

The critical temperature and interaction strength occur at precisely the point
where one of the eigenvalues of $\mathcal{M}^{\mu\nu}_{jl}$ changes sign. Fig.
4(b) shows the critical temperature for a range of fillings around the VHS at
$U=30$meV. The magnetic dome peaks at $T_{c}\sim 7$K near an optimal value of
$w/(v_{F}/L)\simeq 1.54$: surprisingly, the peak does not coincide with the
higher-order VHS [20] at $w/(v_{F}/L)\simeq 1.36$. Fig. 4(c) shows the
evolution of the energy and the number of electrons per superlattice unit cell
where the VHS occurs as a function of $w/(v_{F}/L)$. The maximum critical
temperature for the magnetic transition takes place around $n\sim 2/3$
filling, denoted as the horizontal green line in the inset of Fig. 4(c).

Fig. 4(d) shows the critical interaction $U_{c}$ necessary to induce the
magnetic instability as a function of the filling $n$ for different
temperatures $T$. Lowering $T$ reduces the magnetic transition to arbitrarily
weak repulsive interactions. In contrast, in the absence of the superlattice
potential, spontaneous magnetization of the surface state [5, 13, 14, 12, 15,
16, 17, 18] in the $\rm Bi_{2}\rm Se_{3}$ family requires a critical value of
the Coulomb interaction, on the order of $2$eV [13, 14, 18]. We expect our
results are robust to weak disorder, which will reduce the transition
temperature [60, 61, 62] without eliminating the magnetic state.

###  IV.1 Hexagonal warping

Beyond the linear momentum dependence of the Dirac cone (1) the dispersion
develops a hexagonal warping term
$H_{w}(\mathbf{k})=\lambda(k^{3}_{x}-3k_{x}k^{2}_{y})\sigma_{z}$ [5]. Despite
being small in $\lambda/v_{F}L^{2}$ this correction gives rise to a series of
interesting effects that are experimentally relevant. The hexagonal warping
breaks $m_{y}$ and $C_{2z}T$ symmetries which gives a finite $U_{\perp}$
component in the magnetization pattern. While the average magnetic moment in
the unit cell continues to vanish, $\int d^{2}\mathbf{r}\leavevmode\nobreak\
m_{i}=0$, the explicitly broken $m_{y}$ symmetry gives rise to a finite out-
of-plane toroidal moment $\mathcal{T}_{z}=\int
d^{2}\mathbf{r}\,(xm_{y}-ym_{x})/2$ [63, 64, 65] which typically manifests in
the magnetoelectric susceptibility $\alpha_{xy}=-\alpha_{yx}$ [66, 67].

##  V Conclusions

We have shown that a superlattice potential on the surface of a TI provides a
new route to spontaneously breaking time reversal symmetry on the TI surface.
The magnetic order realizes a meron lattice exhibiting pairs of merons with
opposite topological charge in the unit cell. The periodicity of the meron
lattice is determined by the period of the potential; consequently, the meron
lattice periodicity is completely tunable. Although the magnetic order breaks
time-reversal, it preserves $C_{2z}T$ and thus does not open a gap on the TI
surface.

The meron lattice phase can be measured by imaging the magnetization in real
space through Lorentz transmission electron microscopy [68, 26] and nitrogen
vacancy magnetometry [69] or in reciprocal space via X-ray diffraction [22,
70, 71, 72, 73, 72, 74, 75]. The magnetization can also be measured through
the magneto-optical Kerr effect [76, 77] and reflective magnetic circular
dichroism [78]. In addition, the magnetoelectric susceptibility $\alpha_{xy}$
is observed either by measuring an in-plane magnetization in response to an
electric field, or a current resulting from an in-plane magnetic field. The
reduced density of states resulting from the magnetic order gives a drastic
variation of the electronic compressibility across the transition. Finally, we
also expect that the response to a magnetic field [79, 80, 81, 82, 83, 84, 85]
gives rise to distinctive signatures of the meron state.

This finding has far-reaching consequences, as this unconventional magnetic
state will have implications both on the experimental and on the theoretical
level. First, the topological spin texture realizes an electromagnetic field
on the scale of the superlattice that can be employed by proximity effects to
modify band structure and topological properties of electronic systems [35,
86, 87, 37, 33, 88, 89]. Second, the broken mirror symmetry that results from
hexagonal warping of the Fermi surface implies that the meron lattice might
host a magnetoelectric response with potential application to spintronics [90,
91, 92, 93, 88, 89]. This effect can be enhanced by strain [94, 95, 96] and
lattice relaxation [97]. Finally, the interplay between the spin density wave
ordering and possible superconducting instabilities [20] is an open problem
which is left to future studies.

###### Acknowledgements.

_Acknowledgements.—_ We have benefited from discussions with Yang-Zhi Chou, M.
Michael Denner, Shiang Fang, Jiawei Zang, Andrew J. Millis, Zhentao Wang,
Tiancheng Song, and Justin Wilson. We are grateful to Lucy Reading-Ikkanda for
creating the figure of the experimental setup and the sketch of the magnetic
state. We also acknowledge the support of the Flatiron Institute, a division
of the Simons Foundation. This work was partially supported by the Air Force
Office of Scientific Research under Grant No. FA9550-20-1-0260 (J.C.) and
Grant No. FA9550-20-1-0136 (J.H.P.) and the Alfred P. Sloan Foundation through
a Sloan Research Fellowship (J.H.P.). J.H.P. and J.C. acknowledge hospitality
of the Aspen Center for Physics, where some of this work was developed and
which is supported by National Science Foundation grant PHY1607611.

##  Appendix A The $\mathbf{k}\cdot\mathbf{p}$ expansion around $\mathbf{K}$
and $\mathbf{K}^{\prime}$

The origin of the high order van Hove singularities (HOVHS) [52, 53, 54, 55]
can be derived from the lattice symmetries of the Hamiltonian
$H_{0}(\mathbf{r})$. We define the small deviation from the $\mathbf{K}$
point, $\mathbf{q}=\mathbf{k}-\mathbf{K}$, and the linear combinations
$q_{\pm}=q_{x}\pm iq_{y}$, which transforms under $C_{3z}$ as $q_{\pm}\to
e^{\pm i2\pi/3}q_{\pm}$. At second order in $\mathbf{q}$ the only allowed term
in the energy dispersion is $q_{+}q_{-}=q^{2}$. The next contributions
invariant under $C_{3z}$ are $q^{3}_{+}$ and $q^{3}_{-}$. As a consequence of
the mirror $m_{y}$ symmetry the cubic term is
$(q^{3}_{+}+q^{3}_{-})/2=q^{3}_{x}-3q^{2}_{y}q_{x}$.

Figure 5: Parameters $\alpha$ and $\eta$ of the $\mathbf{k}\cdot\mathbf{p}$
expansion of the dispersion relation around $\mathbf{K}$. The vanishing of
$\alpha$ at $w/(v_{F}/L)=1.36$ highlighted by the green line implies a HOVHS.

Thus, to cubic order in $\mathbf{q}$ the dispersion around $\mathbf{K}$ reads
[20] $\epsilon_{\mathbf{K}}(\mathbf{q})=\epsilon_{0}+\alpha
q^{2}+\eta(q^{3}_{x}-3q_{x}q^{2}_{y})+\cdots$. Applying time-reversal symmetry
we find that the expansion at $\mathbf{K}^{\prime}$ reads
$\epsilon_{\mathbf{K}^{\prime}}(\mathbf{q})\simeq\epsilon_{0}+\alpha
q^{2}-\eta(q^{3}_{x}-3q_{x}q^{2}_{y})$. The values of the coefficients
$\alpha$ and $\eta$ are obtained from the $\mathbf{k}\cdot\mathbf{p}$
perturbation theory at $\mathbf{K}$. Given
$H_{0}(\mathbf{K}+\mathbf{q})-H_{0}(\mathbf{K})$ expanding to the third order
in the small deviation $\mathbf{q}=\mathbf{k}-\mathbf{K}$ yields:

$\alpha=2\sum_{m\neq
n}\frac{\real\left[\matrixelement{u_{n\mathbf{K}}}{\sigma_{+}}{u_{m\bm{K}}}\matrixelement{u_{m\mathbf{K}}}{\sigma_{-}}{u_{n\mathbf{K}}}\right]}{\epsilon_{n\mathbf{K}}-\epsilon_{m\mathbf{K}}},$
(11)

and

$\begin{split}\eta=&\sum_{m,j\neq
n}2\imaginary[\frac{\matrixelement{u_{n\mathbf{K}}}{\sigma_{-}}{u_{m\mathbf{K}}}\matrixelement{u_{m\mathbf{K}}}{\sigma_{-}}{u_{j\mathbf{K}}}}{(\epsilon_{n\mathbf{K}}-\epsilon_{m\mathbf{K}})(\epsilon_{n\mathbf{K}}-\epsilon_{j\mathbf{K}})}\\\
&\matrixelement{u_{j\mathbf{K}}}{\sigma_{-}}{u_{n\mathbf{K}}}\Bigg{]},\end{split}$
(12)

where $\ket{u_{n\mathbf{K}}}$ and $\epsilon_{n\mathbf{K}}$ are the Bloch state
and the eigenvalue, respectively, of $H_{0}(\mathbf{r})$ at $\mathbf{K}$, $n$
is the first positive energy band and $\sigma_{\pm}=\sigma_{x}\pm
i\sigma_{y}$. The evolution of the couplings $\alpha$ and $\eta$ is shown in
Fig. 5. The higher-order van Hove singularity occurs when the quadratic term
$\alpha$ in Eq. (11) vanishes at $w/(v_{F}/L)\simeq 1.36$, highlighted by the
vertical green line in Fig. 5. The result is a power-law divergence in the
density of states $\rho(\epsilon)\sim|\epsilon|^{-1/3}$. By looking at the
solutions of $\nabla_{\mathbf{q}}\epsilon_{\pm}(\mathbf{q})=0$ ($\pm$ for
$\mathbf{K}$ and $\mathbf{K}^{\prime}$ respectively), we find that
$\mathbf{q}=0$ is a maximum for $\alpha<0$, a minimum for $\alpha>0$ and,
finally, a higher-order critical point for $\alpha=0$. Away from the origin
there are three further solutions at $\bm{\kappa}_{\pm,j}=\mp 2\alpha
C^{j-1}_{3z}(1,0)/3\eta$, which are saddle points. Approaching the higher-
order VHS at $\alpha\to 0$, these three saddle points merge at $\mathbf{q}=0$,
i.e., at $\mathbf{K}$ or $\mathbf{K}^{\prime}$.

##  Appendix B Solution of the saddle-point equations

In this section we detail the mean-field solution of the interacting electrons
on the surface of the TI subject to the superlattice potential
$w(\mathbf{r})$. Minimizing $F$ (5) with respect to the filling $n$ and the
magnetic configurations $\bm{m}(\mathbf{r})$ we find the saddle-point
equations:

$\frac{\delta F}{\delta n}=0\Rightarrow
n=\frac{1}{A}\Tr\left[\left(G_{0}^{-1}+X\right)^{-1}\right],$ (13)

and

$\frac{\delta F}{\delta s^{*}_{\mu j}}=0\Rightarrow s_{\mu
j}=\frac{2}{UA}\Tr\left[\frac{\delta X}{\delta s^{*}_{\mu
j}}\left(G^{-1}_{0}+X\right)^{-1}\right].$ (14)

In the main text we show the result of the numerical solution of Eqs. (13) and
(14), which was obtained as follows. Eqs. (13) and (14) take a simple form in
the basis of the eigenstates $\ket{\phi_{n\mathbf{k}}}$ of the Hartree-Fock
Hamiltonian in Eq. (6) with eigenvalues $\bar{\epsilon}_{n\mathbf{k}}$:

$n=\frac{1}{3N}\sum_{n\mathbf{k}}f(\bar{\xi}_{n\mathbf{k}})-n_{0},$ (15)

where $\bar{\xi}_{n\mathbf{k}}=\bar{\epsilon}_{n\mathbf{k}}-\mu^{*}$ and
$\mu^{*}=\mu-Un/2$; and

$s_{\mu
j}=\frac{1}{3N}\sum_{n\mathbf{k}}f(\bar{\xi}_{n\mathbf{k}})\matrixelement{\phi_{n\mathbf{k}}}{O^{\dagger}_{\mu
j}V^{\bm{\delta}_{j}}}{\phi_{n\mathbf{k}}},$ (16)

where the filling $n$ is measured with respect to the charge neutrality point
$n_{0}$, $N$ is the number of points sampling the Brillouin zone, the factor
of three takes into account the size difference between the original and the
magnetic Brillouin zone. Moreover, the function
$f(\epsilon)=1/(e^{\beta\epsilon}+1)$ is the Fermi-Dirac distribution
function, and $V^{\mathbf{G}}$ is the sewing matrix satisfying
$V^{\mathbf{G}}_{\mathbf{Q},\mathbf{Q}^{\prime}}=\delta_{\mathbf{Q},\mathbf{Q}^{\prime}+\mathbf{G}}$,
$V^{\mathbf{G}}\ket{\phi_{n\mathbf{k}}}=\ket{\phi_{n\mathbf{k}+\mathbf{G}}}$.
The self-consistent equations is performed by a find-root algorithm.

##  Appendix C Symmetries of the free energy and classification of magnetic
orderings

In this Appendix we discuss the symmetries of the second order matrix
$\mathcal{M}^{\mu\nu}_{jl}$ and classify the different magnetic orderings
according to their symmetry properties. The second order tensor is constrained
by the symmetries in Table 1. The three-fold rotation about the $z$-axis,
$C_{3z}:s_{\mu j}\to s_{\mu j+1}$, gives
$\mathcal{M}^{\mu\nu}_{jl}=\mathcal{M}^{\mu\nu}_{j+1l+1}$. As a result
$\mathcal{M}^{\mu\nu}_{jl}$ is expressed in terms of two three-dimensional
matrices: $\Lambda_{\mu\nu}=\mathcal{M}^{\mu\nu}_{j,j}$,
$\Omega_{\mu\nu}=\mathcal{M}^{\mu\nu}_{j,j+1}$ while
$\mathcal{M}^{\mu\nu}_{j+1,j}=\Omega^{*}_{\nu\mu}$ since the free energy
$F_{(2)}$ is real. In addition, the mirror symmetry $m_{x}$ [given in Table
(1)] implies that $\Lambda_{\mu\nu}$ is real and symmetric,
$\Omega_{\mu\nu}=\Omega_{\nu\mu}$ but with complex elements
$\Omega^{*}_{\mu\nu}\neq\Omega_{\nu\mu}$. Finally, the remaining $m_{y}$
symmetry imposes that the matrix elements $\Lambda_{\parallel\perp}=0$ and
$\Lambda_{\perp z}=0$, $\Omega_{\parallel\perp}$ and $\Omega_{\perp z}$ are
purely imaginary, while the other components of $\Omega_{\mu\nu}$ are reals.
The latter constraint implies that the $C_{3z}$ symmetric spiral-xy order
decouples from the 120-xy and out-of-plane orderings. As a result by
projecting the second order free energy $F_{(2)}$ on the $C_{3z}$ symmetric
configuration $s_{\mu j}=\xi\,U_{\mu}/\sqrt{3}$ we find the expression in the
main text Eq. (8) where $L_{\mu\nu}=\Lambda_{\mu\nu}+2\real\Omega_{\mu\nu}$.

| $E$ |  2$C_{3}$ |  3$m$
---|---|---|---
$\quad\Gamma_{1}\quad$ |  1 |  1 |  1
$\quad\Gamma_{2}\quad$ |  1 |  1 |  -1
$\quad\Gamma_{3}\quad$ |  2 |  -1 |  0
Table 2: Character table of $C_{3v}$. $E$, $C_{3}$ and $m$ represent the
conjugation classes of the identity, $C_{3z}$ and $m_{y}$, respectively.

We now characterize the spectrum of $\mathcal{M}^{\mu\nu}_{jl}$ and we
classify the magnetic states with respect to their symmetry behaviour under
the point-group symmetries of the model. The magnetic ordering is described by
the real order parameter
$\bm{m}(\mathbf{r})=\sum_{j}\bm{m}_{j}\,e^{i\bm{\delta}_{j}\cdot\mathbf{r}}+c.c.$
which is completely determined by the Fourier amplitudes
$\bm{m}_{j}=\sum_{\mu}s_{\mu j}\bm{v}_{\mu j}$. The nine magnetic
configurations originating from $s_{\mu j}$ are decomposed in one one-
dimensional representation $\Gamma_{1}$, two one-dimensional representations
$\Gamma_{2}$ and three two-dimensional representations $\Gamma_{3}$ given in
Table 2. The one-dimensional representations $\Gamma_{1}$ and $\Gamma_{2}$ are
symmetric under $C_{3z}$ and corresponds to equal superpositions of the
modulations $\bm{\delta}_{1,2,3}$. The parity of $\Gamma_{1}$ under $m_{y}$ is
even consistently with the spiral-xy order. On the other hand, $\Gamma_{2}$ is
odd under $m_{y}$ like the 120-xy and out-of-plane orderings. Due to the
intrinsic spin-orbit coupling in-plane 120-xy and out-of-plane orderings mix
to form the meron lattice state. Finally, we have the two-dimensional
representation $\Gamma_{3}$. Here, we find two eigenstates of $C_{3z}$
characterized by eigenvalues $\omega=e^{2\pi i/3}$ and $\omega^{*}$ and
$m_{y}$ is off-diagonal in this basis. Two dimensional degeneracy implies from
this symmetry.

In the normal state the two valleys $\mathbf{K}$ and $\mathbf{K}^{\prime}$ are
equally occupied giving rise to a time-reversal symmetric state with vanishing
magnetization. It is important to notice that the two valleys are
characterized by non-degenerate Bloch states, $\ket{u_{\mathbf{K}}}$ and
$\ket{u_{\mathbf{K}^{\prime}}}$, [first band above charge neutrality] with
opposite mirror eigenvalues. The meron lattice state belonging to $\Gamma_{2}$
breaks the mirror symmetry $m_{y}$ of $\mathbf{K}$ and $\mathbf{K}^{\prime}$
and couples $\ket{u_{\mathbf{K}}}$ and $\ket{u_{\mathbf{K}^{\prime}}}$. The
resulting Hartree-Fock orbitals at $\bm{\Gamma}_{m}$ of the magnetic moiré
Brillouin zone are eigenstates of the mirror symmetry $m_{x}$ and correspond
to a coherent superposition of the two valleys $\mathbf{K}$ and
$\mathbf{K}^{\prime}$.

## References

  * Fu  _et al._ [2007] L. Fu, C. L. Kane, and E. J. Mele, Topological insulators in three dimensions, Phys. Rev. Lett. 98, 106803 (2007).
  * Moore and Balents [2007] J. E. Moore and L. Balents, Topological invariants of time-reversal-invariant band structures, Phys. Rev. B 75, 121306 (2007).
  * Roy [2009] R. Roy, Topological phases and the quantum spin hall effect in three dimensions, Phys. Rev. B 79, 195322 (2009).
  * König  _et al._ [2007] M. König, S. Wiedmann, C. Brüne, A. Roth, H. Buhmann, L. W. Molenkamp, X.-L. Qi, and S.-C. Zhang, Quantum spin hall insulator state in hgte quantum wells, Science 318, 766 (2007), https://www.science.org/doi/pdf/10.1126/science.1148047 .
  * Fu [2009] L. Fu, Hexagonal warping effects in the surface states of the topological insulator Bi2Te3, Phys. Rev. Lett. 103, 266801 (2009).
  * Hsieh  _et al._ [2009a] D. Hsieh, D. Qian, L. Wray, Y. Xia, Y. S. Hor, R. J. Cava, and M. Z. Hasan, A topological dirac insulator in a quantum spin hall phase (experimental realization of a 3d topological insulator) (2009a), arXiv:0910.2420 [cond-mat.mes-hall] .
  * Hsieh  _et al._ [2009b] D. Hsieh, Y. Y. Xia, D. Qian, L. A. Wray, J. H. Dil, F. Meier, J. Osterwalder, L. Patthey, J. G. Checkelsky, N. P. Ong, A. V. Fedorov, H. I. Lin, A. Bansil, D. Grauer, Y. S. Hor, R. J. Cava, and M. Z. Hasan, A tunable topological insulator in the spin helical dirac transport regime, Nature 460, 1101 (2009b).
  * Chen  _et al._ [2009] Y. L. Chen, J. G. Analytis, J.-H. Chu, Z. K. Liu, S.-K. Mo, X. L. Qi, H. J. Zhang, D. H. Lu, X. Dai, Z. Fang, S. C. Zhang, I. R. Fisher, Z. Hussain, and Z.-X. Shen, Experimental realization of a three-dimensional topological insulator, bi2te3, Science 325, 178 (2009), https://www.science.org/doi/pdf/10.1126/science.1173034 .
  * Xia  _et al._ [2009] Y. Xia, D. Qian, D. Hsieh, L. A. Wray, A. Pal, H. Lin, A. Bansil, D. Grauer, Y. S. Hor, R. J. Cava, and M. Z. Hasan, Observation of a large-gap topological-insulator class with a single dirac cone on the surface, Nature Physics 5, 398 (2009).
  * Brüne  _et al._ [2011] C. Brüne, C. X. Liu, E. G. Novik, E. M. Hankiewicz, H. Buhmann, Y. L. Chen, X. L. Qi, Z. X. Shen, S. C. Zhang, and L. W. Molenkamp, Quantum hall effect from the topological surface states of strained bulk hgte, Phys. Rev. Lett. 106, 126803 (2011).
  * Souma  _et al._ [2011] S. Souma, K. Kosaka, T. Sato, M. Komatsu, A. Takayama, T. Takahashi, M. Kriener, K. Segawa, and Y. Ando, Direct measurement of the out-of-plane spin texture in the dirac-cone surface state of a topological insulator, Phys. Rev. Lett. 106, 216803 (2011).
  * Jiang and Wu [2011] J.-H. Jiang and S. Wu, Spin susceptibility and helical magnetic order at the edges/surfaces of topological insulators due to fermi surface nesting, Phys. Rev. B 83, 205124 (2011).
  * Baum and Stern [2012a] Y. Baum and A. Stern, Magnetic instability on the surface of topological insulators, Phys. Rev. B 85, 121105 (2012a).
  * Baum and Stern [2012b] Y. Baum and A. Stern, Density-waves instability and a skyrmion lattice on the surface of strong topological insulators, Physical Review B 86, 10.1103/physrevb.86.195116 (2012b).
  * Marchand and Franz [2012] D. J. J. Marchand and M. Franz, Lattice model for the surface states of a topological insulator with applications to magnetic and exciton instabilities, Phys. Rev. B 86, 155146 (2012).
  * Schmidt [2012] M. J. Schmidt, Strong correlations at topological insulator surfaces and the breakdown of the bulk-boundary correspondence, Phys. Rev. B 86, 161110 (2012).
  * Sitte  _et al._ [2013] M. Sitte, A. Rosch, and L. Fritz, Interaction effects on almost flat surface bands in topological insulators, Phys. Rev. B 88, 205107 (2013).
  * Mendler  _et al._ [2015] D. Mendler, P. Kotetes, and G. Schön, Magnetic order on a topological insulator surface with warping and proximity-induced superconductivity, Phys. Rev. B 91, 155405 (2015).
  * Cano  _et al._ [2021] J. Cano, S. Fang, J. H. Pixley, and J. H. Wilson, Moiré superlattice on the surface of a topological insulator, Phys. Rev. B 103, 155157 (2021).
  * Wang  _et al._ [2021] T. Wang, N. F. Q. Yuan, and L. Fu, Moiré surface states and enhanced superconductivity in topological insulators, Phys. Rev. X 11, 021024 (2021).
  * Gao  _et al._ [2019] N. Gao, S. G. Je, M. Y. Im, J. W. Choi, M. Yang, Q. Li, T. Y. Wang, S. Lee, H. S. Han, K. S. Lee, W. Chao, C. Hwang, J. Li, and Z. Q. Qiu, Creation and annihilation of topological meron pairs in in-plane magnetized films, Nature Communications 10, 10.1038/s41467-019-13642-z (2019).
  * Muhlbauer  _et al._ [2009] S. Muhlbauer, B. Binz, F. Jonietz, C. Pfleiderer, A. Rosch, A. Neubauer, R. Georgii, and P. Boni, Skyrmion lattice in a chiral magnet, Science 323, 915 (2009).
  * Pappas  _et al._ [2009] C. Pappas, E. Lelièvre-Berna, P. Falus, P. M. Bentley, E. Moskvin, S. Grigoriev, P. Fouquet, and B. Farago, Chiral paramagnetic skyrmion-like phase in MnSi, Phys. Rev. Lett. 102, 197202 (2009).
  * Neubauer  _et al._ [2009] A. Neubauer, C. Pfleiderer, B. Binz, A. Rosch, R. Ritz, P. G. Niklowitz, and P. Böni, Topological hall effect in the a phase of MnSi, Phys. Rev. Lett. 102, 186602 (2009).
  * Nagaosa and Tokura [2013] N. Nagaosa and Y. Tokura, Topological properties and dynamics of magnetic skyrmions, Nature nanotechnology 8, 899 (2013).
  * Heinze  _et al._ [2011] S. Heinze, K. von Bergmann, M. Menzel, J. Brede, A. Kubetzka, R. Wiesendanger, G. Bihlmayer, and S. Blügel, Spontaneous atomic-scale magnetic skyrmion lattice in two dimensions, Nat Phys 7, 713 (2011).
  * Romming  _et al._ [2013] N. Romming, C. Hanneken, M. Menzel, J. E. Bickel, B. Wolter, K. von Bergmann, A. Kubetzka, and R. Wiesendanger, Writing and deleting single magnetic skyrmions, Science 341, 636 (2013), https://www.science.org/doi/pdf/10.1126/science.1240573 .
  * Gross  _et al._ [2018] I. Gross, W. Akhtar, A. Hrabec, J. Sampaio, L. J. Martínez, S. Chouaieb, B. J. Shields, P. Maletinsky, A. Thiaville, S. Rohart, and V. Jacques, Skyrmion morphology in ultrathin magnetic films, Phys. Rev. Materials 2, 024406 (2018).
  * Seki  _et al._ [2012a] S. Seki, X. Z. Yu, S. Ishiwata, and Y. Tokura, Observation of skyrmions in a multiferroic material, Science 336, 198 (2012a), https://www.science.org/doi/pdf/10.1126/science.1214143 .
  * Yu  _et al._ [2012] X. Yu, M. Mostovoy, Y. Tokunaga, W. Zhang, K. Kimoto, Y. Matsui, Y. Kaneko, N. Nagaosa, and Y. Tokura, Magnetic stripes and skyrmions with helicity reversals, Proceedings of the National Academy of Sciences 109, 8856–8860 (2012).
  * Yu  _et al._ [2018] X. Z. Yu, W. Koshibae, Y. Tokunaga, K. Shibata, Y. Taguchi, N. Nagaosa, and Y. Tokura, Transformation between meron and skyrmion topological spin textures in a chiral magnet, Nature 564, 95 (2018).
  * Gao  _et al._ [2020] S. Gao, H. D. Rosales, F. A. G. Albarracín, V. Tsurkan, G. Kaur, T. Fennell, P. Steffens, M. Boehm, P. Čermák, A. Schneidewind, E. Ressouche, D. C. Cabra, C. Rüegg, and O. Zaharko, Fractional antiferromagnetic skyrmion lattice induced by anisotropic couplings, Nature 586, 37 (2020).
  * Divic  _et al._ [2021] S. Divic, H. Ling, T. Pereg-Barnea, and A. Paramekanti, Magnetic skyrmion crystal at a topological insulator surface (2021), arXiv:2103.15841 [cond-mat.mes-hall] .
  * Nomura and Nagaosa [2010] K. Nomura and N. Nagaosa, Electric charging of magnetic textures on the surface of a topological insulator, Phys. Rev. B 82, 161401 (2010).
  * Hurst  _et al._ [2015] H. M. Hurst, D. K. Efimkin, J. Zang, and V. Galitski, Charged skyrmions on the surface of a topological insulator, Phys. Rev. B 91, 060401 (2015).
  * Tiwari  _et al._ [2019] K. L. Tiwari, J. Lavoie, T. Pereg-Barnea, and W. A. Coish, Tunable skyrmion-skyrmion binding on the surface of a topological insulator, Phys. Rev. B 100, 125414 (2019).
  * Paul and Fu [2021] N. Paul and L. Fu, Topological magnetic textures in magnetic topological insulators, Physical Review Research 3, 10.1103/physrevresearch.3.033173 (2021).
  * Forsythe  _et al._ [2018] C. Forsythe, X. Zhou, K. Watanabe, T. Taniguchi, A. Pasupathy, P. Moon, M. Koshino, P. Kim, and C. R. Dean, Band structure engineering of 2d materials using patterned dielectric superlattices, Nature Nanotechnology 13, 566–571 (2018).
  * Li  _et al._ [2021] Y. Li, S. Dietrich, C. Forsythe, T. Taniguchi, K. Watanabe, P. Moon, and C. R. Dean, Anisotropic band flattening in graphene with one-dimensional superlattices, Nature Nanotechnology 16, 525 (2021).
  * Song  _et al._ [2010] C.-L. Song, Y.-L. Wang, Y.-P. Jiang, Y. Zhang, C.-Z. Chang, L. Wang, K. He, X. Chen, J.-F. Jia, Y. Wang, Z. Fang, X. Dai, X.-C. Xie, X.-L. Qi, S.-C. Zhang, Q.-K. Xue, and X. Ma, Topological insulator bi2se3 thin films grown on double-layer graphene by molecular beam epitaxy, Applied Physics Letters 97, 143118 (2010), https://doi.org/10.1063/1.3494595 .
  * Chou  _et al._ [2021] Y.-Z. Chou, J. Cano, and J. H. Pixley, Band manipulation and spin texture in interacting moiré helical edges, Physical Review B 104, 10.1103/physrevb.104.l201113 (2021).
  * Liu  _et al._ [2021] Z. Liu, H. Wang, and J. Wang, Magnetic moiré surface states and flat chern band in topological insulators (2021), arXiv:2106.01630 [cond-mat.mes-hall] .
  * Shi  _et al._ [2019] L.-k. Shi, J. Ma, and J. C. Song, Gate-tunable flat bands in van der waals patterned dielectric superlattices, 2D Materials 7, 015028 (2019).
  * Zhang  _et al._ [2014] J. Zhang, C. Triola, and E. Rossi, Proximity effect in graphene–topological-insulator heterostructures, Phys. Rev. Lett. 112, 096802 (2014).
  * Schouteden  _et al._ [2016] K. Schouteden, Z. Li, T. Chen, F. Song, B. Partoens, C. Van Haesendonck, and K. Park, Moiré superlattices at the topological insulator bi2te3, Scientific Reports 6, 20278 (2016).
  * Vargas  _et al._ [2017] A. Vargas, F. Liu, C. Lane, D. Rubin, I. Bilgin, Z. Hennighausen, M. DeCapua, A. Bansil, and S. Kar, Tunable and laser-reconfigurable 2d heterocrystals obtained by epitaxial stacking of crystallographically incommensurate bi2se3 and mos2 atomic layers, Science Advances 3, e1601741 (2017), https://www.science.org/doi/pdf/10.1126/sciadv.1601741 .
  * Hennighausen  _et al._ [2019a] Z. Hennighausen, C. Lane, A. Benabbas, K. Mendez, M. Eggenberger, P. M. Champion, J. T. Robinson, A. Bansil, and S. Kar, Oxygen-induced in situ manipulation of the interlayer coupling and exciton recombination in bi2se3/mos2 2d heterostructures, ACS Applied Materials & Interfaces 11, 15913 (2019a), https://doi.org/10.1021/acsami.9b02929 .
  * Hennighausen  _et al._ [2019b] Z. Hennighausen, C. Lane, I. G. Buda, V. K. Mathur, A. Bansil, and S. Kar, Evidence of a purely electronic two-dimensional lattice at the interface of tmd/bi2se3 heterostructures, Nanoscale 11, 10.1039/c9nr04412d (2019b).
  * Dunbrack and Cano [2021] A. Dunbrack and J. Cano, Magic angle conditions for twisted 3d topological insulators (2021), arXiv:2112.11464 [cond-mat.str-el] .
  * Tang and Fu [2014] E. Tang and L. Fu, Strain-induced partially flat band, helical snake states and interface superconductivity in topological crystalline insulators, Nature Physics 10, 964 (2014).
  * Liu  _et al._ [2014] Y. Liu, Y. Li, S. Rajput, D. Gilks, L. Lari, P. Galindo, M. Weinert, V. Lazarov, and L. Li, Tuning dirac states by strain in the topological insulator bi2se3, Nature Physics 10, 294 (2014).
  * Yuan and Fu [2020] N. F. Q. Yuan and L. Fu, Classification of critical points in energy bands based on topology, scaling, and symmetry, Phys. Rev. B 101, 125120 (2020).
  * Yuan  _et al._ [2019] N. F. Q. Yuan, H. Isobe, and L. Fu, Magic of high-order van hove singularity, Nature Communications 10, 10.1038/s41467-019-13670-9 (2019).
  * Chandrasekaran  _et al._ [2020] A. Chandrasekaran, A. Shtyk, J. J. Betouras, and C. Chamon, Catastrophe theory classification of fermi surface topological transitions in two dimensions, Phys. Rev. Research 2, 013355 (2020).
  * Guerci  _et al._ [2021] D. Guerci, P. Simon, and C. Mora, Higher-order van hove singularity in magic-angle twisted trilayer graphene (2021), arXiv:2106.14911 [cond-mat.mes-hall] .
  * Mudry [2014] C. Mudry,  _Lecture Notes on Field Theory in Condensed Matter Physics_ (WORLD SCIENTIFIC, 2014) https://www.worldscientific.com/doi/pdf/10.1142/8697 .
  * Fradkin [2013] E. Fradkin,  _Field Theories of Condensed Matter Physics_, 2nd ed. (Cambridge University Press, 2013).
  * Altland and Simons [2010] A. Altland and B. D. Simons,  _Condensed Matter Field Theory_, 2nd ed. (Cambridge University Press, 2010).
  * Stoner [1938] E. C. Stoner, Collective electron ferromagnetism, Proceedings of the Royal Society of London. Series A. Mathematical and Physical Sciences 165, 372 (1938), https://royalsocietypublishing.org/doi/pdf/10.1098/rspa.1938.0066 .
  * Litak  _et al._ [1998] G. Litak, A. Martin, B. Györffy, J. Annett, and K. Wysokiński, Van hove singularity and d-wave pairing in disordered superconductors, Physica C: Superconductivity 309, 257–262 (1998).
  * Wilson  _et al._ [2020] J. H. Wilson, Y. Fu, S. Das Sarma, and J. H. Pixley, Disorder in twisted bilayer graphene, Physical Review Research 2, 10.1103/physrevresearch.2.023325 (2020).
  * Chandrasekaran and Betouras [2022] A. Chandrasekaran and J. J. Betouras, Effect of disorder on density of states and conductivity in higher-order van hove singularities in two-dimensional bands, Phys. Rev. B 105, 075144 (2022).
  * Sawada and Nagaosa [2005] K. Sawada and N. Nagaosa, Optical magnetoelectric effect in multiferroic materials: Evidence for a lorentz force acting on a ray of light, Phys. Rev. Lett. 95, 237402 (2005).
  * Ederer and Spaldin [2007] C. Ederer and N. A. Spaldin, Towards a microscopic theory of toroidal moments in bulk periodic crystals, Phys. Rev. B 76, 214404 (2007).
  * Bhowal and Spaldin [2022] S. Bhowal and N. A. Spaldin, Magnetoelectric classification of skyrmions (2022), arXiv:2201.01667 [cond-mat.mes-hall] .
  * Zhong  _et al._ [2016] S. Zhong, J. E. Moore, and I. Souza, Gyrotropic magnetic effect and the magnetic moment on the fermi surface, Physical Review Letters 116, 10.1103/physrevlett.116.077201 (2016).
  * Ma and Pesin [2015] J. Ma and D. A. Pesin, Chiral magnetic effect and natural optical activity in metals with or without weyl points, Phys. Rev. B 92, 235205 (2015).
  * Yu  _et al._ [2010] X. Yu, Y. Onose, N. Kanazawa, J. Park, J. Han, Y. Matsui, N. Nagaosa, and Y. Tokura, Real-space observation of a two-dimensional skyrmion crystal, Nature 465, 901 (2010).
  * Dovzhenko  _et al._ [2018] Y. Dovzhenko, F. Casola, S. Schlotter, T. X. Zhou, F. Büttner, R. L. Walsworth, G. S. D. Beach, and A. Yacoby, Magnetostatic twists in room-temperature skyrmions explored by nitrogen-vacancy center spin texture reconstruction, Nature Communications 9, 10.1038/s41467-018-05158-9 (2018).
  * Grigoriev  _et al._ [2007] S. V. Grigoriev, V. A. Dyadkin, D. Menzel, J. Schoenes, Y. O. Chetverikov, A. I. Okorokov, H. Eckerlebe, and S. V. Maleyev, Magnetic structure of ${\mathrm{fe}}_{1-x}{\mathrm{co}}_{x}\mathrm{Si}$ in a magnetic field studied via small-angle polarized neutron diffraction, Phys. Rev. B 76, 224424 (2007).
  * Seki  _et al._ [2012b] S. Seki, J.-H. Kim, D. S. Inosov, R. Georgii, B. Keimer, S. Ishiwata, and Y. Tokura, Formation and rotation of skyrmion crystal in the chiral-lattice insulator cu2oseo3, Phys. Rev. B 85, 220406 (2012b).
  * Grigoriev  _et al._ [2009] S. V. Grigoriev, D. Chernyshov, V. A. Dyadkin, V. Dmitriev, S. V. Maleyev, E. V. Moskvin, D. Menzel, J. Schoenes, and H. Eckerlebe, Crystal handedness and spin helix chirality in ${\mathrm{fe}}_{1-x}{\mathrm{co}}_{x}\mathrm{Si}$, Phys. Rev. Lett. 102, 037204 (2009).
  * Wilhelm  _et al._ [2011] H. Wilhelm, M. Baenitz, M. Schmidt, U. K. Rößler, A. A. Leonov, and A. N. Bogdanov, Precursor phenomena at the magnetic ordering of the cubic helimagnet fege, Phys. Rev. Lett. 107, 127203 (2011).
  * Hirschberger  _et al._ [2019] M. Hirschberger, T. Nakajima, S. Gao, L. Peng, A. Kikkawa, T. Kurumaji, M. Kriener, Y. Yamasaki, H. Sagayama, H. Nakao, K. Ohishi, K. Kakurai, Y. Taguchi, X. Yu, T. hisa Arima, and Y. Tokura, Skyrmion phase and competing magnetic orders on a breathing kagomé lattice, Nature Communications 10 (2019).
  * Brearton  _et al._ [2021] R. Brearton, L. Turnbull, J. Verezhak, G. Balakrishnan, P. Hatton, G. Laan, and T. Hesjedal, Deriving the skyrmion hall angle from skyrmion lattice dynamics, Nature Communications 12, 2723 (2021).
  * Huang  _et al._ [2017] B. Huang, G. Clark, E. Navarro-Moratalla, D. R. Klein, R. Cheng, K. L. Seyler, D. Zhong, E. Schmidgall, M. A. McGuire, D. H. Cobden, W. Yao, D. Xiao, P. Jarillo-Herrero, and X. Xu, Layer-dependent ferromagnetism in a van der waals crystal down to the monolayer limit, Nature 546, 270–273 (2017).
  * Liu  _et al._ [2020] J. Liu, A. Singh, B. Kuerbanjiang, C. H. W. Barnes, and T. Hesjedal, Kerr effect anomaly in magnetic topological insulator superlattices, Nanotechnology 31, 434001 (2020).
  * Song  _et al._ [2018] T. Song, X. Cai, M. W.-Y. Tu, X. Zhang, B. Huang, N. P. Wilson, K. L. Seyler, L. Zhu, T. Taniguchi, K. Watanabe, M. A. McGuire, D. H. Cobden, D. Xiao, W. Yao, and X. Xu, Giant tunneling magnetoresistance in spin-filter van der waals heterostructures, Science 360, 1214–1218 (2018).
  * Ye  _et al._ [1999] J. Ye, Y. B. Kim, A. J. Millis, B. I. Shraiman, P. Majumdar, and Z. Tešanović, Berry phase theory of the anomalous hall effect: Application to colossal magnetoresistance manganites, Phys. Rev. Lett. 83, 3737 (1999).
  * Taguchi  _et al._ [2001] Y. Taguchi, Y. Oohara, H. Yoshizawa, N. Nagaosa, and Y. Tokura, Spin chirality, berry phase, and anomalous hall effect in a frustrated ferromagnet, Science 291, 2573 (2001), https://www.science.org/doi/pdf/10.1126/science.1058161 .
  * Bruno  _et al._ [2004] P. Bruno, V. K. Dugaev, and M. Taillefumier, Topological hall effect and berry phase in magnetic nanostructures, Phys. Rev. Lett. 93, 096806 (2004).
  * Onoda  _et al._ [2004] M. Onoda, G. Tatara, and N. Nagaosa, Anomalous hall effect and skyrmion number in real and momentum spaces, Journal of the Physical Society of Japan 73, 2624 (2004).
  * Binz and Vishwanath [2008] B. Binz and A. Vishwanath, Chirality induced anomalous-hall effect in helical spin crystals, Physica B: Condensed Matter 403, 1336 (2008).
  * Nagaosa and Tokura [2012] N. Nagaosa and Y. Tokura, Emergent electromagnetism in solids, Physica Scripta T146, 014020 (2012).
  * Schulz  _et al._ [2012] T. Schulz, R. Ritz, A. Bauer, M. Halder, M. Wagner, C. Franz, C. Pfleiderer, K. Everschor, M. Garst, and A. Rosch, Emergent electrodynamics of skyrmions in a chiral magnet, Nature Physics 8, 301 (2012).
  * Shimizu  _et al._ [2021] K. Shimizu, S. Okumura, Y. Kato, and Y. Motome, Spin moiré engineering of topological magnetism and emergent electromagnetic fields, Physical Review B 103, 10.1103/physrevb.103.184421 (2021).
  * Guan  _et al._ [2021] Y. Guan, O. V. Yazyev, and A. Kruchkov, Unconventional flat chern bands and 2e charges in skyrmionic moiré superlattices (2021), arXiv:2111.10410 [cond-mat.str-el] .
  * Wang  _et al._ [2021] J. Wang, J. Zang, J. Cano, and A. J. Millis, Staggered Pseudo Magnetic Field in Twisted Transition Metal Dichalcogenides: Physical Origin and Experimental Consequences, arXiv e-prints , arXiv:2110.14570 (2021), arXiv:2110.14570 [cond-mat.mes-hall] .
  * Zang  _et al._ [2021] J. Zang, J. Wang, J. Cano, and A. J. Millis, Hartree-fock study of the moiré hubbard model for twisted bilayer transition metal dichalcogenides, Phys. Rev. B 104, 075150 (2021).
  * Fert  _et al._ [2013] A. Fert, V. Cros, and J. Sampaio, Skyrmions on the track, Nature Nanotechnology 8, 152–156 (2013).
  * Fert  _et al._ [2017] A. Fert, N. Reyren, and V. Cros, Magnetic skyrmions: advances in physics and potential applications, Nature Reviews Materials 2, 10.1038/natrevmats.2017.31 (2017).
  * Zhang  _et al._ [2015] X. Zhang, M. Ezawa, and Y. Zhou, Magnetic skyrmion logic gates: conversion, duplication and merging of skyrmions, Scientific Reports 5, 10.1038/srep09400 (2015).
  * Zhang  _et al._ [2020] X. Zhang, Y. Zhou, K. Mee Song, T.-E. Park, J. Xia, M. Ezawa, X. Liu, W. Zhao, G. Zhao, and S. Woo, Skyrmion-electronics: writing, deleting, reading and processing magnetic skyrmions toward spintronic applications, Journal of Physics: Condensed Matter 32, 143001 (2020).
  * Bi  _et al._ [2019] Z. Bi, N. F. Q. Yuan, and L. Fu, Designing flat bands by strain, Physical Review B 100, 10.1103/physrevb.100.035448 (2019).
  * Balents [2019] L. Balents, General continuum model for twisted bilayer graphene and arbitrary smooth deformations, SciPost Physics 7, 10.21468/scipostphys.7.4.048 (2019).
  * He  _et al._ [2020] W.-Y. He, D. Goldhaber-Gordon, and K. T. Law, Giant orbital magnetoelectric effect and current-induced magnetization switching in twisted bilayer graphene, Nature Communications 11, 10.1038/s41467-020-15473-9 (2020).
  * Nam and Koshino [2017] N. N. T. Nam and M. Koshino, Lattice relaxation and energy band modulation in twisted bilayer graphene, Phys. Rev. B 96, 075311 (2017).