# Controllable magnetic domains in twisted trilayer magnets

Kyoung-Min Kim<EMAIL_ADDRESS>Center for Theoretical Physics of Complex
Systems, Institute for Basic Science, Daejeon 34126, Korea  Moon Jip Park
<EMAIL_ADDRESS>Department of Physics, Hanyang Univercity, Seoul 04763,
Republic of Korea

###### Abstract

The use of moiré patterns to manipulate two-dimensional materials has
facilitated new possibilities for controlling material properties. The moiré
patterns in the two-dimensional magnets can cause peculiar spin texture, as
shown by previous studies focused on twisted bilayer systems. In our study, we
develop a theoretical model to investigate the magnetic structure of twisted
trilayer magnets. Unlike the twisted bilayer, the twisted trilayer magnet has
four different local stacking structures distinguished by the interlayer
couplings between the three layers. Our results show that the complex
interlayer coupling effects in the moiré superlattice can lead to the
stabilization of rich magnetic domain structures; these structures can be
significantly manipulated by adjusting the twist angle. Additionally, external
magnetic fields can easily manipulate these domain structures, indicating
potential applications in spintronics devices.

moiré engineering, twisted trilayer, var der Waals magnets

###### pacs:

75.70.Ak, 75.75.+a, 75.10.Hk

##  I Introduction

The manipulation of two-dimensional materials using moiré patterns has
facilitated exciting possibilities for designing unique material properties
[1, 2, 3, 4, 5, 6, 7, 8]. Recent experimental demonstrations of moiré magnets
in twisted transition metal trihalides MX3 (M = Cr, Ru; X = Cl, Br, I) have
shown the potential to manipulate nanoscale magnetic domain structures by
adjusting the twist angle [9, 10, 11]. Additionally, theoretical studies
indicated the emergence of magnetic skyrmions [12, 13, 14, 15, 16, 17] and
topological magnons [18, 19] within the domain structure, which are promising
for the development of spintronics [20, 21, 22] and multiferroic materials
with electrical tunability [23].

The moiré engineering of the van der Waals interface possesses great potential
for extension to various heterostructures and multilayer systems. Examples of
recent progress include nonmagnetic 2D materials, such as graphene [1, 2],
transition metal dichalcogenides [3], and transition metal oxides [4], as well
as magnetic 2D materials, such as transition metal trihalides [9, 10, 11] with
multilayer moiré superlattices.

Figure 1: Schematic of twist-engineering in trilayer magnets. a Twisted
trilayer magnet is formed by rotating the middle layer of an aligned trilayer
magnet. b Twist generation of the four distinct local interlayer couplings
$J_{ij}^{\textrm{tm,mb}}$, each with either ferromagnetic (FM) or
antiferromagnetic (AFM) characteristics. c Interlayer coupling stabilization
of the four different magnetic domain structures. Red (blue) indicates
magnetization in the positive (negative) z-direction. The numbers in
parentheses denote the number of magnetic domains per moiré unit cell in each
layer.

In this study, we present a theoretical model to investigate the effects of
moiré patterns on twisted trilayer magnets (TTMs), as shown in Fig. 1a. In
comparison to previous studies on twisted bilayer magnets (TBMs)[24, 14, 15,
16, 18, 17, 23, 25, 26, 19], the TTM has complex magnetic domain structures
due to the four distinct coupling configurations of the local interlayer
interaction between the layers, as shown in Fig. 1b. By examining the
intriguing interlayer coupling effects in the moiré superlattice, we
demonstrate that the magnetic domain structures can be altered by adjusting
the twist angle and external magnetic field, as depicted in Fig. 1c. Our
findings provide a controllable platform for magnetic textures and enable
further exploration of the potential applications of moiré engineering in
multilayer magnetic systems.

##  II Noncollinear magnetism induced by moiré patterns

Figure 2: Local stacking patterns and interlayer exchange couplings. a Two
adjacent layers showing local stacking patterns that can be either
ferromagnetic (FM) or antiferromagnetic (AFM), with corresponding interlayer
exchange couplings. b Moiré unit cell with six distinct trilayer structures,
each composed of the bilayer stacking patterns shown in a. The first (second)
letters in each trilayer structure denote the bilayer stacking pattern between
the top and middle (middle and bottom) layers. c Six local stacking regions
are grouped into four patches with distinct interlayer coupling
characteristics (such as AFM-AFM), where the first (second) letters indicate
the coupling type of the interlayer exchange between the top-middle (middle-
bottom) layers.

The moiré superlattices of the TTMs are formed by rotating the middle layer of
an aligned trilayer honeycomb lattice around one of its hexagonal centers via
the commensurate angle
$\theta_{m,n}=\arccos{\left(\frac{m^{2}+n^{2}+4mn}{2m^{2}+2n^{2}+2mn}\right)}$,
where $m$ and $n$ are coprime integers [27] [see Fig. 1a]. The Heisenberg spin
model that describes the TTM is given by the following Hamiltonian,

