# A hierarchy of multi-moment skyrmion phases in twisted magnetic bilayers

Sujay Ray, and Tanmoy Das<EMAIL_ADDRESS>Department of Physics, Indian
Institute of Science, Bangalore-560012, India

###### Abstract

The recent discovery of two-dimensional (2D) Van der Waals (VdW) magnets is a
crucial turning point in the quantum magnet research field, since quantum
fluctuations and experimental difficulties often elude stable magnetic orders
in 2D. This opens new doors to delve for novel quantum and topological spin
configurations, which may or may not have direct analogs in bulk counterparts.
Here we study a twisted bilayer geometry of 2D magnets in which long-range
spin-spin interactions naturally commence along the inter-layer Heisenberg
($J_{\perp}$) and dipole-dipole ($J_{\rm D}$) channels. The $J_{\perp}-J_{\rm
D}$ parameter space unveils a hierarchy of distinct skyrmions phases, ranging
from point-, rod-, and ring-shaped topological charge distributions.
Furthermore, we predict a novel topological antiferroelectric phase, where
oppositely-charged antiskyrmion pairs are formed, and the corresponding
topological dipole moments become ordered in a Néel-like state $-$ hence
dubbed topological antiferroelectric state. The results indicate that twisted
magnetic layer provides a new setting to engineer and tune a plethora of novel
and exotic skyrmion phases and their dynamics.

##  I Introduction

Skyrmion is a particle dual to a topological configuration of the $O(3)$
fields (read spin) in a 2+1 dimensionSkyrme ; SkyrmionGenealTheory ;
SkyrmionReview1 ; Skyrmionbook . Such a spin configuration is an allowed
classical solution of the non-linear sigma model. However, the main challenge
lies in stabilizing a skyrmion solution at a saddle-point energy minimum,
requiring distinct magnetic interactions and frustration. Long-range dipole-
dipole interaction, in addition to easy-plane magnetic anisotropy and magnetic
field, was initially proposed to mediate skyrmion solution.SkyrmionDDI
Dzyaloshinskii-Moriya interaction (DMI) brings in the chiral spin-spin
interaction required for a skyrmion solutionSkyrmionGenealTheory ; SkyrmionDMI
; skDMI1 ; skDMI2 ; skDMI3 ; skDMI4 . Competition between DMI and
ferromagnetic exchange interactions makes a spatially varying magnetic texture
to have a lower energy than the mean-field long-range magnetic order.
Geometrical and magnetic frustration can also stabilize skyrmion
structuresmagnetic_frus1 ; magnetic_frus2 ; magnetic_frus3 ; SkyrmionFrus .
This is induced, for example, in a triangular lattice by the competition
between a ferromagnetic nearest neighbour (NN) exchange with an
antiferromagnetic next nearest neighbour (NNN) interaction. Apart from these,
more recently, spin-orbit couplingSkyrmionSOC , Kondo couplingmagnetic_frus4 ,
and magnetic disorder with the application of magnetic pulse skyrmion_disorder
are shown to assist skyrmion solution. Proposals to obtain skyrmions via
quantum Hall substratesSkyrmionQHI  and optical latticesSkyrmionBEC  are also
presented.

By now, there have been several materials realizations of the skyrmions,
mainly in systems with DMI, such as Bloch-type skyrmions in MnSiMnSi ,
Co0.5Fe0.5Si,FeCoSi ; Tokura10 ; MildeKohlerSci13  Cu2OSeO3,Cu2OSeO3
CoZnMn,CoZnMn  and FeGeFeGe , and Néel-type skyrmions in ferromagnetic heavy-
metalsheavymetals , Kagome lattice Fe3Sn2Kagome , magnetic
films,WisendangerReview , and antiskyrmions in Mn2RhSnMn2RhSnParkin ;
Antiskyrmion . Experimental abilities to breed,MildeKohlerSci13 ;
WritingSkyrmion  mobilize,MobileSkyrmion  rotateMochizukiNatMat  skyrmions are
also demonstrated recently.Hofman ; WisendangerReview  Owing to the
topological robustness, skyrmions have numerous potential applications in
quantum information,SkyrmionReviewTech  racetrack memoryRacetrackParkin ;
Racetrack2 , which demand enhanced materials flexibilities and tunabilities.

Obtaining skyrmion solutions in 2D systems will give new opportunities for
science and applications. Recently, long-range magnetic order has been
observed in 2D Van der Waals (VdW) chalcogenides, halides, and related
materials. Intrinsic antiferromagnetic order is observed in monolayer FePS3,
FePS3 ; FePS3a , and in $M$P$X_{3}$($M$=Mn,Fe,Co,Ni; $X$=S,Se)MPX ; MPXa .
Later on, many VdW materials such as Fe3GeTe2FGT ; FGTa ,
Mn$X_{2}$($X$=S,Se)MnX ; MnXa , V$X_{2}$($X$=S,Se,Te)VX ; VXa ; VXb  were
found to be intrinsic ferromagnet. Another exciting family of 2D vdW magnets
is the Cr based materials Cr$X_{3}$ ($X=$I,Br,Cl)CrX ; CrXa , which are
ferromagnets in monolayer, but antiferromagnets in bilayer structure, and the
two orderings are externally tunable.CrXbi ; CrXbia  Theoretical and
experimental efforts to obtain skyrmions in 2D systems are present. Continuum
theory of magnetization field in the non-linear sigma model and Landau-
Lifshitz-Gilbert in monolayer and bilayer Moiré systems show the existence of
Neél type skyrmions.CrI3_DFT ; Moire_magnets ; skTBG ; Substrate ; skFerroDMI
Recently, a skyrmion phase is observed in 2D Fe3GeTe2 on (Co/Pd)n
superlatticeFGT_CoPd  and Fe3GeTe2/$h$-BN heterostructureFGT_hBN  due to their
sizable DMI strength.

To strategize a new mechanism of the skyrmion, it is worth revisiting it’s key
ingredients. Firstly, topological skyrmion configurations in 2+1 dimension
generally belong to the homotopy group
$\pi_{2}(\mathbb{S}^{2})\cong\mathbb{Z}$. The homotopy mapping is exact when
both the coordinate space and the target (spin) space are compact
$\mathbb{S}^{2}$. The constraint $|{\bf S}|=1$ compactifies the spin space.
When a spin configuration has a one-to-one correspondence with the spatial
dimensions, this in turn compactifies the position space
$\mathbb{R}^{2}\rightarrow\mathbb{S}^{2}$. The resulting one-to-one mapping
guarantees the spin configuration to be topological with its skyrmion charge
$Q\in\mathbb{Z}-$ an integer winding number. We can reverse the above
reasoning for a bottom-up approach. If the effective magnetic field ${\bf
B}({\bf r})$, experienced by a local spin due to the surrounding spins and
extrinsic fields, lives on a Bloch sphere $\mathbb{S}^{2}$, then within the
minimal Zeeman-like coupling, the field would lay the ground for a topological
configuration for the spins. The second essential requirement is that the
local field configuration must concomitantly promote a saddle-point energy
minimum to stabilize a skyrmion structure.

Figure 1: Twisted bilayer setup and five distinct topological spin
configurations in real space. (a) Schematic diagram of twisted bilayer
honeycomb ($\theta$ = twist angle) lattice with a Moiré pattern is shown along
with two lattice translation vectors ${\bf{R}}_{1,2}$. A 3D illustration is
shown on the lower right corner, in which a single spin of the upper layer
(bright magenta color) interacts with bottom-layer spins within a circular
region (dark green color) of radius rcut. (b-f) Five distinct spin textures
(enumerated by Phase-I to Phase-V) obtained in the $J_{\perp}-J_{\rm D}$
parameter spare (see Fig. 2). Arrows denote the spin ${\bf s}$ vector, while
red to green color-gradient denotes $s_{z}=1$ to $-1$ values. Dashed thin line
in each panel indicates the direction along which a 2D spin projection is
shown in the corresponding inset. (b),(d),(f) We show three skyrmion
structures with distinct charge density (shown in Fig. 3) and integer
topological charge $Q$. (c) Phase-II corresponds to a novel higher-order
topological phase with streamlines of down spins, and topological dipole
moments of antiskyrmions pairs (see Fig. 4). (e) This is a trivial topological
phase with finite noncollinear ferromagnetic moment.

Guided by these principles, and with the recent discoveries of 2D magnets, we
lay a blueprint for novel and multifaceted skyrmions (and antiskyrmions) in
twisted magnetic bilayers. We construct a Moiré superlattice of spins formed
in twisted bilayer of VdW magnetic layers with ferromagnetic order at the
bottom and $O(3)$ spin dynamics on the top layer. The setup is illustrated in
Fig. 1(a). We include Heisenberg exchange terms $J_{||}$ (for intra-layer) and
$J_{\perp}$ (for inter-layer) interactions, and the inter-layer dipole
interaction $J_{\rm D}$ as shown in DFT calculations to be dominant in such
setup.CrI3_DFT  We carry out Monte Carlo simulations to determine the
microscopic ground state spin configurations at low-temperature, and sweep the
entire $J_{\perp}/J_{||}$ and $J_{D}/J_{||}$ parameter space. We identify
three distinct skyrmion phases with topological charges $Q=\pm 1$. whose
topological charge distributions reveal a previously unknown hierarchy of
point-, rod-, and circle - shape in different topological phases. We predict a
novel topological spin configuration in the vicinity of the $J_{\perp}/J_{\rm
D}\sim-0.4$. We find that near the Moiré lattice sites, a pair of spatially
separated and oppositely charged antiskyrmions is formed and govern a
topological (‘electric’) dipole moment. More interestingly, such dipoles are
found to align anti-ferroelectrically between the nearest-neighbor sites of
the Moiré lattice, and produce a Néel like order for the topological electric
dipole moment. We explain these results with a dual electromagnetic theory,
demonstrating the ‘electric field’ lines for all topological charge
distributions. We also study the ‘x-ray diffraction’ (XRD) pattern of the
topological charges (topological charge-charge correlation function) as well
as the spin-spin correlation functions to elucidate the crystallization and
phase transitions of the topological charge centers and dipoles. These results
expand the list of possible skyrmion and magnetic phases (Néel and Bloch
phase) obtained in continuum models to a hierarchy of skyrmions and its higher
order topological phase.continum_model ; Moire_magnets

