Whirling interlayer fields as a source of stable topological order in moiré CrI3

The moiré engineering of two-dimensional magnets opens unprecedented opportunities to design novel magnetic states with promises for spintronic device applications. The possibility of stabilizing skyrmions in these materials without chiral spin-orbit couplings or dipolar interactions is yet to be explored. Here, we investigate the formation and control of ground state topological spin textures (TSTs) in moiré CrI3\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${Cr}{I}_{3}$$\end{document} using stochastic Landau–Lifshitz–Gilbert simulations. We unveil the emergence of interlayer vortex and antivortex Heisenberg exchange fields, stabilizing spontaneous and field-assisted ground state TSTs with various topologies. The developed study accounts for the full bilayer spin dynamics, thermal fluctuations, and intrinsic spin-orbit couplings. By examining the effect of the Kitaev interaction and the next nearest-neighbor Dzyaloshinskii–Moriya interaction, we propose the latter as the unique spin-orbit coupling mechanism compatible with experiments on monolayer and twisted CrI3\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${Cr}{I}_{3}$$\end{document}. Our findings contribute to the current knowledge about moiré skyrmionics and uncover the nature of spin-orbit coupling in CrI3\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${Cr}{I}_{3}$$\end{document}. The reduced dimensions of 2D magnets expose their magnetic anisotropy, adding a twist into the system provides another degree of freedom to explore. Here, the authors use stochastic Landau–Lifshitz-Gilbert simulations to investigate ground-state topological spin textures induced by interlayer fields in twisted bilayer CrI3.