$\displaystyle H$ $\displaystyle=\sum_{\textrm{l=t,m,b}}\left[-J\sum_{\langle
i,j\rangle}\mathbf{m}_{i}^{\textrm{l}}\cdot\mathbf{m}_{j}^{\textrm{l}}-D_{z}\sum_{i}\big{(}{m}_{z,i}^{\textrm{l}}\big{)}^{2}\right]$
$\displaystyle+\sum_{i,j}\Big{[}J_{ij}^{\textrm{tm}}\mathbf{m}_{i}^{\textrm{t}}\cdot\mathbf{m}_{j}^{\textrm{m}}+J_{ij}^{\textrm{mb}}\mathbf{m}_{i}^{\textrm{m}}\cdot\mathbf{m}_{j}^{\textrm{b}}\Big{]},$
(1)

where $\mathbf{m}_{i}^{\textrm{l}}$ represents the magnetization vector for
$i$-th spins in the top, middle, and bottom layers
($\textrm{l}=\textrm{t},\textrm{m},\textrm{b}$), respectively. $J>0$ is the
intralayer FM exchange interactions between nearest-neighbor spins.
$J_{ij}^{\textrm{l}_{1},\textrm{l}_{2}}$ represents the interlayer exchange
interactions between layer l1 and l2. $D_{z}>0$ indicates the single-ion
anisotropy energy.

In CrI3, the stacking-dependent interactions alternatively exhibit FM and AFM
interactions depending on the local stacking structures [28, 29]. To account
for this behavior, we use the following coupling function, derived from ab
initio calculations on laterally shifted bilayer CrI3 [18],

$\displaystyle J_{ij}^{\textrm{tm},\textrm{mb}}$
$\displaystyle=J_{0}\exp\big{(}-|r-d_{0}|/l_{0}\big{)}$
$\displaystyle+J_{1}^{s}\exp\big{(}-|r-r_{*}|/l_{1}^{s}\big{)}\sum_{a=1}^{3}\sin(\mathbf{q}_{a}^{s}\cdot\mathbf{r}_{ij})$
$\displaystyle+J_{1}^{c}\exp\big{(}-|r-r_{*}|/l_{1}^{c}\big{)}\sum_{a=1}^{3}\cos(\mathbf{q}_{a}^{c}\cdot\mathbf{r}_{ij}),$
(2)

where $\mathbf{r}_{ij}=\mathbf{r}_{i}-\mathbf{r}_{j}$ represents the in-plane
displacement between two interlayer spins, $d_{0}\approx 6.7\textrm{\AA}$ is
the out-of-plane distance, and $r=\sqrt{|\mathbf{r}_{ij}|^{2}+d_{0}^{2}}$. The
wave vectors
[$\mathbf{q}_{a}^{s,c}=q_{s,c}(\cos\varphi_{a}^{s,c},\sin\varphi_{a}^{s,c})$
for $a=1,2,3$, where $\varphi_{a}^{s}=(3\pi/2)(a-1)$ and
$\varphi_{a}^{c}=(3\pi/2)(a-1)+\pi/6$], and length parameters ($l_{0}$,
$l_{1}^{s,c}$) describe the leading oscillatory and decaying behavior of the
interlayer coupling, respectively. The energy scales ($J_{0},J_{1}^{s,c}$)
determine the overall strength of the coupling. The specific values of the
parameters are listed in Tab.1.

The local stacking between neighboring layers exhibits continuous modulation
between six distinct stacking configurations (AA, AB, BA, AB′, AC′,
AB${}_{1}^{\prime}$). According to ab-initio calculations, AA, AB, BA, and AB′
stacking exhibit FM character, while AC′ and AB${}_{1}^{\prime}$ exhibit AFM
character [See Fig.2a] [18]. In principle, the moiré superlattice of the TTM
is composed of combinations of the six stacking structures. However, the
stacking between the top-middle and the middle-bottom layers is not
independent since the inversion symmetry ties the stacking structures between
the top-middle and middle-bottom layers. As a result, the TTM possesses six
local stacking structures: AB-BA, BA-AB, AA-AA,
AB${}_{1}^{\prime}$-AB${}_{1}^{\prime}$, AC′-AB′, and AB′-AC′ [Fig. 2b], each
of which has FM-FM, FM-FM, FM-FM, AFM-AFM, AFM-FM, or FM-AFM couplings for the
top-middle and the middle-bottom layers, respectively. In turn, those local
stacking regions can be grouped into four distinct local patches according to
their interlayer coupling [See Fig.2c].

##  III Magnetic domain structure via twist engineering