Figure 2: Phase diagram of five distinct topological and quantum phases. The
inter-layer exchange coupling $J_{\perp}$ and dipole interaction $J_{\rm D}$
are varied with respect to the in-plane exchange coupling $J_{||}$. $J_{\rm
D}<0$ is an unphysical value, however, mathematically it leaves the phase
diagram symmetric when the sign of $J_{\perp}$ is also reversed. Red and blue
color distinguishes skyrmion charges of $Q=\pm 1$. Phase-II, although features
a net $Q=0$, exhibits a novel topological dipole moment of antiskyrmions which
produces an antiferroelectric order state. The horizontal dashed line
indicates the realistic parameter range as deduced in a DFT
calculation.CrI3_DFT

##  II Realization

We consider a single layer honeycomb magnet (magenta color) placed on a single
layer magnetic substrate (green color) of the same lattice structure and
lattice constant, as shown in Fig. 1(a). The distance between the layers $d$
is taken to be same as the lattice constant of the honeycomb lattice. The spin
in the substrate layer is fixed to be a collinear ferromagnetic state. This
can be achieved with a strong bulk ferromagnetic material as studied in the
literatureSubstrate . We primarily focus on small relative twist angles
$\theta$ which give the commensurate Moiré superlattices. We focus on
$\theta=1.61^{o}$, which gives a hexagonal Moiré lattice with $a=35.6a_{0}$,
where $a_{0}$ is the lattice constant of the single layer system. A critical
number of atoms in a Moiré supercell, determined by the twist angle, is
important to stabilize a skyrmion. Above this critical value, the obtained
topological phase diagram remains essentially invariant to twist angles and
number of atoms, except that the skyrmion radius grows with the Moiré cell
dimension. So, there is an upper critical value of the twist angle ($\sim
2^{0}$) above which the Moiré unit cell becomes small enough and the magnetic
unit cell is no longer commensurate with the Moiré unit cell.

The general computational strategy is as follows. We split the full
Hamiltonian into two parts: $H=H_{1}+H_{2}$, where $H_{1}$ and $H_{2}$ are the
intra-layer, and inter-layer parts. The intra-layer Hamiltonian consists of a
nearest-neighbor Heisenberg exchange and a spin asymmetry term, for both
layers. The inter-layer term consists of many neighbors Heisenberg exchange
term and the dipole-dipole interaction term. We then integrate out the bottom
layer’s spin to obtain an effective Hamiltonian for the top layer as $H_{\rm
top}\sim H_{1}+{\bf B}\cdot{\bf s}$, where ${\bf B}$ is the effective magnetic
field exerted from the bottom layer on the top layers spin ${\bf s}$. The
Moiré periodicity is imposed by expanding these terms in the plane wave basis
of the Moiré supercell. We solve $H_{\rm top}$ within the Monte Carlo
simulation.

We now give the details of the model. We label the spin variables for the top
and bottom (subtrate) layers by ${\bf s}$, and ${\bf S}$, respectively. The 2D
VdW systems in single and bilayer setups are observed to show an in-plane
ferromagnetic (and out-of-plane ferro- or anti-ferromagnetic) order with the
spin quantization axis to be out-of-plane ($z$-direction).CrI3 ; Chalcogenides
Such a magnetic ground state is reproduced by the model:

$\displaystyle H_{1}({\bf s})=-J_{||}\sum_{<ij>}{\bf s}_{i}.{\bf
s}_{j}-K\sum_{i}\left(s_{iz}\right)^{2},$ (1)

where $i,j$ lattice sites within a Moiré supercell. The first term is the
nearest neighbor Heisenberg interaction with coupling constant $J_{||}>0$ for
a ferromagnetic phase, and $K$ gives the $z$-axis asymmetry, breaking the
$O(3)$ spin degeneracy. $H_{1}({\bf S})$ term gives the corresponding
Hamiltonian for the bottom layer, with $J_{||}$ and $K$ kept fixed.

The inter-layer interaction $H_{2}$ is the crucial part. Depending on the
twist angle, especially at small twist angles, a spin at the top layer
interacts with several neighboring bottom layer spins, and hence the inter-
layer interaction involves terms beyond the nearest-neighbor exchange
interaction. The inter-layer interaction is mainly dominated by several
nearest-neighbor exchange interaction $H_{\rm ex}$ and dipole-dipole
interaction ($H_{\rm D}$) as $H_{2}({\bf S},{\bf s})=H_{\rm ex}+H_{\rm D}$
where

$\displaystyle H_{\rm ex}$ $\displaystyle=$
$\displaystyle-\sum_{ij}J_{\perp}({\bf r}_{ij}){\bf S}_{i}.{\bf s}_{j},$ (2)
$\displaystyle H_{\rm D}$ $\displaystyle=$ $\displaystyle J_{\rm
D}\sum_{ij}\frac{1}{r_{ij}^{3}}\left[{\bf S}_{i}.{\bf s}_{j}-3({\bf
S}_{i}.\hat{{\bf r}}_{ij}).({\bf s}_{j}.\hat{{\bf r}}_{ij})\right].$ (3)

$i$, $j$ are the site indices in the bottom and top layer, respectively. ${\bf
r}_{ij}={\bf r}_{i}-{\bf r}_{j}$ with $\hat{{\bf r}}_{ij}$ being the
corresponding unit vector. Since both interactions are long-ranged, we need to
set a cutoff radius $r_{\rm cut}$, Fig. 1(a) (inset). Due to higher power of
$r_{ij}$ in the denominator, the result converges quickly by $r_{\rm
cut}<20a_{0}$ which is much smaller than the Moiré latiice size $\sim
35a_{0}$.

Next, we integrate out the bottom layer’s spins ${\bf S}$, and define an
effective magnetic field at the top layer at ${\bf r}_{i}$ as ${\bf B}({\bf
r}_{i})={\bf B}_{\rm ex}({\bf r}_{i})+{\bf B}_{\rm D}({\bf r}_{i})$, where
${\bf B}_{\rm ex}$, and ${\bf B}_{\rm D}$ distinguish the contributions from
the exchange and dipole-dipole interaction terms as ${\bf B}_{\rm ex}({\bf
r}_{i})=-\frac{\partial H_{\rm ex}}{\partial{\bf
s}_{i}}=\frac{J_{\perp}}{2}\sum_{a}e^{i{\bf G}_{a}.{\bf r}_{i}}$, and ${\bf
B}_{\rm D}({\bf r}_{i})=-\frac{\partial H_{\rm D}}{\partial{\bf s}_{i}}$. We
set ${\bf S}_{j}={\hat{z}}$ for all unit spins at the bottom layer. ${\bf
G}_{a}(a=1$ \- $6)$ are the six minimal reciprocal lattice vectors of the
Moiré superlattice. We include the six possible smallest reciprocal lattice
vectors in the Moiré lattice and the results do not change with the inclusion
of negligibly small contribution of higher reciprocal lattice vectors. Dipole
interaction is known to be an useful ingredient for (generally bubble-type)
skyrmions and antiskyrmions,SkyrmionDDI ; SkyrmionReview1 ; Antiskyrmion  but
are significantly weaker in strength in real materials. In twisted bilayer
system, however, the intra-layer dipole interaction is considerably enhanced,
and is found here to be detrimental to the bubble or Bloch skyrmion phases
(see below), while promoting novel and distinct skyrmion phases.

In this algorithm, the bottom layer’s effects can be cast into a Zeeman-like
term at the top layer. The full Hamiltonian hence takes the form

$H=H_{1}-\sum_{i}{\bf B}({\bf r}_{i})\cdot{\bf s}_{i}.$ (4)

In this way, the total Hamiltonian simplifies to the form used in the
introduction for discussing how the in-plane spin is enforced to lie on a
topological compact space $\mathbb{S}^{2}$ by tailoring the ‘long-range’
Zeeman coupling. In our Monte Carlo simulation we minimize the local
Hamiltonian for a single spin with temperature annealing as well as parameter
annealing, and check for the total energy convergence of the system in each
case. In all our calculations, the minimization of the local Hamiltonian
corresponds to minimization of the total energy. Details of our Monte Carlo
simulation can be found in the Appendix B.