T opological magnetism is an intriguing field on the frontier of condensed matter physics with great promises for future information technology [1][2][3][4] . Skyrmion 5 , a particle-like spinwhirling vortex, was the first class of topological spin textures (TSTs) to be realized experimentally in the chiral magnet MnSi 6 . While skyrmions remain a prominent example of TSTs, several alternatives have been predicted and observed in recent years, such as antiskyrmions 7 , biskyrmions 8,9 , skyrmioniums 10,11 , and bimerons 12,13 . Generally, TSTs are robustly stable with particle-like properties due to their topological protection. They carry a quantized topological charge and can interact via attractive and repulsive forces. The topological charge quantifies the real-space nontrivial topology of the spins, and it is defined as [14][15][16] where SðrÞ denotes the spin density field. In Eq. 1, the term ∂ x S ∂ y S represents the vorticity of the spin texture and is determined by the in-plane components of the boundary spins. The vorticity classifies out-of-plane spin textures into skyrmions and antiskyrmions, with vortex and antivortex profiles, respectively 17 . Based on their whirling profile, skyrmions can be further classified as Bloch-type and Néel-type. Precisely, in a Bloch-type (Néel-type) skyrmion, the spins rotate perpendicular to (along) the radial direction when moving from the core to the periphery. Other significant properties of the TST morphology are the polarity (alignment of the core spins) and helicity (the global rotation angle around the out-of-plane axis). Chiral interactions in magnetic films stabilize TSTs with fixed chirality (fixed vorticity and helicity) 6,[17][18][19][20][21][22][23] . The most prominent example of chiral interactions is the NN DMI, arising in noncentrosymmetric lattices lacking space-inversion symmetry. Moreover, TSTs can form in centrosymmetric materials due to nonchiral interactions 8,9,[24][25][26][27][28][29][30][31][32][33] , such as dipolar interactions 24,25,29 , magnetic anisotropy 28,32 , quantum fluctuations 34 , and static random fields 32,33 . In the nonchiral magnetic films, the vorticity and helicity act as additional degrees of freedom, relevant for spintronic and topological applications 8,[29][30][31] . An external magnetic field usually assists the formation of TSTs in chiral and nonchiral magnetic films. On the other hand, materials hosting spontaneous TSTs are rare. So far, theoretical research has predicted spontaneous TSTs in itinerant magnets with high-order spin interactions 19,21,[35][36][37] and NiI 2 monolayers 38 with anisotropic exchange interaction.
Furthermore, 2D magnets offer unique opportunities to engineer TSTs via the modulated interlayer coupling in the moiré superlattices of twisted or mismatched bilayers. The mechanism was initially discussed in a heterostructure formed of a nonchiral ferromagnetic (FM) monolayer on an antiferromagnetic (AFM) substrate 51 . Bloch-type skyrmions emerge from the registrydependent interlayer exchange and dipolar couplings in the heterostructure. Several theoretical studies followed, exploring TSTs in mismatched and twisted magnetic bilayers [52][53][54][55][56] . Akram et al. 54 and Hejazi et al. 55 reported skyrmions in the chiral FM-AFM heterostructure. Ray et al. 56 explored the TSTs in nonchiral moiré magnets with interlayer dipolar interactions. Moreover, Xiao et al. 53 predicted metastable skyrmions as excited states in twisted chromium trihalides CrX 3 bilayers ðX ¼ I; BrÞ without NN DMI or dipolar interactions. The effect of NN DMI in twisted CrX 3 ðX ¼ I; Br; ClÞ was addressed in a recent work 57 , reporting Néeltype skyrmions in CrI 3 and CrBr 3 .
The theoretical predictions on moiré magnets are reinforced by recent experiments on twisted CrI 3 bilayers 58,59 . Xu et al. 58 reported coexisting FM-AFM states in twisted bilayer CrI 3 at small twist angles (1 ≲ θ ≲ 4 ). The FM (AFM) domains form at the local rhombohedral (monoclinic) stackings in the moiré supercell. However, the AFM domains are found to disappear at relatively larger twist angles (θ ≳ 4 ), resulting in a pure FM ground state. On the other hand, Song et al. 59 targeted slightly twisted CrI 3 bilayers with θ ≲ 0:5 . Due to the lattice reconstruction at tiny twist angles, the authors observed a substantially disordered moiré pattern and a complicated magnetic structure.
From the theoretical side, the methods used so far to determine the magnetic textures in moiré magnets are the Landau-Lifshitz-Gilbert (LLG) approach 51,53,54,58 , the continuum low-energy field theory 52,55 , and Monte Carlo simulations with continuum interlayer exchange 56 . Nevertheless, an atomistic theoretical model that incorporates thermal fluctuations has not been reported yet. Moreover, existing studies on moiré magnets (except Akram et al. 57 ) oversimplify the moiré bilayer as a monolayer with a spatially modulated external magnetic field. This approximation assumes fixed spins in one of the layers. Akram et al. 57 included the spin dynamics in both layers but with an interlayer interaction taken only between nearest neighbors. Additionally, previous studies on moiré CrI 3 52,53,56,57 adopt the Heisenberg model, excluding spin-orbit interactions like the intrinsic next NN DMI (IDMI) and the Kitaev interaction. However, the pure Heisenberg model is inadequate for CrI 3 since it fails to explain the gapped magnon spectrum [60][61][62] in this material. Meanwhile, both the Heisenberg-IDMI 61,62 and Heisenberg-Kitaev 60,61 models can reproduce the spectrum, creating wide controversy on the true microscopic spin Hamiltonian for CrI 3 63 . Hence, it is important to investigate the Heisenberg-IDMI and Heisenberg-Kitaev models for moiré CrI 3 and check if they can reproduce the recent experimental observations 58 .
In this work, we predict the emergence of whirling interlayer exchange fields in twisted moiré CrI 3 . At long moiré periodicity, the interlayer interaction dominates, and the spins align with the nontrivial moiré fields, inducing stable TSTs. Hence, the whirling moiré fields uncover a mechanism to stabilize ground state TSTs without the need for NN DMI or dipolar interactions.
The presented study mimics realistic experiments by cooling the bilayer system from initial paramagnetic states. We use the stochastic Landau-Lifshitz-Gilbert (sLLG) equation to account for thermal effects. We include both layers in the spin dynamics simulations after developing an atomistic approach that accounts for the effective interlayer coupling beyond the NN approximation. Further, the study covers and compares the chiral and nonchiral Heisenberg, Heisenberg-IDMI, and Heisenberg-Kitaev models.
In the absence of the chiral NN DMI, the vorticity and helicity of the interlayer fields act as degrees of freedom. As a result, the stochastic time evolution of the interlayer fields does not follow deterministic rules in the nonchiral models. Instead, the interlayer fields' profiles at low temperatures depend crucially on the initial paramagnetic state. As a result, nontrivial moiré fields with various chiralities can emerge at small twist angles to stabilize a zoo of TSTs. On the contrary, a sizeable chiral NN DMI locks the moiré interlayer field's chirality, resulting in Néel-type skyrmions for any initial paramagnetic state.
The study explores the field-assisted and spontaneous formation of the TSTs. Cooling with an external magnetic field traps the TSTs at the monoclinic AFM regions of the moiré supercell. Conversely, spontaneous TSTs can emerge in the FM regions, AFM regions, or a combination of FM and AFM regions. Moreover, spontaneous spin textures can merge to form magnetic strips with chiral domain walls. Both the spontaneous and the field-assisted TSTs can be drastically manipulated by a new magnetic field applied at 0 K.
Finally, we show that moiré engineering provides insights into the fundamental interactions underlying 2D magnets. Specifically, we consider the recent experimental results 58 on moiré CrI 3 to address the controversy regarding the correct microscopic model for CrI 3 . We explore the twist and temperature (T) dependence of the averaged magnetization (M) in the Heisenberg-IDMI and Heisenberg-Kitaev models. The Heisenberg-IDMI model is found to display a twist-dependent M À T curve that evolves towards an FM ground state at large angles (θ ≳ 4:3 ), in agreement with the experimental results 58 . However, the Heisenberg-Kitaev model fails to describe the twist-dependent ground state. Therefore, we conclude that the Heisenberg-IDMI model is the unique model that can reproduce the experimental results in monolayer and moiré CrI 3 .