Due to the coexistence of the AFM and FM interlayer couplings, the magnetic
ground state forms a magnetic domain in one of the layers where the domain has
an AFM spin configuration from its neighboring layer. In the TBM, the three
$C_{3z}$ symmetric magnetic domains reside in one of the two identical layers,
resulting in double degeneracy. In contrast, the TTMs possess multiple
possibilities with distinct domain energies depending on the location of the
domain. Specifically, the AFM-AFM patch can potentially stabilize the magnetic
domains in the middle layer. Alternatively, the magnetic domain can span
across a wider region, including the AFM-FM. FM-AFM patch, and the AFM-AFM
patch, in either the top or bottom layer. In Fig. 1c, we classify the four
possible distinct domain configurations, where the local interlayer coupling
and the domain configurations are consistent with each other.

$\textrm{J}_{0}(\textrm{meV})$ | $\textrm{J}^{\textrm{s}}_{1}(^{\prime\prime})$ | $\textrm{J}^{\textrm{c}}_{1}(^{\prime\prime})$ | $r_{*}(\textrm{\AA})$ | $l_{0}(^{\prime\prime})$ | $l^{\textrm{s}}_{1}(^{\prime\prime})$ | $l^{\textrm{c}}_{1}(^{\prime\prime})$ | $q_{s}(\textrm{\AA}^{-1})$ | $q_{c}(^{\prime\prime})$
---|---|---|---|---|---|---|---|---
$-0.1$ | $-0.5$ | $0.1$ | $7.3$ | $0.1$ | $0.3$ | $0.6$ | $0.7$ | $2.0$
Table 1: The parameters for the coupling function in Eq. (2), which was
derived from ab initio calculations on laterally shifted bilayer CrI3 [18].
Figure 3: Magnetic domain structures in twisted trilayer magnets. a Phase
diagram illustrating four magnetic phases: ferromagnetic (FM), noncollinear
domain (NCD), type-I magnetic domain (MD-I), and type-II magnetic domain (MD-
II) phases. The out-of-plane net magnetization
($M_{z}=\frac{1}{N}\sum_{i}m_{z,i}$) is shown as a function of twist angle
($\theta$) in the top (magenta), middle (blue), and bottom (green) layers, and
all layers (black), respectively. b Magnetic energy per spin relative to the
FM energy as a function of $\theta$. Each curve shows the intralayer gradient
energy ($E_{\textrm{intra}}$), interlayer coupling energy
($E_{\textrm{inter}}$), and total energy
($E_{\textrm{tot}}=E_{\textrm{intra}}+E_{\textrm{inter}}$), respectively. c
Magnetic configuration in the NCD phase, where the out-of-plane local
magnetization ($m_{z,i}$) is depicted using a color scale. The $3\times 3$
moiré superlattice unit cells are shown, with a green line indicating a single
unit cell. The second to fourth panels depict the magnified area (red line) in
each layer in the first panel, where the arrows denote the in-plane components
of the local magnetization. d Magnetic configuration in the MD-I phase. e
Magnetic configuration in the MD-II phase. f Phase diagram showing the
transition from the MD-I to MD-II phases as a function of intralayer exchange
($J$) and single-ion anisotropy ($D_{z}$) with the out-of-plane net
magnetization ($M_{z}$). For a-e, we use $J=2$ meV and $D_{z}=0.2$ meV. For
c-e, the values of $\theta$ are denoted by red markers in a and b. For f, we
use $\theta=2.45^{\circ}$ [$(m,n)=(13,14)$].

As the twist angle increases, the size of the moiré superlattice decreases.
The magnetic domain structure is determined by the competition between the
intralayer gradient energy ($E_{\textrm{intra}}$) and the interlayer coupling
energy ($E_{\textrm{inter}}$). In the case of small superlattices, the
intralayer gradient energy outweighs the interlayer coupling energy, resulting
in a higher energy cost for magnetic domain formation. Therefore, it becomes
energetically unfavorable for the system to form multiple magnetic domains. As
a result, at large twist angles ($\theta\gtrsim 4^{\circ}$), the ground state
exhibits a uniform FM order.

As the twist angle decreases, a continuous transition to the intermediate
noncollinear domain (NCD) phase occurs before the formation of the magnetic
domain. In the NCD phase, the spins in all layers are continuously tilted into
the $x-y$ plane, resulting in a decrease in the net magnetization:
$M_{z}=\frac{1}{3N}\sum_{\textrm{l=t,m,b}}\sum_{i=1}^{N}m_{z,i}^{\textrm{l}}$
[See Fig. 3a for the phase diagram]. This intermediate phase can be considered
a precursor to the formation of the magnetic domain and is characterized by
the continuous rotation of the magnetization vector, without a well-defined
magnetic domain. The presence of the NCD phase is related to the interlayer
exchange coupling between the top and middle layers, which competes with the
intralayer exchange coupling and leads to the rotation of the spins.