Our model and parameters are justified as follows. Note that the dipole-dipole
interaction $J_{\rm{D}}$ is only included in the inter-layer Hamiltonian. This
is because, in VdW magnets, DFT calculation has indicated that the intra-layer
exchange is almost 10 times larger than the inter-layer exchangeCrI3_DFT .
Therefore, any small amount of inter-layer dipole-dipole interaction can have
considerable effect on the spin configuration, whereas the intra-layer dipole
interaction is always subsided by the large exchange interaction. Furthermore,
the DMI in the existing 2D VdW magnets is small compared to strong intra-layer
exchange interaction. Therefore, realistic small values of the intra-layer
dipole interaction and DMI interaction do not impact our phase diagram (see
Discussion section). The realistic parameter range in CrI3 is indicated by
horizontal dash line in Fig. Fig. 2, as deduced from a DFT calculationCrI3_DFT
. Interestingly, the phase diagram essentially remains the same below and
above this parameter range, emphasizing that all the phases are attainable in
the existing materials.

##  III Results

Figure 3: Local variation of polarity, vorticity, and topological charge
densities in the five different topological phases. (a-e) The $z$-component of
the polarity density $n_{z}({\bf r})$ (see main text) is plotted in blue to
red colormap, denoting down to up spin components of $s_{z}$. (f-j) We plot
the vorticity density $v_{z}({\bf r})$ with blue to red colormap. The black
solid line marks the $n_{z}=0$ domain wall boundary. (k-o) Corresponding
topological charge density $\rho({\bf r})$ is plotted here. In Phase-II, IV,
although charge centers are formed but the net charge $Q=0$. In Phase-II,
fractional charge centers are confined by streamline domain wall, giving a
novel topological dipole moment. There is no domain wall of $n_{z}$ in Phase-
IV and hence its a trivial phase. (Different color scales are for the five
phases in the middle and bottom panels.)

###  III.1 Phase diagram

Our simulation yields a plethora of quantum and topological phases; five such
distinct configurations are identified in Fig. 1(b)-(f). We do not however
find any homogeneous mean-field solutions at finite twist angles, while at
zero twist angles we only find either a ferromagnetic or an antiferromagnetic
order depending on the sign of $J_{||}$. The corresponding phase diagram is
presented in Fig. 2 for $J_{\perp}/J_{||}$ and $J_{\rm D}/J_{||}$ with
$J_{||}>0$. In the phase diagram, the red and blue shaded areas denote
distinct skyrmion phases with topological invariant $Q=\pm 1$, respectively.
The white regions represent a higher-order dipoler antiskyrmion phase (Phase-
II), and a trivial phase (Phase-IV). From Fig. 1(b),(d), and (f), it is
evident that the three skyrmion phases, denoted by Phase-I, III, and V, are
characteristically exclusive, which will be distinguished below in multiple
ways. Phase-II has zero net topological charge, but possess fractional charge
centers of opposite sign, and hence a topological dipole moment which are
antiferroelectrically ordered. The Phase-IV bears no topological or exotic
quantum order (except a finite magnetization due to a collinear spin
ordering).

The particle dual of the skyrmion is a topological charge denoted by $Q=\int
d^{2}{\bf r}\rho({\bf r})$, where $\rho({\bf r})$ is the topological charge
density. With its corresponding current density $\mathcal{J}_{x,y}$, we
concisely define the three-component density operators
$\mathcal{J}_{\mu}=(\rho,\mathcal{J}_{x},\mathcal{J}_{y})$ as

$\displaystyle\mathcal{J}_{\mu}({\bf
r})=\frac{\epsilon_{\mu\nu\tau}}{8\pi}{\bf n}\cdot\partial_{\nu}{\bf
n}\times\partial_{\tau}{\bf n},$ (5)

where $\mu,\nu,\tau=0,x,y$ are time-space indices, and the ${\bf r}$
dependence on the unit vector field ${\bf n}={\bf s}/|{\bf s}|$ is implied. We
define the vortex density as ${\bf v}({\bf
r})=\epsilon_{\nu\tau}\partial_{\nu}{\bf n}\times\partial_{\tau}{\bf n}$. In
our layered geometry and with the $z$-axis asymmetry, it is natural to expect
that that the vorticity of the spin-texture commences in the $xy$-plane, i.e,
$v_{z}$ dominates. Then the corresponding polarity density is simply governed
by $n_{z}({\bf r})$. The $z$-components of the polarity density $n_{z}({\bf
r})$, the vortex density $v_{z}({\bf r})$ and the charge density $\rho({\bf
r})$ are investigated in Fig. 3 in three different rows for the five distinct
phases (different columns).

The mechanism of skyrmions and antiskyrmions is retrieved as follows. It is
known from the topological band theoryDasRMP  that the polarity field
(equivalent to the Dirac mass for fermion fields) forms a nodal (closed)
contour, across which $n_{z}({\bf r})$ changes sign $-$ this is called the
domain wall (see top row in Fig. 3). The vorticity field $v_{z}({\bf r})$
inside the domain wall acquires singularity $-$ either point- or rod-, or
ring- shaped singularity $-$ and cannot be removed with a trivial gauge
transformation (see middle row in Fig. 3). Hence the topological charge
density $\rho({\bf r})$ becomes confined within the domain wall (see bottom
row in Fig. 3). The homotopy mapping of the ${\bf n}({\bf r})$ field on
$\mathbb{S}^{2}$ in the ${\bf r}$-space quantizes the topological charge
$Q\in\mathbb{Z}$, where the integration of ${\bf r}$ is performed within a
single domain wall of the polarity density [see black solid lines in Fig.
3(f-o)].

###  III.2 Skyrmion hierarchy

Phase-I is a topological phase with $Q=-1$, and is present in most of the
$J_{\perp}/J_{||}<0$ region. It is destabilized at a small value of
$J_{\perp}$ by stronger dipole interaction $J_{\rm D}$. The distributions of
$n_{z}$, $v_{z}$, and $\rho$ for Phase-I are shown in the left-most column in
Fig. 3. The polarity density map $n_{z}({\bf r})$ demarcates a sharp and
circular domain wall boundary, which reminds us of a magnetic bubble observed
in astronomical space, as well in magnetic systems.MagBubble  The vorticity
and charge densities of this phase, however, reveal much richer structures
unknown before. In Fig. 3(f), we find that the nodal ring of the polarity
density (black line) encloses a circular vortex density $v_{z}({\bf r}$)
structure. In fact, $v_{z}({\bf r})$ is positive (negative) outside (inside)
the domain wall, and share the same nodal ring as that of $n_{z}({\bf r})$. In
what follows, the charge density $\rho({\bf r})$ also acquires a singular ring
geometry, confined by the domain wall boundary, see Fig. 3(k). This phase is
also topologically equivalent to the Dirac nodal ring stateDiracNodal  in the
electronic structure in which the topology is defined via Berry gauge
connection.

As we switch $J_{\perp}\rightarrow-J_{\perp}$, keeping all other interactions
fixed, we obtain a characteristically distinct skyrmion texture (denoted by
Phase-V) with opposite charge ${\it Q}$. Unlike the sharp domain wall in
Phase-I, $n_{z}$ varies smoothly with ${\bf r}$, and forms a (nearly)
elliptical domain wall in Phase-V. Moreover, the Phase-V has a point-like
topological charge center sitting at the Moiré supercell center. Hence as
opposed to ring-singularity in Phase-I, Phase-V acquires point- singularity.
The asymmetry between the Phase-I and Phase-V at $\pm J_{\perp}$ for fixed
$J_{\rm D}$, results from the competition between $J_{\perp}$ and $J_{\rm D}$.
The phase diagram is reversal between $\pm J_{\perp}$ for $J_{\rm
D}\rightarrow-J_{\rm D}$.

The skyrmion Phase-III occurs in the vicinity of $J_{\perp}\sim 0$, and is
mainly stabilized by the long-range (out-of-plane) dipole-dipole interaction
$J_{\rm D}$. The magnetic domains are elliptical in shape, and concentrated at
two sides of the Moiré supercell. As seen in Fig. 3(h) and (m), the $n_{z}$
nodal contour confines a fixed-sign vorticity field (positive), and hence the
topological charge distribution (negative since $n_{z}<0$) becomes quantized.
The topological charge density is distributed inside the elliptical domain
wall and gives a rod-line singularity. Such a rod-like topological charge
distribution repeats periodically [see Fig. 3(m)] and give a nematic or
smectic crystal.

In all three skyrmion phases, each Moiré supercell contains a single skyrmion
configuration. Therefore, a suitable characteristic length scale associated
with different skyrmions can be defined by the domain wall contour ${\bf
r}_{\rm d}$ where $n_{z}({\bf r}_{\rm d})=0$. This condition is very much
satisfied where the effective magnetic field due to the bottom layer along the
$z$-direction vanishes, i.e., $B^{z}_{ex}({\bf r}_{\rm d})+B^{z}_{D}({\bf
r}_{\bf d})=0$. From the expression for $B_{\rm ex}$, and $B_{\bf D}$ given in
the Appendix A, we see that ${\bf r}_{\rm d}$ depends on the ratio
$J_{\perp}/J_{\rm D}$, and the inter-layer distance $d$ for a given Moiré
lattice. Depending on the ratio $J_{\perp}/J_{\rm D}$, the condition can turn
into an equation of a circle or an ellipse, as we also find numerically. In
Phase-III, the domain wall takes an elliptic shape in Fig. 3(h) and the two
topological charge centers we find here in Fig. 3(m) sits at the two focal
points.