Results
Spin Hamiltonians. We adopt a generic Hamiltonian for twisted CrI 3 bilayers including exchange, DM, and Kitaev interactions, The vector S li denotes the classical spin on a site i in layer l, with position vector r li . We set l ¼ 1; 2 for the bottom and top layers, respectively. J and A are the NN intralayer Heisenberg coupling and the easy axis magnetic anisotropy, respectively. The third term in H is the nonchiral Néel-type IDMI 61,62,[64][65][66][67] between next NN illustrated in Supplementary Fig. 1a. In contrast, the fourth term accounts for possible NN chiral DMI ( Supplementary Fig. 1b) induced by the broken inversion symmetry in the twisted system. J ? ðr ij Þ ¼ J ? ðr 2i À r 1j Þ is the distance-dependent interlayer coupling between spins S 2i and S 1i . The sixth term is the Zeeman coupling by an external magnetic field B normal to the bilayer. The last term is the Kitaev Hamiltonian 60,67-69 , parametrized by K and Γ. All terms in H are expressed in the coordinate axes illustrated in Fig. 1a, except the Kitaev Hamiltonian, which is expressed in the octahedral coordinate axes 60,61 . Details on the octahedral coordinate axes can be found, for example, in Chen et al. 61 . The triplet ðλ; μ; νÞ in the summation represents any permutation of the octahedral coordinates. Table 1 reports the proper parametrization of H to reproduce the Heisenberg, Heisenberg-IDMI, and Heisenberg-Kitaev models. Further, we will denote models with sizeable (respectively negligible) NN DMI as chiral (nonchiral).
Here, we adopt the DFT results obtained by Xiao et al. 53 and develop a method to determine the atomistic interlayer exchange coupling accordingly. This approach will later allow us to simulate the time evolution of the interlayer exchange fields. Figure 1b presents the moiré interlayer exchange energy E int (Xiao et al. 53 ), characterized by three AFM monoclinic regions labeled I, II, and III. Without loss of generality, we choose a spin S 2i at position r 2i in layer 2 (top layer). The spatially modulated effective interlayer coupling can be expressed as 53 Here, S ¼ 3=2 is the spin of the Cr magnetic atom. Next, we assume J ? ðr 2i À r 1j Þ decays exponentially 54,55 as a function of the distance, that is where r 0 denotes the interlayer separation. Then, the decay factor δ 2i (and consequently J ? ðr 2i À r 1j Þ) can be determined numerically for every spin S 2i in the moiré supercell by solving the equation, In particular, we solve Eq. 4 for a large cutoff interlayer interaction radius r 2i À r 1j ≤ a (a is the CrI 3 lattice constant) to ensure adequate distribution of the effective interlayer interaction over the interlayer sites, while the exponential term cancels irrelevant contributions. This is desired to avoid biased interlayer fields in the atomistic sLLG simulation. The numerical approach eventually yields the distancedependent interlayer coupling J ? for all interacting spins in layers 1 and 2. The effective moiré interlayer field on a specific site i in layer 2 can be deduced as ∑ j J ? ðr 2i À r 1j ÞS 1j . Similarly, for a specific site j in layer 1, the moiré field is ∑ i J ? ðr 2i À r 1j ÞS 2i . Finally, we note that Eq. 3 assumes an isotropic form for J ? ðr 2i À r 1j Þ. Further improvement to include anisotropic terms such as an interlayer spin-orbit coupling might be interesting, requiring future DFT investigations.
The stochastic LLG approach. We used the Vampire atomistic spin dynamics software 75 to determine the time evolution of the spin textures using the sLLG equation. Generally, atomistic simulations better account for the complex stacking-dependent magnetism in moiré CrI 3 compared to the continuum approach. The input files for the Vampire simulations were prepared using the Mathematica software 76 .
The sLLG equation reads where γ and λ are the gyromagnetic ratio and Gilbert damping, respectively. H eff li is the net magnetic field on spin S li , The first term in H eff li can be derived from Eqs. (2), (3), and (4), while H th li is the effective thermal field included in Vampire using the Langevin Dynamics method 77 .
To mimic real experiments, we launch the sLLG simulations with random initial spins at a high temperature (50 K), followed by gradual cooling to 0 K. The system is cooled with 3 10 7 total time steps. We use 1 10 À16 s for the time step and a cooling time of 1 ns. The spin configurations are determined in both layers at different temperatures, down to the ground state (0 K), using the Heun integration scheme 77 and imposing periodic boundary conditions. The Heun integration scheme is preferred for stochastic spin dynamics due to its computational efficiency and accuracy in reproducing the correct ground state 75 .
We analyzed the ground state spin textures at commensurate twist angles in the range 0:65 ≤ θ ≤ 6 , which is relevant to the recent experimental work 58 . We will focus more on the Heisenberg and the Heisenberg-IDMI models since the Heisenberg-Kitaev model was found inconsistent with the experimental results on moiré CrI 3 .
Cooling with an applied magnetic field. We start the discussion by considering the field-assisted TSTs in the nonchiral Heisenberg and Heisenberg-IDMI models. For slight twists, the interlayer interaction dominates the intralayer exchange 52,58 and induces three magnetic bubbles in the monoclinic AFM regions of the moiré. Consequently, simulating the time evolution of the interlayer exchange fields is crucial to investigate the possible emergence of stable TSTs.
The monolayer approximation freezes S 1j along þẑ and yields a collinear interlayer field on the top layer, aligned along ±ẑ. Consequently, the trivial interlayer field in the monolayer approximation does not favor the formation of stable TSTs in the absence of the NN DMI 53 . Here, we reveal a different picture where the moiré interlayer fields acquire nontrivial profiles, stabilizing ground state TSTs.
In our approach, the orientation of the interlayer field on S 2i , ∑ j J ? ðr 2i À r 1j ÞS 1j , is determined by the sign of J ? ðr 2i À r 1j Þ and the directions of the contributing spins S 1j . Similarly, for the interlayer field on layer 1. Consequently, in the stochastic description of the moiré spin dynamics, the interlayer fields act as dynamic (time-dependent) random fields during the cooling process before reaching static configurations near 0 K (Supplementary Movie 1). The thermal fluctuations dominate the magnetic interactions down to the ordering temperature and play an important role in the time evolution of the moiré interlayer fields. However, the thermal fluctuations diminish at low temperatures, and the competing magnetic interactions gradually dominate below the ordering temperature. The combined effects of the thermal fluctuations and the competing magnetic interactions shape the moiré field during the cooling process, and the moiré fields can converge to nontrivial textures near 0 K. Figure 1c illustrates an example of antivortex interlayer fields emerging in the AFM regions of the moiré. The figure simulates the moiré field on the bottom layer of the nonchiral Heisenberg model with θ ¼ 1:35 and B ¼ 1T. Since the interlayer interaction dominates the intralayer exchange at long moiré periodicity, the spins align with the interlayer field to minimize the energy. Consequently, the interlayer field shapes the spin textures' morphology in the AFM regions. Specifically, the spins inherit the vorticity and helicity of the antivortex interlayer field (Fig. 1c), stabilizing three ground state antiskyrmions (Fig. 1d).
In the nonchiral Heisenberg and Heisenberg-IDMI models, the vorticity and helicity are degrees of freedom. The time evolution of the interlayer fields does not follow deterministic rules that can predict the stochastic simulation's results a priori. The profile of the final interlayer fields and the corresponding magnetic ground state depend crucially on the initial paramagnetic state, chosen randomly in our simulations. Any combination of the degrees of freedom is allowed and can be probed by different initial states. Accordingly, topological and trivial MBs can coexist in the moiré bilayer ( Supplementary Fig. 2a, b), with arbitrary distribution over the layers. For completeness, examples of trivial interlayer fields stabilizing non-topological MBs are presented in Supplementary Fig. 2c, d. Generally, the spins interacting with a topological or trivial MB in the opposite layer are (slightly) twisted relative to þẑ.
The crucial dependence of the ground state on the initial paramagnetic configuration is elaborated in Supplementary Fig. 3. Therefore, a reliable study requires exhaustive numerical experiments, including several initial random spin configurations. We studied bilayers with commensurate angles 0:65 ≤ θ ≤ 6 and external magnetic fields in the range 100mT ≤ B ≤ 1:5T, varied in steps of 100mT. We additionally included the magnetic fields 0:25T, 0:75T, and 1:25T. We tested six distinct random initial states for each twist angle and magnetic field, generating 108 different simulations for a given twist angle. Similar to previous theoretical works 53,56,57 , we present results for simulations over a single moiré supercell. Nevertheless, we have extensively investigated samples with multiple moiré supercells. We found that including additional moiré supercells in the simulation of field-assisted TSTs is equivalent to changing the initial random spin configurations on a single moiré supercell. Meanwhile, we stress that the moiré-periodicity of the magnetic ground state is broken in multi-moiré supercell samples. In particular, since adjacent moiré supercells have distinct initial random spin configurations, they converge to ground states with different types of MBs (topological or trivial) in the AFM regions.
The low-temperature interlayer fields manifest in various profiles, inducing multi-flavored TSTs, such as antiskyrmions (Fig. 1c, d), Néel-type skyrmions (Fig. 2a, b), Bloch-type skyrmions (Fig. 2c, d), and topological defects with Q j j> 1 (Fig. 2e, f). The TSTs are observed in the range θ ≲ 2:13 as a rough estimation. The high topological charge (Q ¼ 2) spin texture in Fig. 2e can be interpreted as a magnetically stable bound state of two antiskyrmions with opposite helicities. Bound states of two skyrmions with opposite helicities (Q ¼ À2) are also possible in the moiré bilayer, and an example is presented in Supplementary Fig. 4e. The formation of such skyrmionic molecules stems from the helicity degree of freedom and can be realized only in nonchiral magnets. They are analogous to biskyrmions and bi-antiskyrmions observed in nonchiral frustrated magnetic films 8,29-31 , but with zero separation between the antiskyrmions 31 .
The rich spectrum of TSTs in the nonchiral Heisenberg and Heisenberg-IDMI models is further elaborated in Table 2. In particular, we present the ground state with the maximum number of TSTs from six trials performed for each angle and magnetic field. For example, we choose to present the result of simulation 1 from Supplementary Fig. 3 because it displays three TSTs for the Heisenberg model with θ ¼ 1:35 and B ¼ 0:25T. This criterion is occasionally dropped in Table 2 to report TSTs with a high topological charge. Therefore, the results of Table 2 are to be interpreted as insightful rather than deterministic results.
It can be noticed from Table 2 that TSTs can be realized even at relatively large angles (e.g., θ ¼ 2:13 ) by varying the magnetic field. Moreover, simulations based on the nonchiral Heisenberg and Heisenberg-IDMI models yield different results for the same initial random state, indicating that the IDMI affects the spin dynamics in the bilayer system. Nevertheless, the overall results and topological spectra are comparable for the nonchiral Heisenberg and Heisenberg-IDMI models, suggesting that the IDMI interaction does not have characteristic signatures on the morphology of the TSTs.
We verified the robustness of the TSTs when the external magnetic field is turned off at 0 K and observed stability in the TSTs' morphology with no effect on the topological charge ( Supplementary Fig. 4a-d). Consequently, the textured interlayer interaction stabilizes the topological order without a permanent external field, which is desired for skyrmion-based spintronic devices. In experiments, adjacent moiré supercells will have distinct initial states and converge to different ground states. Nevertheless, our results suggest that successive heating-cooling trials and magnetic field manipulation can experimentally establish a topologically rich ground state, with skyrmions, antiskyrmions, high topological charge spin textures, and a minimal number of trivial magnetic bubbles. Such ground states with coexisting distinct TSTs are challenging to realize in conventional materials and are in-demand for memory and logic applications 78 .
Unsurprisingly, a sizeable NN DMI locks the chirality of the TSTs in the Heisenberg and the Heisenberg-IDMI models, hence producing deterministic results. Specifically, the chiral models present Néel-type skyrmions with a fixed topological charge Q ¼ À1 for any initial paramagnetic spin configuration (Supplementary Fig. 1c, d). The Néel-type skyrmions can form in any layer, and the moiré-periodicity is broken in multi-moiré supercell samples. Our calculations assumed a sizeable NN DMI d ¼ 0:15 meV, stabilizing Néel-type skyrmions in the range θ ≲ 3:15 .
The MBs in the chiral/nonchiral Heisenberg and Heisenberg-IDMI models are found to disappear at relatively large angles (θ ≳ 4:2 ), and the bilayer converges systematically towards an FM state as we increase the twist angle. Further discussions are presented in the last part of the Results section.
We proceed to discuss the possibility of manipulating the topological ground state at 0 K. Applying a new magnetic field in the direction of the MBs' core spins promotes a controllable outward motion of the domain walls without affecting the topological charge. Consequently, it is possible to inflate and couple the spin textures initially trapped in the AFM regions to realize skyrmion pairs ( Supplementary Fig. 5a, b), antiskyrmion pairs ( Supplementary Fig. 5d, e), and TST-(trivial) MB pairs ( Supplementary Fig. 5g, h). Moreover, a sufficiently large Numerical values for the Hamiltonian's parameters in the Heisenberg, Heisenberg-IDMI, and Heisenberg-Kitaev models for a monolayer CrI 3 . IDMI stands for the intrinsic Dzyaloshinskii-Moriya interaction with strength D. The parameters J and A are the intralayer Heisenberg coupling and the easy axis magnetic anisotropy, respectively. K and Γ are parameters that describe the Kitaev interaction. We note that the nearest-neighbor Dzyaloshinskii-Moriya interaction parameter d is zero in the monolayer case due to the centrosymmetric symmetry of the honeycomb lattice.
magnetic field inflates the MBs drastically and eventually induces a global magnetization reversal ( Supplementary Fig. 5c, f, i and Supplementary Movie 2). As a result, a final topological ground state is achieved with reversed spins compared to the initial state. Generally, the MBs jump to the opposite layer in the final ground state, and the spin reversal can modify the TST's vorticity and helicity ( Supplementary Fig. 5). We note the related discussion in Xiao et al. 53 , based on the monolayer approximation that cannot capture the complete picture presented here. In particular, the magnetization reversal erases the MBs in the monolayer approach since they cannot form in the opposite layer (assumed with fixed spins). The stabilization mechanism for TSTs via whirling moiré interlayer fields is general and applies to the nonchiral Heisenberg-Kitaev model. The TSTs' charges in the Heisenberg-Kitaev model are comparable to the previous models (Table 2). However, unlike the IDMI, the Kitaev interaction has clear signatures on the ground state TSTs' chirality. The Kitaev interaction induces a canting-like effect, where neighboring peripheral spins can cross to form pairs of frustrated spins (Supplementary Fig. 6a-d). Therefore, the nonchiral Heisenberg-Kitaev model displays TSTs with unconventional morphologies compared to the Heisenberg and the Heisenberg-IDMI models. The canting-like effect persists in the chiral Heisenberg-Kitaev model, and ideal Néeltype skyrmions were not observed even for a sizeable NN DMI ( Supplementary Fig. 6e, f).
The Heisenberg-Kitaev model does not transfer to an FM bilayer, and the MBs survive at large angles. The in-plane Heisenberg interaction is weak (Table 1) and is not expected to dominate the local AFM coupling even at large angles. Therefore, the observed behavior suggests that the interlayer AFM interaction is dominant over the Kitaev in-plane interaction throughout the range 0:65 ≤ θ ≤ 6 .
Cooling without an applied field. Whirling interlayer fields can also stabilize spontaneous TSTs without an external magnetic field in the three models. Unlike the field-assisted TSTs, which are confined to the AFM regions, the spontaneous TSTs can form in any region of the moiré supercell. In particular, the spontaneous TSTs can emerge in the AFM regions, FM regions, or a combination of the AFM and FM regions, depending on the initial spin configuration (Fig. 3). The merging of spin textures over the AFM and FM regions can generate giant spontaneous TSTs at slight twists ( Supplementary Fig. 7). Moreover, the spin textures can merge across the entire moiré supercell (Fig. 3) to form magnetic strips in one or both layers. These wavy strip-shaped magnetic domains are separated by chiral domain walls and can develop in different directions depending on the particularly merged MBs. Figure 3 illustrates the crucial dependence of the spontaneous spin textures on the initial random configuration. Therefore, similar to the discussion in the previous section, the moiréperiodicity of the spontaneous spin textures is broken since adjacent moiré supercells have different initial states. Nevertheless, the single-moiré supercell simulations remain faithful to reproduce the main features of the spontaneous spin textures. Our intensive numerical investigation of multi-moiré supercell samples did not disclose substantially new information, except that the moiré-periodicity is broken in such samples.
In the chiral/nonchiral Heisenberg and Heisenberg-IDMI models, the various merging scenarios are promoted at small angles (θ ≲ 2:44 ), whereas the spin textures are confined to the AFM regions at larger angles. Moreover, our previous discussion regarding the degrees of freedom in the nonchiral models remains valid for spontaneous TSTs. As a result, intriguing merged TSTs profiles can form at small angles in the nonchiral Heisenberg and Heisenberg-IDMI models, like skyrmion -(trivial) MBs (Fig. 4b, d), antiskyrmion -(trivial) MBs (Fig. 4e), skyrmion--antiskyrmion pairs, and (anti)skyrmion clusters with high topological charges ( Fig. 4a and Table 3). The Q j j>1 spontaneous TSTs are allowed by the helicity degree of freedom in the nonchiral models. For example, the Q ¼ 3 TST in Fig. 4a constitutes a bound state of three skyrmions with oppositely swirling spins. Such spontaneous skyrmionic clusters have been predicted only in itinerant magnets 37,79 and semiconducting NiI 2 38 . On the other hand, a sizeable NN DMI ðd ¼ 0:15meVÞ locks the chirality and induces spontaneous Néel-type skyrmions for θ ≲ 3:15 in the chiral Heisenberg and Heisenberg-IDMI models (Table 3 and Supplementary Fig. 8).
The spontaneous TSTs can be drastically tuned at 0 K. A magnetic field applied opposite to the core spins decouples the merged spin textures and confines them back to the AFM regions (Fig. 4b, c; e, f; Supplementary Movie 3). Conversely, the magnetization can be reversed by a magnetic field applied parallel to the core spins, trapping the reversed state's MBs in the AFM regions of the opposite layer (Supplementary Movie 4).
The merging and manipulation scenarios described above are also observed in the chiral/nonchiral Heisenberg-Kitaev models. The merging remains possible in these models throughout the inspected twist angle range 0:65 ≤ θ ≤ 6 . Moreover, the Kitaev interaction induces a canting-like effect for the spontaneous TSTs ( Supplementary Fig. 10), similar to their field-assisted counterparts.
Further insights on the spontaneous TSTs in all models are presented in Table 3. The results are selected from six simulations with distinct initial random states for each twist angle. We applied the same criterion as Table 2, with a preference for TSTs not trapped in the AFM regions.
Finally, we compare our methods and results with Akram et al. 57 , who studied the TSTs in the chiral Heisenberg model, neglecting the thermal effects and adopting a continuum approximation of the bilayer Hamiltonian. The authors used a NN interlayer approach and performed LLG simulations at 0 K starting from various initial FM states. The present work includes the thermal effects and uses an atomistic Hamiltonian. Further, we present an approach for the interlayer field beyond the NN and initiate the sLLG simulations from random spins. The differences in results and conclusions are summarized below.
Akram et al. 57 predicted three topological magnetic phases and concluded firm rules to determine the magnetic phase as a function of the twist angle. The first phase appears only at minimal angles, with Néel-type skyrmions confined to the AFM regions. However, in our study, simulations from various initial random states showed that this phase could be realized in the chiral Heisenberg model at any angle 0:65 < θ ≲ 3:15 . For a short-range above the minimal angles, Akram et al. 57 reported a second phase with a single merged-skyrmion scenario. In particular, one skyrmion is formed over the three AFM regions in a layer, leaving the opposite layer FM. This merge avoids the FM regions and differs substantially from the various merged skyrmions reported in our study (Table 3). At relatively larger angles, the previous work 57 revealed a third magnetic phase with strips in one of the layers and a Néel-type skyrmion trapped in an AFM region of the opposite layer. In our work, the Néel-type skyrmion in this phase forms over a combination of AFM and FM regions. Moreover, we demonstrated that distinct initial states Table 2 Extended results on the field-assisted topological spin textures (TSTs).