As the twist angle decreases further, the spin tilt angle increases, resulting
in a higher intralayer gradient energy $E_{\textrm{intra}}$ [Fig. 3b]. Despite
this increase in $E_{\textrm{intra}}$, the interlayer coupling energy
$E_{\textrm{inter}}$ results in a reduction of the total energy
$E_{\textrm{tot}}=E_{\textrm{intra}}+E_{\textrm{inter}}$, stabilizing the
magnetic domain (MD) phases where the spins in the layers flip inside the
domain. Among the four possible categorizations of the MDs [Fig. 1c], two are
energetically stabilized (MD-I and MD-II). In MD-I, a localized magnetic
domain is found in the AFM-AFM patch of the middle layer; in the MD-II phase,
the magnetic domains form in the top and bottom layers instead of the middle
layer, doubling the number of domains from $3$ to $6$ per unit cell [Fig. 3e].
The number and size of domains are increased and include the AFM-FM and FM-AFM
patches, as well as AFM-AFM patches, which in turn leads to a significant
increase in the intralayer gradient energy compared to that in the MD-I phase
[Fig. 3b]. However, the dominant interlayer coupling energy again results in a
reduction in the total energy, leading to the stabilization of these extensive
magnetic domains. Furthermore, this phase transition to the MD-II phase is a
first-order transition and involves an abrupt reconstruction of the domain
structure [Fig. 3a].

Figure 4: Manipulating the magnetic domain structures with external magnetic
fields. a Effect of the external magnetic fields ($B_{\textrm{ext}}$) in the
out-of-plane direction on a ferromagnetic (FM) ground state, transforming it
into a noncollinear domain (NCD) or an interlayer antiferromagnetic (I-AFM)
state. The net magnetization ($M_{z}$) is shown as a function of
$B_{\textrm{ext}}$, where $B_{\textrm{ext}}M_{z}>0$
($B_{\textrm{ext}}M_{z}<0$) indicates that the magnetic field is in the same
(opposite) direction as the net magnetization. b Transformation of an NCD
ground state into an FM or I-AFM state under the influence of a magnetic
field. c Transformation of a type-I magnetic domain (MD-I) ground state into
NCD or FM states. d Transformation of a type-II magnetic domain (MD-II) ground
state into any of the previously mentioned magnetic states. The schematic
images in each panel illustrate the corresponding magnetic structures. The
values of $\theta$ used for a-d are $6.01^{\circ}$, $3.89^{\circ}$,
$1.89^{\circ}$, and $1.30^{\circ}$, corresponding to
(m,n)=(5,6),(9,10),(17,18),and (25,26), respectively. The parameters
$J=2\mathrm{meV}$ and $D_{z}=0.2\mathrm{meV}$ are used for all plots.

##  IV Manipulation of domain structure by applying external magnetic fields

The magnetic domain structure of twisted trilayer magnets can be manipulated
by applying external magnetic fields, without the need for fine-tuning the
twist angle. Specifically, a uniform magnetic field applied in the out-of-
plane direction ($B_{\textrm{ext}}$) can induce transformations between
different magnetic states. To show this, we calculate stabilized magnetic
structures modified from ground states in each magnetic phase by including the
Zeeman coupling term,

$H_{Z}=-B_{\textrm{ext}}\sum_{\textrm{l=t,m,b}}\sum_{i}m_{z,i}^{\textrm{l}},$
(3)

in addition to the magnetic Hamiltonian in Eq. (1).

We manipulate each magnetic phase through the applications of magnetic fields.
First, as shown in Fig. 4a, an FM ground state can be transformed into an NCD
state by applying a magnetic field in the opposite direction as the net
magnetization ($M_{z}$). With increasing field strength, the NCD state can
further evolve into a metastable state without domains, called an interlayer
antiferromagnetic (I-AFM) state. In this state, the spins in the middle layer
flip toward the field direction, while the other layers remain unchanged. As
the field strength increases further, spins in the other layers also flip,
resulting in the stabilization of an FM state in the direction of the magnetic
field.

Similarly, an NCD ground state can be transformed into an FM or I-AFM state by
applying a magnetic field in the opposite direction [Fig. 4b]. Conversely, it
can be continuously transformed into an FM state by applying a magnetic field
in the same direction as the net magnetization.

Moreover, an MD-I ground state can be continuously (or discontinuously with
spin flipping) transformed into an NCD or FM state by applying a magnetic
field in the same (opposite) direction as the net magnetization [Fig. 4c].
However, in this case, unlike the FM and NCD ground states, the I-AFM state is
not stabilized.