It is then easy to grasp that a skyrmion phase transition occurs when domain
wall radius ${\bf r}_{\bf d}$ either shrinks to zero or expands to the Moiré
cell boundary ${\bf R}_{1,2}$. The phase transition between Phase-I to Phase-
II occurs when ${\bf r}_{\rm d}={\bf R}_{1}$ or ${\bf R}_{2}$. The rotational
symmetry breaking renders a small domain wall to form in Phase-III with
opposite polarity at one of the Moiré supercell’s site for small values of
$J_{\perp}$, as seen in Fig. 3(c). This small domain wall then shrinks to zero
with the sign reversal of $J_{\perp}$ which disfavors the domain wall of
negative polarity. Curiously, there still exists a finite vortex structure and
a finite charge density in the trivial Phase-IV. But owing to the absence of a
compact domain wall of the polarity density, the net topological winding
number vanishes. Hence, the Phase-IV corresponds to a quasi-uniform
ferromagnertic (or antiferromagnetic phase if $J_{||}<0$) as seen in the
untwisted CrI3 bilayer samples.CrI3_DFT  Finally, large $J_{\perp}$ creates
another compact domain wall at the center of the Moiré supercell in which the
topological charge is $Q=-1$.

Figure 4: Antiferroelectric Phase-II.. (a-b) Zoomed-in view of the topological
charge density in blue to red color map is plotted in the background, with
arrows denoting spins ${\bf s}$ in (a), and emergent electric field ${\bf E}$
vector in (b). The spin texture clarifies the formation of topological charge
at the Moiré lattice side, and a pair of antiskyrmions with opposite charge
centers near the Moiré lattice corners. The electric field lines in (b)
confirm the formation of topological dipole moment (long arrows) between the
antiskymion pairs. (c-d) Real space view of the topological charges in many
Moiré unit cells for Phase-II and Phase-III, respectively. (c) We clearly
observe the Néel analog of the ordering of the topological dipole moments,
giving an antiferroelectric phase. (d) As we move from Phase-II to Phase-III,
the antiskyrmion pairs are annihilated and integer topological charges become
confined by a compact domain wall on different lattice sides. Inset: The
structure factor of the corresponding topological charge density, showing no
charge ordering in Phase-II, as opposed to Bragg peaks in Phase-III.

###  III.3 Topological antiferroelectric phase

Phase-II is very intriguing and novel, and requires separate discussions. The
naive spin-texture of this phase (Fig. 1(c)) is reminiscence of a spin spiral
phase. However, unlike in the other trivial spiral phases, here several new
types of topological charge centers are formed as shown in Fig. 3(l) and 4(a).
Firstly, we observe that the polarity density has a streamline flow diagram
with 1D domain wall, see Fig. 3(b). But it fails to commence a compact
geometry to produce full skyrmion charge centers. However, there exists five
sharp charge centers (three inside the $n_{z}=-1$ (blue) region and two
outside). This structure periodically repeats in a smectic pattern. These
charge centers have different origins from the previous three skyrmion phases
with compact polarity density, and result from splitting of the vorticity by
streamline flow of polarity density.

A zoomed-in view of the spin texture on top of the topological charge density,
as shown in Fig. 4(a), unravels the mechanism of these charge centers. The
charge center at the middle (blue-colored charge density) is a ‘meron’-like
structure, but with a fractional charge of $Q\sim 2/9\pm 0.025$. The other
four charge centers form in pairs with opposite sign of charges at the Moiré
zone corners. The corresponding spin textures reveal that they are
antiskyrmions,Antiskyrmion ; Mn2RhSnParkin  with fractional charges $Q\sim\pm
2/9\pm 0.025$. We note that although there are five distinct meron like charge
centers with $\pm 2/9$, they are shared between neighboring Moiré supercells.
Apart from these concentrated point charges, another +2/9 (red) is spread over
the Moiré cell which is not clearly visible in the colormap of Fig. 4(a),(b).
Hence, the total charge (point-charge and background-charge) within a given
Moire supercell vanishes.

Each pair of oppositely charged antiskyrmions acts as a topological electric
dipole, sitting at each lattice site of the Moiré lattice. The dipole moment
is estimated to be $(0.05\pm 0.005)a$, where $a$ is the Moiré lattice
constant. More interestingly, the dipoles are aligned anti-ferroelectrically
between the nearest neighbors. This gives a Néel-like ordering of the
topological dipole moments, in close analogy to the Néel order of magnetic
moments in a honeycomb lattice.

Comparisons of the spin-textures between Phase-II (anti-ferroelectric) and
Phase-III (skyrmions) throw light on the phase transition between them. The
phase transition occurs when the antiskyrmion pairs coalesce. This also
results in the closing of the streamline flow of the polarity density to form
a compact domain wall (see Fig. 3(c)). Hence a topological winding number
description becomes appropriate in Phase-III.

###  III.4 Electromagnetic duality

To visualize the formation of the dipole moment in Phase-II, we write down a
gauge-dual theory, and calculate the topological electric field lines. In an
electronic quantum Hall analog, we know that a physical charge is attached to
a magnetic flux via the Chern-Simons coupling $-$ also known as the Stréda
formula.Streda  Using this analogy, we affix an emergent gauge field
$a_{\mu}=(\phi,{\bf a})$ with the topological charge density $\rho$ as
$\mathcal{J}_{\mu}=\epsilon_{\mu\nu\tau}\mathcal{F}_{\nu\tau}$. Here
$\mathcal{F}_{\nu\tau}=\partial_{\nu}a_{\tau}$ is the corresponding curvature
field tensor (an emergent electro-magnetic field, but not the same ${\bf B}$
field seen by local spins). The emergent ‘electric field’ is read as
$\mathbf{E}=-{\bf\nabla}\phi-\partial_{t}{\bf a}$. This ‘electric field’
follows the Gauss’ law, and acquires distinct spatial dependence according to
the topological charge distributions. Using the Gauss law, the electric field
lines can be found numerically from $\int{\bf E}.d{\bf S}=\rho$. The electric
field lines shown by black arrows in Fig. 4(b) confirm the existence of the
dipole moment and their antiferroelectric ordering. We extend the calculations
of electric field lines to all the other phases, and find that the the field
line configurations are consistent with the ring-, rod-, and point-like charge
centers as obtained in Phase-I, Phase-II, and Phase-V, respectively (see Fig.
9).

###  III.5 Topological correlation function

In what follows, the topological dipole moment indeed provides a measure of
the order parameter in Phase-II: it continuously disappears in going from
Phase-II to Phase-III. To elucidate this further, we compare larger views of
the Moiré lattice in Phase-II and Phase-III in Figs. 4(c) and 4(d). We also
calculate the topological charge susceptibility as $\chi_{c}({\bf r})=\int
d^{2}{\bf r}^{\prime}\rho({\bf r^{\prime}})\rho({\bf r+r^{\prime}})$. The
corresponding Fourier transformation gives the structure factor $S_{c}({\bf
q})$, which mimics the x-ray diffraction (XRD) pattern as seen by transmission
electron microscopy (TEM), and may be indirectly probed by Lorentz
TEM.Tokura10  The corresponding structure factors are plotted in the inset
figures to Figs. 4(c) and 4(d). As expected, we observe Bragg peaks in Phase-
III as the topological charges form a triangular lattice. The lattice also
features a broken spatial rotation symmetry, and gives a smectic-like skyrmion
lattice. In phase-II the charge centers do not exhibit any Bragg peak up to
the third Brillouin zone. As the appearance of the new Bragg like peaks at the
antiferromagnetic wavevector in the spin-spin correlation function, the
antiferroelectric phase of the topological dipole moment may also be
visualized in the dipole-dipole correlation function.

##  IV Discussions and Conclusions

So far, our discussions were primarily devoted to delineate the mechanism of
skyrmion charges and antiskyrmion dipoles. Although not our primary focus
here, it is worthwhile to outline few possible mechanism to destroy
topological configurations, and the corresponding phase transitions. We
discussed that the size of the compact domain wall ($r_{\rm d}$) as it reduces
to zero or expands all the way to the Moiré lattice vector as a function of
$J_{\perp}/J_{D}$ and the bilayer thickness $d$, it destroys the skyrmion
configuration. Moreover, energetic makes another dominant factor to
destabilizing the saddle-point minima of skyrmions. Once fluctuations are
included, the topological charge centers oscillate, creating ‘phonon’ like
excitations, which melt the skyrmion crystals.MeltSkyrmion  We have studied
the short range (nearest-neighbor) topological charge correlation function,
and find that it exhibits a similar divergence behavior at all phase
transition points (see Appendix C). In addition, we observe that the phase
transition from Phase-II to Phase-III occurs via the coalescence of the
antiskyrmion pairs and vanishing dipole moment. This is reminiscence of the
Kosterlitz-Thouless (KT) like transition, but generalized to the $O(3)$ field.