Model
Magnetic Field (T) The table elaborates the rich spectrum of field-assisted TSTs in the nonchiral models, presenting selected results from several simulations with distinct random initial states as described in the text. The symbols θ and B stand for the twist angle and the magnitude of the external magnetic field, respectively. We have calculated the topological charge Q from Eq. (1), using the numerical algorithm detailed by M. Augustin et al. 50 (see their supplementary materials for full details). The subscript of Q specifies the antiferromagnetic region to which the TST is confined. The superscripts B and T stand for the bottom (l ¼ 1) and top (l ¼ 2) layers, respectively. IDMI stands for the intrinsic Dzyaloshinskii-Moriya interaction. In the Heisenberg and the Heisenberg-IDMI models, the symbol next to the topological charge indicates the TST's nature. We used the following conventions: antiskyrmion*, Bloch-type skyrmion**, Néel-type skyrmion • , antiskyrmion pair •• , and skyrmion pair • *. In the Heisenberg-Kitaev model, the TSTs are identified as vortices ðQ < 0Þ and antivortices ðQ > 0Þ. The dashes indicate the absence of TSTs.
could generate merged skyrmions or magnetic strips at any angle 0:65 < θ ≲ 2:44 , without firm rules. Further, we observed additional magnetic phases not captured previously, such as magnetic strips in both layers and skyrmions trapped in the FM regions.
The twist-dependent averaged magnetization. The sLLG approach offers valuable insights into the variation of the averaged magnetization with the temperature and the twist angle. Figure 5a and b present the M À T curves for selected twist angles in the nonchiral Heisenberg and Heisenberg-IDMI models, respectively. The M À T curves are comparable in the two models. The ordering temperature is virtually independent of the twist (Fig. 5a, b), with a value near 25 K. The ground state averaged magnetization (at T ¼ 0 K) varies smoothly with the twist, and the moiré bilayer approaches a pure FM ground state at large angles. This behavior follows the gradual alignment of the MB's spins along the positive z-axis in the nonchiral Heisenberg and Heisenberg-IDMI models. Above 4.3°, all spins acquire positive S z components (Fig. 5d, e) and the MBs disappear from the moiré superlattice. However, the normalized averaged magnetization remains slightly below unity (Fig. 5a, b) due to residual tilted spins near the cores of the AFM regions.
The Heisenberg-Kitaev model shows fundamentally different behavior. In particular, the MBs persist throughout the range 0:65 ≤ θ ≤ 6 , leaving the M À T curve almost independent of the twist (Fig. 5c, f). Consequently, the ground state averaged magnetization does not approach the FM limit. Moreover, the Heisenberg-Kitaev model reveals a lower ordering temperature ($ 15 K). As a test, we added a sizable easy-axis anisotropy (A ¼ 0:15 meV) to the Heisenberg-Kitaev model, and we did not observe significantly different results. We conclude that the Heisenberg-Kitaev model cannot reproduce the experimentally observed twist-dependent magnetic ground state 58 .