Finally, an MD-II ground state can be transformed into an MD-I state by
applying a magnetic field in the same or opposite direction [Fig. 4d]. With
increasing field strength, the MD-I state can further change into an NCD or FM
state. As a result, all different magnetic structures of twisted trilayer
magnets can be achieved by simply applying a magnetic field to an MD-II ground
state, which offers wide flexibility in controlling the magnetic domain
structure without requiring fine-tuning of the twist angle.

##  V Discussion

Upon further investigation, we discovered that each magnetic phase could
exhibit additional magnetic domain structures with varying numbers of small
and large magnetic domains as excited states [see Fig. 5-7 in Appendix C], in
addition to the ground magnetic configurations [Fig. 3c-e]. These excited
states are notable for carrying skyrmion spin textures around the domain wall
[12, 14, 15, 13, 16, 17], which can be stabilized due to the domain structure
alone, without requiring other magnetic interactions that favor a winding spin
texture, such as the Dzyaloshinskii-Moriya interaction [30, 31, 32]. These
findings have significant implications since they demonstrate the possibility
of manipulating magnetic domain structures to stabilize new states with unique
spin textures, facilitating the development of novel spintronic devices with
enhanced functionalities [20, 21, 22].

In addition, our findings suggest that the magnetic domains in twisted
trilayer magnets can be easily manipulated via external magnetic fields,
eliminating the need for precise control over the twist angle, which can be
challenging in experimental techniques [9, 10, 11]. By applying a magnetic
field, we can induce a transition between the different magnetic domain
phases, enabling the precise control over the domain structure. Our study
provides a promising avenue for designing magnetic devices based on twisted
trilayer magnets, where the magnetic domains can be manipulated by external
magnetic fields, enabling greater flexibility and control.

###### Acknowledgements.

We acknowledge financial support from the Institute for Basic Science in the
Republic of Korea through the project IBS-R024-D1. This work was also
supported by the National Research Foundation of Korea(NRF) grant funded by
the Korea government(MSIT) (RS-2023-00252085,RS-2023-00218998). M.J.P. would
like to acknowledge the hospitality of the Center for Theoretical Physics of
Complex Systems during the visit. Lastly, we would like to express our
gratitude to Hee Chul Park for the helpful discussions.