Can we probe the topological electric field, dipole moment and the KT
transition of the topological charge of the skyrmions? The electric field is
the mediator of the force between two charge particles. Since a skyrmion
charge is a topological charge, it cannot be destroyed without deforming all
the spins in a skyrmion. This prohibits two skyrmions of the same charge to
come close to each other $-$ as if they experience a Coulomb repulsion between
them. This phenomenon can be associated with an electric field. Much like how
we measure an electric field by adding a test charge, here one can think of
adding a test skyrmion in a skyrmion background, and study its dynamics. It
will be found that the test skyrmion of same (opposite) charge will be
repelled (attracted) from the skyrmion background. Similarly, when two
skyrmions/antiskyrmions of opposite charges are spatially separated as in
Phase-II, one can associate a dipole moment in the usual way. With tuning, the
two opposite charges can either annihilate each other or the two charges can
become unbound from each other. The second phenomena is analogous to the KT
transition as seen in the vortex case. In the case of the KT transition in
vortices, pairs of vortices of opposite charges are energetically favorable at
low-temperature, which forms dipole moment and the material behaves as a
dielectric. With increasing temperature, the vortex pairs split and the vortex
charges become unbound, giving a plasma like phase. The phase transition
between them is denoted by the KT transition. In our present case, we can
speculate that a similar splitting of the dipole may occur with increasing
temperature in Phase-II, giving us a unique opportunity to explore a possible
KT transition of skyrmions.

Finally, the spin-spin correlation also plays an important role. We have
calculated the transverse and longitudinal spin-spin correlation functions in
the static limit. The transverse component does not have any $q\sim 0$ mode,
and remains nearly unchanged across all the phases. The longitudinal
susceptibility shows Bragg peaks in four phases, but not in Phase-I. Because,
Phase-I is non magnetic. These additional results and discussions are given in
Appendix.

One may wonder how sensitive is the the phase diagram to DMI and SOC terms. We
have checked that DMI brings in very little change to the spin configurations,
and their topological properties are robust as long as the DMI strength is
considerably weaker than $J_{\perp}$, $J_{\rm D}$. With an eye to synthesize
twisted bilayers of VdW magnets,CrI3 ; Chalcogenides  is it known that the
spins are local in nature, and the materials are charge insulators. Thus the
SOC does not play an important role to the skyrmion configurations.

##  Appendix A Effective inter-layer magnetic field

Figure 5: Zoomed views of the spin structure in five different phases.

As mentioned in the main text, the effect of $H_{\rm ex}$ and $H_{\rm D}$ can
be described by effective magnetic field ${\bf B}({\bf r}_{i})$ coming from
exchange (${\bf B}_{\rm ex}({\bf r}_{i})$) and dipole (${\bf B}_{\rm D}({\bf
r}_{i})$). Taking all the bottom layer’s spin S to be along the
$+z$-direction, we can explicitly write the field components as

$\displaystyle{\bf B}_{\rm ex}({\bf r}_{i})$ $\displaystyle=$ $\displaystyle
J_{\perp}\sum_{a=1,2,3}\cos({\bf G}_{a}.{\bf r}_{i})\hat{{\bf z}},$ (6)
$\displaystyle{\bf B}_{\rm D}({\bf r}_{i})$ $\displaystyle=$
$\displaystyle\sum_{j}\frac{J_{D}}{r_{ij}^{3}}\left[1-\frac{3d^{2}}{r_{ij}^{2}}\right]\hat{{\bf
z}}$ (7) $\displaystyle-\frac{3J_{D}d}{r_{ij}^{5}}\left[(x_{i}-x_{j})\hat{{\bf
x}}+(y_{i}-y_{j})\hat{{\bf y}}\right].$

The index $i$ denotes a top layer spin. The sum over $j$ denotes a sum over
bottom layer spins and is restricted up to $r_{\rm cut}=N_{\rm cut}a_{0}$. Due
to the higher power of $r_{ij}$ in the denominator, it is easy to check that
the summation converges very rapidly ($\sim N_{\rm cut}<20$), much before the
Moire supercell lattice vector $R_{1,2}\sim 35a_{0}$. ${\bf G}_{a}$ denotes
three lowest Moiré reciprocal lattice vectors, $d$ is the inter-layer
distance, and $x_{i}$, $y_{i}$ are $x$ and $y$ coordinates of $i^{\rm th}$
spin. ${\bf B}_{\rm ex}$ is plotted in fig. 6(a) and (c) for $J_{\perp}<0$
(ferromagnetic) and $J_{\perp}>0$ (antiferromagnetic) respectively.

Clearly, a compact domain wall forms at the nodal contour of $s_{z}({\bf r})$
where the total magnetic field roughly vanishes, i.e., $B_{z}({\bf r}_{\rm
d})=0$. (This approximation works better where the in-plane spin exchange
$J_{||}<<|B|$, so that the second term in Eq. (4) dominates). Then the
condition simplifies to

$\frac{J_{\perp}}{J_{\rm D}}\beta_{ex}({\bf r}_{d})+\beta_{D}({\bf r}_{d})=0.$
(8)

where ${\bf r}_{d}$ is the locii of the domain wall boundary and $\beta_{ex}$
and $\beta_{D}$ are given by

$\displaystyle\beta_{ex}({\bf r}_{d})$ $\displaystyle=$
$\displaystyle\sum_{a}\cos({\bf G}_{a}.{\bf r}_{d})$ (9)
$\displaystyle\beta_{D}({\bf r}_{d})$ $\displaystyle=$
$\displaystyle\sum_{j}\frac{1}{|{\bf r}_{d}-{\bf
r}_{j}|^{3}}\left[1-\frac{3d^{2}}{|{\bf r}_{d}-{\bf r}_{j}|^{2}}\right].$ (10)

It is not easy to find an analytical expression for this nodal contour from
Eqs. (6), and (7). But it’s clear that the value of ${\bf r}_{d}$ depends on
the $J_{\perp}/J_{\rm D}$ ratio and the bilayer thickness $d$ for a given
Moiré lattice denoted by ${\bf G}$. The equation of the nodal line can be a
circle or an ellipse depending on these parameters. As we see form the
numerical simulation, the domain wall is circular for large values of
$J_{\perp}$, while it takes an elliptical form in Phase-III for small values
of $J_{\perp}$.

Again in the limit of $J_{||}<<|B|$, a phase transition is defined at
$J_{\perp}/J_{\rm D}=-\beta_{D}({\bf r}_{d})/\beta_{ex}({\bf r}_{d})$ (from
Eq. (8)) for different values of ${\bf r}_{d}$. In phase I, the radius of
circular patch increases as $J_{\perp}\rightarrow 0$, and at the transition
point to phase II two circular regions merge together at ${\bf r}_{d}={\bf
R}_{1}/2$. On the other hand, ${\bf r}_{d}$ decreases in phase V as
$J_{\perp}\rightarrow 0$, which implies ${\bf r}_{d}=0$ at the transition from
the phase V to phase IV. Similarly transitions from phase II to III and III to
IV are given by ${\bf r}_{d}={\bf R}_{2}/2$ and ${\bf r}_{d}=({\bf R}_{1}+{\bf
R}_{2})/2$ respectively. By numerically evaluating $\beta_{ex}$ and
$\beta_{D}$ at various ${\bf r_{d}}$, we extract the approximate critical
value of the ratio $J_{\perp}/J_{D}$ for different transitions which are
listed in the table below:

Transitions ${\bf r}_{d}$ $J_{\perp}/J_{D}$ Phase I to phase II ${\bf
R}_{1}/2$ -0.57 Phase II to phase III ${\bf R}_{2}/2$ -0.4 Phase III to phase
IV $({\bf R}_{1}+{\bf R}_{2})/2$ 0.3 Phase IV to phase V ${\bf 0}$ 1.6

A direct comparison of transition from phase ‘I to II’ and phase ‘IV to V’
reveals that $|J_{\perp}/J_{D}|$ is larger on the positive side ($1.6$ on
positive side and $-0.57$ on negative side). This can be understood from Fig.
6. At the corner of the Moiré supercell $B_{\rm ex}$ and $B_{\rm D}$ have the
same sign in the $J_{\perp}<0$ region. So their effects add up. This is also
responsible for a thin domain wall of magnetic bubbles in phase I. In the
positive $J_{\perp}$ region, however, $B_{\rm ex}$ and $B_{\rm D}$ have
opposite sign. So a larger value of $J_{\perp}$ is needed to overcome the
effect $J_{\bf D}$ in positive side. This leads to a wider domain wall in
phase V. We should distinguish the effects of the dipole-dipole interaction
$H_{\rm D}$ of the present study compared to other studies. In earlier
studies,SkyrmionDDI  such an interaction is involved for the same intra-layer
spin (between ${\bf s}$ and ${\bf s}$ variables) which tends to produce
magnetic bubble phases.SkyrmionReview1  On the other hand, here the dipole
interaction is between the inter-layer (${\bf s}$ and ${\bf S}$), and is
detrimental to the bubble or Bloch skyrmion phases, while promote the
streamline flow of the polarity density.

Figure 6: Effective magnetic fields along the $z$-direction in the Moiré unit
cell. (a) Field due to the inter-layer exchange interaction ($B_{ex}$) for
$J_{\perp}<0$ and (c) for $J_{\perp}>0$. (b) Field due to the inter-layer
dipole interaction ($B_{D}$) for $J_{D}>0$. ${\bf{R}}_{1}$ and ${\bf{R}}_{2}$
denote the two translation lattice vectors.

##  Appendix B Detail of Monte Carlo method

As mentioned in the main text, the Hamiltonian for the top layer can be
written as