Discussion
Technological implementation of TSTs crucially depends on discovering new topological magnetic materials and novel mechanisms for their stabilization. Moiré magnets are ideal candidates in this direction, which justifies the current tremendous interest in their fundamental and applied physics. Indeed, moiré skyrmionics is still in its early stages, promising vast opportunities for impactful discoveries. The recent experiments 58,59 on moiré CrI 3 constitute a significant advancement towards the experimental observation of TSTs in moiré magnets, which would require accurate mapping of the directions of spins in the moiré superlattice.
Theoretical models of moiré CrI 3 with sizeable NN DMI or dipolar interactions are expected to produce skyrmionic structures since these interactions are conventional sources for TSTs.
Nevertheless, the dipolar interactions might be negligibly weak 53 , while a significant NN DMI is not guaranteed and awaits future DFT or experimental confirmation. Consequently, it is essential to discover sources of topological magnetic textures that emerge exclusively from the moiré magnetism and go beyond the conventional skyrmion sources. This will pass through theoretical modeling that minimizes the approximations, conjugated with simulations that include the most relevant effects. The TST extends over regions II, III, and the AB-stacked ferromagnetic (FM) region. It hosts more than 2800 spins, carries a topological charge Q ¼ 3, and represents a bound state of three skyrmions with different helicities. b The nonchiral Heisenberg-IDMI model with θ ¼ 2:28 hosts a TST combining a Bloch-type skyrmion (region II) and a topologically trivial magnetic bubble (MB) in region III. Applying a magnetic field À0:5ẑT at 0 K separates the MB and the skyrmion (c) while changing the skyrmion's morphology to a Néel-type skyrmion. d, e TSTs in the Heisenberg model for θ ¼ 2 (d) and θ ¼ 2:13 (e). They correspond to a Néel-type skyrmion merged with a topologically trivial MB (d), and an antiskyrmion merged with two topologically trivial MBs (e). f The merged spin configuration in (e) can be detached by a magnetic field À0:5ẑT at 0K. The decoupling inverts the helicity of the antiskyrmion confined to region I.
Selected spontaneous ground state TSTs in the chiral and nonchiral versions of the three models. The criterion for selecting the results is described in the text. IDMI stands for the intrinsic Dzyaloshinskii-Moriya interaction. The subscript of the topological charge Q specifies the (merged) regions over which the TST is localized. I À III stands for the monoclinic regions I, II, and III of the moiré supercell.
In our study, we investigated the ground state TSTs in moiré CrI 3 , developing a study that accounts for the thermal effects, presents an atomistic approach for the interlayer coupling, and includes the IDM and Kitaev interactions on top of the Heisenberg in-plane exchange. We uncovered a stabilization mechanism for topological magnetic textures emerging from moiré interlayer fields with nontrivial textures. At large moiré periodicity, the whirling interlayer fields stabilize various types of ground state TSTs in the nonchiral models. Including the spin dynamics in both layers is crucial to account for this stabilization mechanism. Further, we showed that the monolayer approximation does not accurately describe the ground state manipulation.
The extra degrees of freedom (vorticity and helicity) characterizing the TSTs in nonchiral magnetic films attracted significant attention due to their technological relevance. These degrees of freedom are present in the nonchiral models of the moiré CrI 3 , which allow the formation of skyrmionic clusters with high topological charges. A magnetic field is required temporarily during the cooling process to trap the TSTs in the AFM regions. Moreover, moiré CrI 3 broadens the class of materials that can host spontaneous TSTs. The moiré interlayer field constitutes a source for spontaneous TSTs, similar to the high-order spin interactions in itinerant magnets and the anisotropic exchange in NiI 2 38 .
The Heisenberg model cannot describe the magnetic excitations in monolayer CrI 3 . Accordingly, the pure Heisenberg model is not suitable to simulate the moiré magnetism in CrI 3 . We have tested the Heisenberg-IDMI and the Heisenberg-Kitaev models and concluded that only the Heisenberg-IDMI is consistent with the twist-dependent ground state, observed experimentally 58 in moiré CrI 3 . Therefore, our study suggests that the Heisenberg-IDMI is the correct model for CrI 3 , in agreement with a very recent experimental study on topological magnons in monolayer CrI 3 80 . Note that our study does not account for the experimentally observed lattice relaxation at tiny twist angles θ ≲ 0:5 , which are excluded in the present study. Moreover, the lattice relaxation is found to be negligible 58 in experimental samples with θ ≳ 1 . Further investigations are required to determine the relevance of lattice relaxation in the twist angle range 0:5 < θ < 1 .
Beyond CrI 3 , our results suggest that Kitaev magnets might constitute better candidates for moiré skyrmionics than Heisenberg magnets. We demonstrated that a large Kitaev interaction could not dominate the stacking-dependent interlayer interaction. As a result, moiré Kitaev magnets can support nontrivial ground states over a broader twist angle range than Heisenberg magnets. These observations motivate attention towards Kitaev magnets for further advancement in moiré skyrmionics.

Methods
All methods are included in the Results section.

Data availability
The data that support the findings of this study are available from the corresponding author upon reasonable request. Supplementary Movies 1-4 are available in the supplementary information. The twisted bilayers are polarized by a magnetic field B ¼ 1ẑT. In (a) and (b), the averaged magnetization approaches the pure FM ground state in the large twist angle interval θ ≳ 4:3 . However, in (c), the averaged magnetization is almost independent of the twist angle. d, e The minimal S z (z component of the spin) value in the moiré supercell as a function of the twist angle for the Heisenberg (d) and the Heisenberg-IDMI (e) models at T ¼ 0K. A steep increase in S z is observed in the interval θ ≳ 3 . At θ ≳ 4:3 , the magnetic bubbles (MBs) disappear in both models. f Persistent ground state MBs in the top layer of the nonchiral Heisenberg-Kitaev model at θ ¼ 6 .