## References

  * Balents  _et al._ [2020] L. Balents, C. R. Dean, D. K. Efetov, and A. F. Young, Superconductivity and strong correlations in moiré flat bands, Nature Physics 16, 725 (2020).
  * Andrei and MacDonald [2020] E. Y. Andrei and A. H. MacDonald, Graphene bilayers with a twist, Nature Materials 19, 1265 (2020).
  * Huang  _et al._ [2022] D. Huang, J. Choi, C.-K. Shih, and X. Li, Excitons in semiconductor moiré superlattices, Nature Nanotechnology 17, 227 (2022).
  * Chen  _et al._ [2020] X. Chen, X. Fan, L. Li, N. Zhang, Z. Niu, T. Guo, S. Xu, H. Xu, D. Wang, H. Zhang, A. S. McLeod, Z. Luo, Q. Lu, A. J. Millis, D. N. Basov, M. Liu, and C. Zeng, Moiré engineering of electronic phenomena in correlated oxides, Nature Physics 16, 631 (2020).
  * Kennes  _et al._ [2021] D. M. Kennes, M. Claassen, L. Xian, A. Georges, A. J. Millis, J. Hone, C. R. Dean, D. N. Basov, A. N. Pasupathy, and A. Rubio, Moiré heterostructures as a condensed-matter quantum simulator, Nature Physics 17, 155 (2021).
  * He  _et al._ [2021] F. He, Y. Zhou, Z. Ye, S.-H. Cho, J. Jeong, X. Meng, and Y. Wang, Moiré patterns in 2d materials: A review, ACS Nano 15, 5944 (2021).
  * Chen  _et al._ [2021] J. Chen, X. Lin, M. Chen, T. Low, H. Chen, and S. Dai, A perspective of twisted photonic structures, Applied Physics Letters 119, 240501 (2021), https://doi.org/10.1063/5.0070163 .
  * Rakib  _et al._ [2022] T. Rakib, P. Pochet, E. Ertekin, and H. T. Johnson, Moiré engineering in van der waals heterostructures, Journal of Applied Physics 132, 120901 (2022), https://doi.org/10.1063/5.0105405 .
  * Song  _et al._ [2021] T. Song, Q.-C. Sun, E. Anderson, C. Wang, J. Qian, T. Taniguchi, K. Watanabe, M. A. McGuire, R. Stöhr, D. Xiao, T. Cao, J. Wrachtrup, and X. Xu, Direct visualization of magnetic domains and moiré magnetism in twisted 2d magnets, Science 374, 1140 (2021), https://www.science.org/doi/pdf/10.1126/science.abj7478 .
  * Xu  _et al._ [2022] Y. Xu, A. Ray, Y.-T. Shao, S. Jiang, K. Lee, D. Weber, J. E. Goldberger, K. Watanabe, T. Taniguchi, D. A. Muller, K. F. Mak, and J. Shan, Coexisting ferromagnetic–antiferromagnetic state in twisted bilayer cri3, Nature Nanotechnology 17, 143 (2022).
  * Xie  _et al._ [2022] H. Xie, X. Luo, G. Ye, Z. Ye, H. Ge, S. H. Sung, E. Rennich, S. Yan, Y. Fu, S. Tian, H. Lei, R. Hovden, K. Sun, R. He, and L. Zhao, Twist engineering of the two-dimensional magnetism in double bilayer chromium triiodide homostructures, Nature Physics 18, 30 (2022).
  * Tong  _et al._ [2018] Q. Tong, F. Liu, J. Xiao, and W. Yao, Skyrmions in the moiré of van der waals 2d magnets, Nano Letters 18, 7194 (2018), pMID: 30285450, https://doi.org/10.1021/acs.nanolett.8b03315 .
  * Akram and Erten [2021] M. Akram and O. Erten, Skyrmions in twisted van der waals magnets, Phys. Rev. B 103, L140406 (2021).
  * Xiao  _et al._ [2021] F. Xiao, K. Chen, and Q. Tong, Magnetization textures in twisted bilayer $\mathrm{Cr}{\mathrm{x}}_{3}$ ($\mathrm{X}$=br, i), Phys. Rev. Research 3, 013027 (2021).
  * Akram  _et al._ [2021] M. Akram, H. LaBollita, D. Dey, J. Kapeghian, O. Erten, and A. S. Botana, Moiré skyrmions and chiral magnetic phases in twisted crx3 (x = i, br, and cl) bilayers, Nano Letters 21, 6633 (2021).
  * Ray and Das [2021] S. Ray and T. Das, Hierarchy of multi-order skyrmion phases in twisted magnetic bilayers, Phys. Rev. B 104, 014410 (2021).
  * Ghader  _et al._ [2022] D. Ghader, B. Jabakhanji, and A. Stroppa, Whirling interlayer fields as a source of stable topological order in moiré cri3, Communications Physics 5, 192 (2022).
  * Kim  _et al._ [2022] K.-M. Kim, D. H. Kiem, G. Bednik, M. J. Han, and M. J. Park, Theory of moire magnets and topological magnons: Applications to twisted bilayer cri3 (2022), arXiv:2206.05264 .
  * Hua  _et al._ [2023] C.-B. Hua, F. Xiao, Z.-R. Liu, J.-H. Sun, J.-H. Gao, C.-Z. Chen, Q. Tong, B. Zhou, and D.-H. Xu, Magnon corner states in twisted bilayer honeycomb magnets, Phys. Rev. B 107, L020404 (2023).
  * Hirohata  _et al._ [2020] A. Hirohata, K. Yamada, Y. Nakatani, I.-L. Prejbeanu, B. Diény, P. Pirro, and B. Hillebrands, Review on spintronics: Principles and device applications, Journal of Magnetism and Magnetic Materials 509, 166711 (2020).
  * Marrows and Zeissler [2021] C. H. Marrows and K. Zeissler, Perspective on skyrmion spintronics, Applied Physics Letters 119, 250502 (2021), https://doi.org/10.1063/5.0072735 .
  * Coey [2012] J. Coey,  _Magnetism and magnetic materials_ (Cambridge University Press, Cambridge, England, 2012).
  * Fumega and Lado [2023] A. O. Fumega and J. L. Lado, Moiré-driven multiferroic order in twisted crcl3, crbr3 and cri3 bilayers, 2D Materials 10, 025026 (2023).
  * Hejazi  _et al._ [2020] K. Hejazi, Z.-X. Luo, and L. Balents, Noncollinear phases in moiré magnets, Proceedings of the National Academy of Sciences 117, 10721 (2020), https://www.pnas.org/content/117/20/10721.full.pdf .
  * Wang  _et al._ [2020] C. Wang, Y. Gao, H. Lv, X. Xu, and D. Xiao, Stacking domain wall magnons in twisted van der waals magnets, Phys. Rev. Lett. 125, 247201 (2020).
  * Zheng [2023] F. Zheng, Magnetic skyrmion lattices in a novel 2d-twisted bilayer magnet, Advanced Functional Materials 33, 2206923 (2023), https://onlinelibrary.wiley.com/doi/pdf/10.1002/adfm.202206923 .
  * Moon and Koshino [2012] P. Moon and M. Koshino, Energy spectrum and quantum hall effect in twisted bilayer graphene, Phys. Rev. B 85, 195458 (2012).
  * 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 (2017).
  * Sivadas  _et al._ [2018] N. Sivadas, S. Okamoto, X. Xu, C. J. Fennie, and D. Xiao, Stacking-dependent magnetism in bilayer cri3, Nano Letters 18, 7658 (2018), pMID: 30408960, https://doi.org/10.1021/acs.nanolett.8b03321 .
  * Fert  _et al._ [2013] A. Fert, V. Cros, and J. Sampaio, Skyrmions on the track, Nature Nanotechnology 8, 152 (2013).
  * Nagaosa and Tokura [2013] N. Nagaosa and Y. Tokura, Topological properties and dynamics of magnetic skyrmions, Nature Nanotechnology 8, 899 (2013).
  * Fert  _et al._ [2017] A. Fert, N. Reyren, and V. Cros, Magnetic skyrmions: advances in physics and potential applications, Nature Reviews Materials 2, 17031 (2017).