$\displaystyle H=H_{1}+\sum_{i}{\bf{B}}({\bf{r}}_{i}).{\bf{s}}_{i}.$ (11)

Here $H_{1}$ is the intra-layer Hamiltonian and ${\bf{B}}({\bf{r}}_{i})$ is
the effective magnetic field due to inter-layer interaction. The lattice site
${\bf{r_{i}}}$ spans the Moiré unit cell. In our calculation at the
commensurate angle $1.61^{0}$ we have 2522 basis site per unit cell. We set
all the classical spins to have unit length ($S_{i}=1$) so that all the spins
can be specified with two parameters $S_{z}$ and $\phi$ where $\phi$ is the
angle that the component of spin on the $xy$-plane ($s^{xy}$) makes with the
$x$ axis. From these two parameters, we can extract all the three components
of the spins as

$\displaystyle s^{xy}_{i}=\sqrt{1-(s^{z}_{i}})^{2},$ $\displaystyle
s^{x}=s^{xy}\cos(\phi),$ $\displaystyle s^{y}=s^{xy}\sin(\phi).$ (12)

We initialize the simulation with all spins pointing upward i.e.,
$S^{z}_{i}=1$, and then the next spin is chosen randomly from all the lattice
points. The update algorithm for that spin is given by

$\displaystyle s^{z}_{i}=s^{z}_{i}+\gamma ds^{z},$
$\displaystyle\phi_{i}=\phi_{i}+\gamma d\phi.$ (13)

And if $|s^{z}_{i}|>1$ then

$\displaystyle s^{z}_{i}=2\mp(s^{z}_{i}+\gamma ds^{z}),$
$\displaystyle\phi_{i}=\phi_{i}+\gamma d\phi+\pi.$ (14)

Here the $\pm$ signs are for $s^{z}_{i}>1$ and $s^{z}_{i}<-1$, respectively.
$\gamma$ is a random number between 1 and -1, and $ds^{z}$ and $d\phi$ are
ranges of $s^{z}$ and $\phi$.

At each Monte Carlo step we calculate the local Hamiltonian

$\displaystyle H({\bf r}_{i})_{loc}$ $\displaystyle=$ $\displaystyle
H_{1}({\bf r}_{i})-\sum_{i}{\bf{B}}({\bf{r}}_{i}).{\bf{s}}_{i},$ (15)
$\displaystyle H_{1}({\bf r}_{i})$ $\displaystyle=$
$\displaystyle\sum_{<i,j>}J_{||}{\bf{s}}_{i}.{\bf{s}}_{j}.$ (16)

and each configuration is accepted with a Boltzmann probability $e^{[H({\bf
r}_{i})_{\rm loc}^{\rm new}-H({\bf r}_{i})_{\rm loc}^{\rm old}]/k_{B}T}$ at
temperature $T$. To find the low temperature Monte Carlo ground state we
perform temperature annealing as well as parameter annealing, and check for
convergence of the total ground state energy. We use $dS^{z}=0.4$ and
$d\phi=0.4\pi$. At low temperature we decrease the value of $ds^{z}(=0.1)$ and
$d\phi(=0.1\pi)$ to increase the acceptance ratio in Monte Carlo steps. We
equilibriate the system with $10^{8}$ steps and take the low energy
configuration for another $10^{8}$ steps.

##  Appendix C Correlation functions and structure factors

Figure 7: Structure factors. (Left) Left row gives the structure factor of the
topological charge density (Fourier transformation of the charge-charge
correlation function). (Middle) Transverse spin-spin correlation functions.
(Right) Longitudinal correlation function in the momentum space. In all three
cases, we give various quadrants of the Brillouin zone with dashed lines.

From Fig. 2, it is evident that there are finite topological charge
distributions even in the trivial topological phases where the polarity
density $n_{z}$ fails to create domain walls. Therefore, to obtain a
microscopic nature of how the quantum fluctuations make it possible to form a
topological configurations in Phase-I, III and V, we investigate the following
correlation functions: Topological charge susceptibility

$\displaystyle\chi_{c}({\bf r})=\int d^{2}{\bf r}^{\prime}\rho({\bf
r^{\prime}})\rho({\bf r+r^{\prime}}).$ (17)

The spin-non-flip (or skyrmion polarity) and spin-flip correlation functions

$\displaystyle\chi_{s}^{zz}({\bf r})$ $\displaystyle=$ $\displaystyle\int
d^{2}{\bf r}^{\prime}n_{z}({\bf r^{\prime}})n_{z}({\bf r+r^{\prime}}),$
$\displaystyle\chi_{s}^{\pm}({\bf r})$ $\displaystyle=$ $\displaystyle\int
d^{2}{\bf r}^{\prime}n_{+}({\bf r^{\prime}})n_{-}({\bf r+r^{\prime}}),$ (18)

where $n_{\pm}=n_{x}\pm in_{y}$. The Fourier transformation of the correlation
function gives the structure factor $S_{c/s}({\bf q})$. For the topological
charge density, the topological structure factor mimics the XRD or TEM maps,
except that this is not directly measurable and indirectly it can be mapped by
Larentz TEM. The two spin structure factors are however measurable via small-
angle neutron scattering (SANS) experiments.

First we notice that in the non-trivial phases, the topological structure
factor shows distinct Bragg peaks at ${\bf q}\rightarrow 0$, signifying the
transnational symmetry breaking and formation of skyrmion lattice. (When the
dynamics is added, this peak will disperse away from ${\bf q}=0$ point as an
typical acoustic Goldstone mode). In the trivial phase (Phase-II and Phase-
IV), there is no Bragg peaks, but some weak intensity at higher Brillouin zone
of the original Moiré supercell) which are often observed in dilute gas or
liquid phases.

The spin-non-flip and spin slip components of the spin-structure factors have
the usual behavior in a skryrmions structure, while the stark difference
between the Phase-I and Phase-IV for the $\chi_{s}^{zz}$ should be noted. In
the case of a ring charge (Phase-I) there is no ${\bf q}\rightarrow 0$ mode
for the polarity density since here we have a sharp domain wall boundary
between the up and down spin states. For all other cases, the ${\bf
q}\rightarrow 0$ mode exists, suggesting a finite value of the magnetic moment
in these cases on the top layer. On the other, in Phase-II, the total
magnetization vanishes in a Moiré unit cell, however, if one define Néel
magnetization between the up spin domain and down-spin domain, we have an
antiferromagnetic like domain wall ordering. In all cases, the spin-flip
structure factor is similar, while an additional spatial rotational symmetry
breaking is observed in Phases II and III as expected. We find that that all
the peaks are present in the second Brillouin zone which is due to the fact
that the Moire cell here forms an honeycomb lattice which has two sublattices.

Figure 8: plot of (a) Average $\chi_{c}$ and (b) average magnetization
($m_{z}$) and derivative of $m_{z}$ as a function of interlayer exchange
($J_{\perp}/J_{||}$) at a fixed $J_{D}=0.3J_{||}$. Discontinuity in $m_{z}$ at
the phase transition can be seen from the jumps in derivative of $m_{z}$ (blue
curve in (b))

To further elucidate the phase transition, we analyze the topological
correlation function at a fixed distance ($r_{i}-r_{j}$=constant):

$\displaystyle\chi_{c}=\frac{1}{N}\sum_{<i,j>}\rho(\bf{r}_{i})\rho(\bf{r}_{j}).$
(19)

We plot the result as a function of inter-layer exchange ($J_{\perp}$) in Fig.
8(a). This gives the short-range correlation of the topological charge
density.

We find the correlation has peaks at the phase transition points as shown in
Fig. 8(a). This result is consistent with jumps in the effective magnetic
field $M_{z}=\left\langle(1/N)\sum_{i}S_{i}^{z}\right\rangle$ along
$J_{\perp}$, as shown in Fig. 8(b). Near the phase transition points, the
landscape of effective magnetic field changes over the Moiré unit cell. As a
result near the domain wall in each phase (except Phase IV where there is no
domain wall) ${\bf B}$ field becomes small over a large area. Therefore, the
Hamiltonian is mostly dominated by strong in-plane ferromagnetic exchange
($J_{||}$), and produces a strong nearest neighbour correlation.

We further calculate scalar chirality (not
shown)$\chi_{sc}=\frac{1}{N}\left\langle\sum_{<i,j,k>}{\bf S}_{i}.\left({\bf
S}_{j}\times{\bf S}_{k}\right)\right\rangle$, where $\left\langle
i,j,k\right\rangle$ denotes three spins forming the smallest triangle in the
honeycomb lattice. We find that the scalar chirality vanishes in Phase II and
Phase IV, where the magnetization has a sharp jump at their phase boundary
(Fig. 8(b)).

##  Appendix D Dipole and quadrupole moments in Phase II

Figure 9: (a-e)We have plotted the emergent topological electric field lines
on top of the topological charge density in the vicinity of the charge centers
(see Fig. 10 for larger view) in all five phases. In (f) we show the schematic
distribution of point charges for Pase II with origin of the coordinates at
the cahrge $q_{3}$.

Field lines for topological charge in Phase II is shown in Fig. 9(c). The
field lines indicate that the charge distribution can be approximated with
point charges ($q_{i}$, $i=1,2,3,4,5$) as shown in Fig. 9(f). As the total
charge of this configuration is zero we investigate the higher moments (dipole
and quadrupole) of this charge distribution. We calculate dipole moment for
charge pair $(q_{1},q_{2})$ and $(q_{4},q_{5})$ and taking the origin of our
co-ordinate at one of the charge center (see Fig. 9(f)) we calculate the
quadrupole moments which are given by