##  Appendix A Iterative optimization method

To obtain a stable magnetic configuration, we iteratively update a trial
configuration by subtracting each magnetization vector $\mathbf{m}_{i}$ with
the gradient vector of the magnetic Hamiltonian $H$. The update rule is given
by,

$\mathbf{m}_{i}^{(s+1)}=\mathbf{m}_{i}^{(s)}-c\mathbf{h}_{i}^{(s)},$ (4)

where $\mathbf{h}_{i}\equiv\frac{\partial
H[{\mathbf{m}_{i}^{(s)}}]}{\partial\mathbf{m}_{i}^{(s)}}$ depends on the
magnetic configuration ${\mathbf{m}_{i}^{(s)}}$ in each iteration step $(s)$,
and $c$ controls the update rate. To preserve the norm of the magnetization
vectors, we normalize each vector after each update, i.e.,
$\mathbf{m}_{i}^{(s+1)}\rightarrow\mathbf{m}_{i}^{(s+1)}/\left|\mathbf{m}_{i}^{(s+1)}\right|$.
We use periodic boundary conditions during the iterative updates to ensure
consistency along the moiré unit cell boundary. We repeated this process more
than 50 times with different trial initial configurations to determine the
magnetic ground state, which is obtained by selecting the configuration with
the lowest energy.

##  Appendix B Estimation of intralayer gradient and interlayer coupling
energy

We calculated the intralayer gradient energy $E_{\textrm{intra}}$ and
interlayer coupling energy $E_{\textrm{inter}}$ presented in Fig. 3b of the
main text by using the following equations:

$\displaystyle E_{\textrm{intra}}$
$\displaystyle=\frac{1}{N}(E_{J}+E_{D_{z}}-E_{J}^{\textrm{FM}}-E_{D_{z}}^{\textrm{FM}}),$
(5) $\displaystyle E_{\textrm{inter}}$
$\displaystyle=\frac{1}{N}(E_{J_{\perp}}-E_{J_{\perp}}^{\textrm{FM}}),$ (6)

where $E_{J}$, $E_{D_{z}}$ and $E_{J_{\perp}}$ denote the magnetic energy
contribution from each magnetic interaction, defined as

$\displaystyle E_{J}$ $\displaystyle=-J\sum_{\textrm{l=t,m,b}}\sum_{\langle
i,j\rangle}\mathbf{m}_{i}^{\textrm{l}}\cdot\mathbf{m}_{j}^{\textrm{l}},$ (7)
$\displaystyle E_{D_{z}}$
$\displaystyle=-D_{z}\sum_{\textrm{l=t,m,b}}\sum_{i}\big{(}{m}_{z,i}^{\textrm{l}}\big{)}^{2},$
(8) $\displaystyle E_{J_{\perp}}$
$\displaystyle=\sum_{i,j}\Big{[}J_{ij}^{\textrm{tm}}\mathbf{m}_{i}^{\textrm{t}}\cdot\mathbf{m}_{j}^{\textrm{m}}+J_{ij}^{\textrm{mb}}\mathbf{m}_{i}^{\textrm{m}}\cdot\mathbf{m}_{j}^{\textrm{b}}\Big{]},$
(9)

and $E_{J}^{\textrm{FM}}$, $E_{D_{z}}^{\textrm{FM}}$, and
$E_{J_{\perp}}^{\textrm{FM}}$ are the corresponding magnetic energies for a
ferromagnetic state $\\{\mathbf{m}_{i}^{\textrm{l}}\\}=(0,0,1)$, respectively.
$N$ denotes the total number of spins in the moiré unit cell, which is given
by $N=\sum_{\textrm{l=t,m,b}}\sum_{i}$.

##  Appendix C Stabilized magnetic states with skyrmion spin textures