$\displaystyle Q$ $\displaystyle=$ $\displaystyle\left[\begin{array}[]{ c c
}Q_{xx}&Q_{xy}\\\ Q_{yx}&Q_{yy}\end{array}\right]\ \ =\left[\begin{array}[]{ c
c }\sum_{i}\rho_{i}x_{i}x_{i}&\sum_{i}\rho_{i}x_{i}y_{i}\\\
\sum_{i}\rho_{i}y_{i}x_{i}&\sum_{i}\rho_{i}y_{i}y_{i}\end{array}\right]$ (24)

Values of the dipole and quadrupole moments for Phase II are given in the
following table.

quantities |  $J_{D}=0.53$ |  $J_{D}=1.0$
---|---|---
$q_{1}$ |  -0.23 |  -0.23
$q_{2}$ |  0.22 |  0.22
$q_{3}$ |  -0.25 |  -0.24
$q_{4}$ |  0.22 |  0.24
$q_{5}$ |  -0.20 |  -0.22
$(x_{1},y_{1})$ in umit of moiré lattice vector R=35.5 |  (0.20 , 0.07) |  (0.20 , 0.07)
$(x_{2},-y_{2})$ in umit of moiré lattice vector R=35.5 |  (0.35 , -0.09) |  (0.35 , -0.09)
${\bf p}_{1}$ |  0.033$\hat{\bf x}$ \+ 0.035$\hat{\bf y}$ |  0.033$\hat{\bf x}$ \+ 0.035$\hat{\bf y}$
${\bf p}_{2}$ |  -0.033$\hat{\bf x}$ \- 0.035$\hat{\bf y}$ |  -0.034$\hat{\bf x}$ \- 0.037$\hat{\bf y}$
$Q_{xx}$ |  -0.035 |  -0.037
$Q_{yy}$ |  -0.001 |  -0.001
$Q_{xy}$ |  0.021 |  0.020

We caution the readers that the estimate of the fractional charge is subjected
to the integration contour around each charge centers, since there is no
compact of domain wall boundary here.

Figure 10: Plots of topological charge density in many Moir’e unit cells.