Figure 5: Stabilized magnetic states and skyrmion spin textures in the
noncollinear domain phase. a-c Out-of-plane magnetization (a), magnetic energy
per spin (b), and skyrmion number (c), where the state index denotes
stabilized magnetic states from different random initial configurations sorted
by energy. In a, the numbers in parentheses indicate the number of magnetic
domains in each layer (top, middle, and bottom, respectively). c Magnetic
configuration in the excited state with (0,1,0) domains, where the local
magnetization in the out-of-plane direction ($m_{z,i}$) is depicted by a color
scale. The $3\times 3$ moiré superlattice unit cells are shown, where the
green line denotes a single unit cell. The second to fourth panels show the
magnified area (red line) in each layer in the first panel, where the arrows
denote the in-plane components of the magnetization. e Magnetic configuration
in the excited state with (0,2,0) domains. f Magnetic configuration in the
excited state with (0,3,0) domains. For a-f, we use $J=2$ meV, $D_{z}=0.2$
meV, and $\theta=2.88^{\circ}$ [$(m,n)=(11,12)$].

Figure 5 illustrates stabilized magnetic states obtained from various random
initial configurations in the noncollinear domain (NCD) phase. In addition to
the ground state with no domains [(0,0,0) state], magnetic domain states
carrying different numbers of domains [(0,1,0), (0,2,0), (0,3,0)] are also
identified, each having a reduced net magnetization [Fig. 5a]. The creation of
domains increases the magnetic energy, indicating that they are metastable
states [Fig. 5b]. These magnetic domains also contain skyrmion spin textures
around their domain wall boundaries [Fig. 5c-f], where an increase in the
skyrmion number results in a slight rise in energy.

Figure 6: Stabilized magnetic states and skyrmion spin textures in the type-I
magnetic domain phase. a-c Out-of-plane magnetization (a), magnetic energy per
spin (b), and skyrmion number (c), where the state index denotes stabilized
magnetic states from different random initial configurations sorted by energy.
In a, the numbers in parentheses indicate the number of magnetic domains in
each layer (top, middle, and bottom, respectively). c Magnetic configuration
in the excited state with (0,3,0) domains, where the local magnetization in
the out-of-plane direction ($m_{z,i}$) is depicted by a color scale. The
$3\times 3$ moiré superlattice unit cells are shown, where the green line
denotes a single unit cell. The second to fourth panels show the magnified
area (red line) in each layer in the first panel, where the arrows denote the
in-plane components of the magnetization. e Magnetic configuration in the
excited state with (1,2,1) domains. f Magnetic configuration in the excited
state with (3,0,3) domains. For a-f, we use $J=2$ meV, $D_{z}=0.2$ meV, and
$\theta=2.13^{\circ}$ [$(m,n)=(15,16)$].

Similarly, Figure 6 illustrates stabilized magnetic states obtained from
various random initial configurations in the type-I magnetic domain (MD-I)
phase. In addition to the ground state with (0,3,0) domains, magnetic domain
states carrying different numbers of domains [(0,1,0), (0,2,0), (3,0,3)] are
identified, each having a reduced net magnetization [Fig. 6a]. The creation of
domains leads to an increase in magnetic energy, indicating that they are
metastable states [Fig. 6b]. Additionally, these magnetic domains contain
skyrmion spin textures around their domain wall boundaries [Fig. 6c-f], where
an increase in the skyrmion number results in a slight rise in energy.

Figure 7: Stabilized magnetic states and skyrmion spin textures in the type-II
magnetic domain phase. a-c Out-of-plane magnetization (a), magnetic energy per
spin (b), and skyrmion number (c), where the state index denotes stabilized
magnetic states from different random initial configurations sorted by energy.
In a, the numbers in parentheses indicate the number of magnetic domains in
each layer (top, middle, and bottom, respectively). c Magnetic configuration
in the excited state with (3,0,3) domains, where the local magnetization in
the out-of-plane direction ($m_{z,i}$) is depicted by a color scale. The
$3\times 3$ moiré superlattice unit cells are shown, where the green line
denotes a single unit cell. The second to fourth panels show the magnified
area (red line) in each layer in the first panel, where the arrows denote the
in-plane components of the magnetization. e Magnetic configuration in the
excited state with (2,1,2) domains. f Magnetic configuration in the excited
state with (1,2,1) domains. For a-f, we use $J=2$ meV, $D_{z}=0.2$ meV, and
$\theta=1.20^{\circ}$ [$(m,n)=(27,28)$].

Figure 7 illustrates stabilized magnetic states obtained from various random
initial configurations in the type-II magnetic domain (MD-II) phase. In
addition to the ground state with (3,0,3) domains, the magnetic domain states
carrying different numbers of domains [(2,1,2), (1,2,1), (0,3,0)] are also
identified, each having a reduced net magnetization [Fig. 7a]. The creation of
domains leads to an increase in magnetic energy, indicating that they are
metastable states [Fig. 7b]. Additionally, these magnetic domains contain
skyrmion spin textures around their domain wall boundaries [Fig. 7c-f], where
an increase in the skyrmion number results in a slight rise in energy.