## References

  * (1) T. H. R. Skyrme, A non-linear field theory. Proc. Royal Soc. Lond. A: Math. Phys. Eng. Sci. 260, 127 (1961); T. H. R. Skyrme, A unified field theory of mesons and baryons. Nucl. Phys. 31, 556 (1962).
  * (2) A. N. Bogdanov and A. Hubert, Sov.Phys.JETP 68, 101-103 (1989); J. Magn. Magn. Mater. 138 , 255-269 (1994); U. K. Roessler, A. N. Bogdanov, and C. Pfleiderer, Nature 442 , 797-801 (2006).
  * (3) N. Nagaosa and Y. Tokura, Nature Nanotechnology 8, 899-911 (2013).
  * (4) J. H. Han, Skyrmions in Condensed Matter, Springer Tracts in Modern Physics, Volume 278
  * (5) Y. S. Lin, J. Grundy, and E. A. Giess, Appl. Phys. Lett. 23, 485–487 (1973); A. P. Malozemoff, and J. C. Slonczewski, Magnetic Domain Walls in Bubble Materials 306–314 (Academic, 1979).; T. Garel, and S. Doniach, Phys. Rev. B 26, 325–329 (1982); T. Suzuki, J. Magn. Magn. Mater. 31–34, 1009–1010 (1983); B. Göbel, J. Henk, I. Mertig, Sci. Reports 9, 9521 (2019). M. Ezawa, Phys. Rev. Lett. 105, 197202 (2010).
  * (6) S. Rohart, and A. Thiaville, Phys. Rev. B 88, 184422 (2013).
  * (7) M. Hoffmann, B. Zimmermann, G. P. Muller, D. Schurhoff, N. S. Kiselev, C. Melcher, and S. Blugel. Nat. comm. 8, 308 (2017).
  * (8) E. Y. Vedmedenko, J. A. Arregi, P. Riego, A. Berger. Phys. Rev. Lett. 122, 257 (2019)
  * (9) S. Rohart and A. Thiaville, Phys. Rev. B Condens. Matter Mater. Phys. 88, 184422 (2013).
  * (10) B. Dupé, G. Bihlmayer, M. Böttcher, S. Blügel, and S. Heinze, Nat. Commun. 7, 11779 (2016).
  * (11) T. Okubo, S. Chung and H. Kawamura, Phys. Rev. Lett. 108, 017206 (2012).
  * (12) A. O. Leonov, and M. Mostovoy, Nat. Commun. 6, 8275 (2015).
  * (13) C. D. Batista, S. Z. Lin, S. Hayami, and Y. Kamiya, Rep. Prog. Phys. 79, 084504 (2016).
  * (14) T. Okubo, S. Chung, and H. Kawamura, Multiple-q states and the skyrmion lattice of the triangular-lattice Heisenberg antiferromagnet under magnetic fields. Phys. Rev. Lett. 108, 017206 (2012); A. Leonov, and M. and Mostovoy, Multiply periodic states and isolated skyrmions in an anisotropic frustrated magnet. Nat. Commun. 6, 8275 (2015); S. Hayami, S. -Z. Lin, and Cristian D. Batista Phys. Rev. B 93, 184413 (2016); C. D. Batista, S. -Z. Lin, S. Hayami, and K. Yoshitomo, Reports on Progress in Physics, 79, 084504 (2016).
  * (15) S. Banerjee, O. Erten, and M. randeria, Nature Physics 9, 626-630 (2013).
  * (16) R. Ozawa, S. Hayami, and Y. Motome, Phys. Rev. Lett. 118, 147205 (2017).
  * (17) I. Gross, W. Akhtar, A. Hrabec, J. Sampaio, L. J. Martinez, S. Chouaieb, B. J. Shields, P. Maletinsky, A. Thiaville, S. Rohart, and V. Jacques, Phys. Rev. Mater. 2, 024406 (2017).
  * (18) S. L. Sondhi, A. Karlhede, S. A. Kivelson, and E. H. Rezayi, Phys. Rev. B 47, 16419 (1993).
  * (19) L. S. Leslie, A. Hansen, K. C. Wright, B. M. Deutsch, and N. P. Bigelow, Phys. Rev. Lett. 103, 250401 (2009).
  * (20) S. Mühlbauer, B. Binz, F. Jonietz, C. Pfleiderer, A. Rosch, A. Neubauer, R. Georgii, and P. Böni, Science 323, 915-919 (2009).
  * (21) W. Münzer, A. Neubauer, T. Adams, S. Mühlbauer, C. Franz, F. Jonietz, R. Georgii, P. Böni, B. Pedersen, M. Schmidt, A. Rosch, and C. Pfleiderer Phys Rev. B 81, 041203(R) (2010).
  * (22) X. Z. Yu, Y. Onose, N. Kanazawa, J. H. Park, J. H. Han, Y. Matsui, N. Nagaosa and Y. Tokura, Nature 465, 901–904 (2010); D. Morikawa, K. Shibata, N. Kanazawa, X. Z. Yu, and Y. Tokura, Phys. Rev. B 88 024408 (2013).
  * (23) P. Milde, D. Köhler, J. Seidel, L. M. Eng, A. Bauer, A. Chacon, J. Kindervater, S. Mühlbauer, C. Pfleiderer, S. Buhrandt, C. Schütte, and A. Rosch, Science 340, 1076-1080 (2013).
  * (24) S. Seki, X. Z. Yu, S. Ishiwata, and Y. Tokura, Science 336, 198-201 (2012).
  * (25) Y. Tokunaga, X. Z. Yu, J. S. White, H. M. Ronnow, D. Morikawa, Y. Taguchi, and Y. Tokura, Nat. Commun. 6, 7638 (2015).
  * (26) X.Zhao, C. Jin, C. Wang, H. Du, J. Zang, M. Tian, R. Che, and Y. Zhang, Proc. Natl. Acad. Sci. 113, 4918 (2016).
  * (27) S. Heinze, K. von Bergmann, M. Menzel, J. Brede, A. Kubetzka, R. Wiesendanger, G. Bihlmayer, S. Blügel, Nat. Phys. 7, 713-718 (2011);
  * (28) Z. Hou, W. Ren, B. Ding, G. Xu, Y. Wang, B. Yang, Q. Zhang, Y. Zhang, E. Liu, F. Xu, W. Wang, G. Wu, X. -X. Zhang, B. Shen, Z. Zhang, Advanced Materials 29, 1701144 (2017).
  * (29) R. Wiesendanger, Nanoscale magnetic skyrmions in metallic films and multilayers: a new twist for spintronics. Nat. Rev. Mat. 1, 16044 (2016).
  * (30) A. K. Nayak, V. Kumar, T. Ma, P. Werner, E. Pippel, R. Sahoo, F. Damay, U. K. Rü$\beta$ler, C. Felser and S. S. P. Parkin, Nature 548, 561-566 (2017).
  * (31) W. Koshibae, and N. Nagaosa, Nat. Commun., 7, 10542 (2016); M. Hoffmann, B. Zimmermann, G. P. Müller, D. Schürhoff, N. S. Kiselev, C. Melcher, and S. Blügel Nat. Commun. 8, 308 (2017)
  * (32) N. Romming, C. Hanneken, M. Menzel, J. E. Bickel, B. Wolter, K. von Bergmann, A. Kubetzka, and R. Wiesendanger, Science 341, 636-639 (2013).
  * (33) J. Sampaio, V. Cros, S. Rohart, A. Thiaville, and A. Fert, Nature Nanotechnology 8, 839–844 (2013).
  * (34) M. Mochizuki, X. Z. Yu, S. Seki, N. Kanazawa, W. Koshibae, J. Zang, M. Mostovoy, Y. Tokura, and N. Nagaosa, Nature Materials 13, 241–246 (2014).
  * (35) W. Jiang, P. Upadhyaya, W. Zhang, G. Yu, M. B. Jungfleisch, F. Y. Fradin, J. E. Pearson, Y. Tserkovnyak, K. L. Wang, O. Heinonen, S. G. E. te Velthuis, and A. Hoffmann, Science 349, 283-286 (2015).
  * (36) A. Fert, N. Reyren, and V. Cros, Nature Reviews Materials, 2, 17031 (2017); A. Fert, V. Cros, and J. Sampaio, Skyrmions on the track. Nat. Nanotechnol. 8, 152-156 (2013).
  * (37) S. S. P. Parkin, M. Hayashi, and L. Thomas, Science 320, 190–194 (2008).
  * (38) R. Tomasello, E. Martinez, R. Zivieri, L. Torres, M. Carpentieri, and G. Finocchio, Sci. Rep. 4, 6784 (2014); P. F. Bessarab, G. P. Müller, I. S. Lobanov, F. N. Rybakov, N. S. Kiselev, H. Jónsson, V. M. Uzdin, S. Blügel, L. Bergqvist, and A. Delin, Sci. Reports 8, 3433 (2018).
  * (39) X. Wang, K. Du, Y. Y. F. Liu, P. Hu, J. Zhang, Q. Zhang, M. H. S. Owen, X. Lu, C. K. Gan, and P. Sengupta, 2D Materials 3, 3, 031009 (2016).
  * (40) J. U. Lee, S. Lee, J. H. Ryoo, S. Kang, T. Y. Kim, P. Kim, C. H. Park, J. G. Park, and H. Cheong, Nano Lett. 16, 12, 7433-7438 (2016).
  * (41) P. A. Joy, and S. Vasudevan, Phys. Rev. B 46, 5425-5433 (1992).
  * (42) C. C. Mayorga-Martinez, Z. Sofer, D. Sedmidubsky, S. Huber, A. Y. S. Eng, and M. Pumera, ACS Appl. Mater. Interfaces 9, 12563–12573 (2017).
  * (43) Z. Fei, B. Huang, P. Malinowski, W. Wang, T. Song, J. Sanchez, W. Yao, D. Xiao, X. Zhu, A. F. May, W. Wu, D. H. Cobden, J. -H. Chu, and X. Xu, Nat.Materials, 17, 778–782 (2018).
  * (44) H. J. Deiseroth, K. Aleksandrov, C. Reiner, L. Kienle, and R. K. Kreme, Eur. J. Inorg. Chem., 1561–1567 (2006).
  * (45) M. Kan, S. Adhikari, and Q. Sun, Phys. Chem. Chem. Phys. 16, 4990-4994 (2016).
  * (46) C. Ataca, andH. Sahin, and S. Stable, J. Phys. Chem. C 116, 8983-8999 (2012).
  * (47) Y. Ma, Y. Dai, M. Guo, C. Niu, Y. Zhu, B. Huang, ACS Nano 6, 1695-1701 (2013).
  * (48) M. Bonilla, S. Kolekar, Y. Ma, H. C. Diaz, V. Kalappattil, R. Das, T. Eggers, H. R. Gutierrez, M. -H. Phan, and M. Batzill, Nat. Nanotechnology 13, 289-293 (2018).
  * (49) G. Duvjir, B. K. Choi, I. Jang, S. Ulstrup, S. Kang, T. T. Ly, S. Kim, Y. H. Choi, C. Jozwiak, A. Bostwick, E. Rotenberg, J. -G. Park, R. Sankar, K. -S. Kim, J. Kim, and Y. Jun Chang, Nano. Lett. 18, 5432-5438 (2018).
  * (50) C. Gong, L. Li, Z. Li, H. Ji, A. Stern, Y. Xia, T. Cao, W. Bao, C. Wang, Y. Wang, Z. Q. Qiu, R. J. Cava, S. G. Louie, J. Xia, X. Zhang, Nature 546, 265-269 (2017).
  * (51) B. Huang, G. Clark, E. N.-Moratalla, D. R. Klein, R. Cheng, K. L. Seyler, D. Zhong, E. Schmidgall, M. A. McGuire, D. H. Cobden, W. Yao, D. Xiao, P. J.-Herrero, and X. Xu, Nature 546, 270-273 (2017).
  * (52) B. Huang, G. Clark, D. R. Klein, D. MacNeill, E. N. -Moratalla, K. L. Seyler, N. Wilson, M. A. McGuire, D. H. Cobden, D. Xiao, W. Yao, P. J. -Herrero, and X. Xu, Nat. Nanotechnology 13, 544–548 (2018).
  * (53) S. Jiang, L. Li, Z. Wang, K. F. Mak, and J. Shan, Nat. Nanotechnology 13, 549–553 (2018).
  * (54) N. Sivadas, S. Okamoto, X. Xu, C. J. Fennie, D. Xiao, Nano Lett. 18, 7658-7664 (2018).
  * (55) K. Hejazi, Z. -X. Luo, and L. Balents, PNAS 117, 10721-10726 (2020).
  * (56) T. Bömerich, L. Heinen, and A. Rosch, Phys. Rev. B 102, 100408(R) (2020).
  * (57) Q. Tong, F. Liu, J. Xiao, and W. Yao Nano Lett. 18, 7194-7199 (2018).
  * (58) M. Akram, and O. Erten, Phys. Rev. B 103, L14040 (2021).
  * (59) M. Yang, Q. Li, R. V. Chopdekar, R. Dhall, J. Turner, J. D. Carlström, C. Ophus, C. Klewe, P. Shafer, A. T. N’Diaye, J. W. Choi, G. Chen, Y. Z. Yu, C. Hwang, F. Wang, and Z, Q. Qiu, Science Advances bf 6, eabb5157 (2020).
  * (60) T.-E. Park, L. Peng, J. Liang, A. Hallal, F. S. Yasin, X. Zhang, S. J. Kim, K. M. Song, K. Kim, M. Weigand, G. Schuetz, S. Finizio, J. Raabe, K. Garcia, J. Xia, Y. Zhou, M. Ezawa, X. Liu, J. Chang, H. C. Koo, Y. D. Kim, M. Chshiev, A. Fert, H. Yang, X. Yu, and S. Woo, Phys. Rev. B 103, 104410 (2021).
  * (61) T. Bömerich, L. Heinen, and A. Rosch, arXiv:2004.13684v2;
  * (62) M. A. McGuire, H. Dixit, V. R. Cooper, and B. C. Sales, Chem. Mater. 27, 612–620 (2015); B. Huang, G. Clark, E. N. -Moratalla, D. R. Klein, R. Cheng, K. L. Seyler, D. Zhong, E. Schmidgall, M. A. McGuire, D. H. Cobden, W. Yao, D. Xiao, P. J. -Herrero, and X. Xu, Nature 546, 270–273 (2017).
  * (63) C. Gong, L. Li, Z. Li, H. Ji, A. Stern, Y. Xia, T. Cao, W. Bao, C. Wang, Y. Wang, Z. Q. Qiu, R. J. Cava, S. G. Louie, J. Xia, and X. Zhang, Nature 546, 265–269 (2017); M. Bonilla, S. Kolekar, Y. Ma, H. C. Diaz, V. Kalappattil, R. Das, T. Eggers, H. R. Gutierrez, M. -H. Phan, and M. Batzill, Nat. Nanotech. 13, 289–293 (2018); D. J. ÓHara, T. Zhu, A. H. Trout, A. S. Ahmed, Y. K. Luo, C. H. Lee, M. R. Brenner, S. Rajan, J. A. Gupta, D. W. McComb, and R. K. Kawakami, Nano Lett. 18, 3125−3131 (2018).
  * (64) A. Bansil, H. Lin, T. Das, Rev. Mod. Phys. 88, 021004 (2016). M. Z. Hasan, C. L. Kane, Rev. Mod. Phys. 82, 3045 (2010). X. L. Qi, S. C. Zhang, Rev. Mod. Phys. 83, 83, 1057 (2011).
  * (65) A. Malozemoff, and J. Slonczewski, Magnetic Domain Walls in Bubble Materials (Academic press, 1979); A. H. Eschenfelder, Magnetic bubble technology (Springer-Verlag Berlin Heidelberg New York, 1980).
  * (66) T. Bzdušek, Q. Wu, A. Rüegg, M. Sigrist, and A. A. Soluyanov, Nature 538, 75–78 (2016).
  * (67) P. Streda, J. Phys. C: Solid State Phys., 15, L717-L721 (1982).
  * (68) P. Huang, et al. Nat. Nanotech. 15, 761-767 (2020).