Blog

Spin Hamiltonians

We adopt a generic Hamiltonian for twisted \({Cr}{I}_{3}\) bilayers including exchange, DM, and Kitaev interactions,

The vector \({{{{{{\bf{S}}}}}}}_{{li}}\) denotes the classical spin on a site i in layer l, with position vector \({{{{{{\bf{r}}}}}}}_{{li}}\). We set \(l={{{{\mathrm{1,2}}}}}\) for the bottom and top layers, respectively. \(J\) and \({{{{{\mathscr{A}}}}}}\) are the NN intralayer Heisenberg coupling and the easy axis magnetic anisotropy, respectively. The third term in \({{{{{\mathscr{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}_{\perp }({{{{{{\bf{r}}}}}}}_{{ij}})={J}_{\perp }({{{{{{\bf{r}}}}}}}_{2i}-{{{{{{\bf{r}}}}}}}_{1j})\) is the distance-dependent interlayer coupling between spins \({{{{{{\bf{S}}}}}}}_{2i}\) and \({{{{{{\bf{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,68,69, parametrized by K and Γ. All terms in \({{{{{\mathscr{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.⁶¹. The triplet \((\lambda ,\mu ,\nu )\) in the summation represents any permutation of the octahedral coordinates.

Table 1 reports the proper parametrization of \({{{{{\mathscr{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).

The atomistic interlayer coupling approach

\({Cr}{I}_{3}\) is the most prominent example of stacking-dependent magnetism in 2D magnets and constitutes an excellent candidate to discover the intriguing physics underlying moiré magnets. In the twisted bilayer, the \({Cr}{I}_{3}\) moiré superlattice hosts regions with local rhombohedral (AB, BA, and AA) and monoclinic (\({{{{{\mathscr{M}}}}}}\)) stackings (Fig. 1a). The interlayer exchange is FM (AFM) at the rhombohedral (monoclinic) stackings^{70,71,72,73,74}. Several authors calculated the stacking-dependent interlayer exchange energy in CrI₃^53,57,72. Here, we adopt the DFT results obtained by Xiao et al.⁵³ 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.⁵³), characterized by three AFM monoclinic regions labeled I, II, and III. Without loss of generality, we choose a spin \({{{{{{\bf{S}}}}}}}_{2i}\) at position \({{{{{{\bf{r}}}}}}}_{2i}\) in layer 2 (top layer). The spatially modulated effective interlayer coupling can be expressed as⁵³\({J}_{\perp }^{{eff}}({{{{{{\bf{r}}}}}}}_{2i})=\mathop{\sum}\limits_{j}{J}_{\perp }({{{{{{\bf{r}}}}}}}_{2i}-{{{{{{\bf{r}}}}}}}_{1j})=-{E}_{{int}}({{{{{{\bf{r}}}}}}}_{2i})/4{S}^{2}\). Here, \(S=3/2\) is the spin of the \({Cr}\) magnetic atom. Next, we assume \({J}_{\perp }({{{{{{\bf{r}}}}}}}_{2i}-{{{{{{\bf{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 \({\delta }_{2i}\) (and consequently \({J}_{\perp }({{{{{{\bf{r}}}}}}}_{2i}-{{{{{{\bf{r}}}}}}}_{1j})\)) can be determined numerically for every spin \({{{{{{\bf{S}}}}}}}_{2i}\) in the moiré supercell by solving the equation,

In particular, we solve Eq. 4 for a large cutoff interlayer interaction radius \(\left|{{{{{{\bf{r}}}}}}}_{2i}-{{{{{{\bf{r}}}}}}}_{1j}\right|\le a\) (a is the \({Cr}{I}_{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 distance-dependent interlayer coupling \({J}_{\perp }\) 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 \(\mathop{\sum}\limits_{j}{J}_{\perp }({{{{{{\bf{r}}}}}}}_{2i}-{{{{{{\bf{r}}}}}}}_{1j}){{{{{{\bf{S}}}}}}}_{1j}\). Similarly, for a specific site j in layer 1, the moiré field is \(\mathop{\sum}\limits_{i}{J}_{\perp }({{{{{{\bf{r}}}}}}}_{2i}-{{{{{{\bf{r}}}}}}}_{1j}){{{{{{\bf{S}}}}}}}_{2i}\).

Finally, we note that Eq. 3 assumes an isotropic form for \({J}_{\perp }({{{{{{\bf{r}}}}}}}_{2i}-{{{{{{\bf{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⁷⁵ 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é \({Cr}{I}_{3}\) compared to the continuum approach. The input files for the Vampire simulations were prepared using the Mathematica software⁷⁶.

The sLLG equation reads

where γ and λ are the gyromagnetic ratio and Gilbert damping, respectively. \({{{{{{\bf{H}}}}}}}_{{li}}^{{eff}}\) is the net magnetic field on spin \({{{{{{\bf{S}}}}}}}_{{li}}\),

The first term in \({{{{{{\bf{H}}}}}}}_{{li}}^{{eff}}\) can be derived from Eqs. (2), (3), and (4), while \({{{{{{\bf{H}}}}}}}_{{li}}^{{th}}\) is the effective thermal field included in Vampire using the Langevin Dynamics method⁷⁷.

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\times {10}^{7}\) total time steps. We use \(1\times {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⁷⁷ 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⁷⁵.

We analyzed the ground state spin textures at commensurate twist angles in the range \(0.65^\circ \le \theta \le 6^\circ\), which is relevant to the recent experimental work⁵⁸. 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é \({Cr}{I}_{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 \({{{{{{\bf{S}}}}}}}_{1j}\) along \(+\hat{{{{{{\bf{z}}}}}}}\) and yields a collinear interlayer field on the top layer, aligned along \(\pm \hat{{{{{{\bf{z}}}}}}}\). Consequently, the trivial interlayer field in the monolayer approximation does not favor the formation of stable TSTs in the absence of the NN DMI⁵³. 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 \({{{{{{\bf{S}}}}}}}_{2i}\), \(\mathop{\sum}\limits_{j}{J}_{\perp }({{{{{{\bf{r}}}}}}}_{2i}-{{{{{{\bf{r}}}}}}}_{1j}){{{{{{\bf{S}}}}}}}_{1j}\), is determined by the sign of \({J}_{\perp }({{{{{{\bf{r}}}}}}}_{2i}-{{{{{{\bf{r}}}}}}}_{1j})\) and the directions of the contributing spins \({{{{{{\bf{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 \(\theta =1.35^\circ\) 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 \(+\hat{{{{{{\bf{z}}}}}}}\).

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^\circ \le \theta \le 6^\circ\) and external magnetic fields in the range \(100{mT}\le B\le 1.5{T}\), varied in steps of \(100{mT}\). We additionally included the magnetic fields \(0.25{T}\), \(0.75{T}\), and \(1.25{T}\). 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 \(\left|Q\right| \, > \, 1\) (Fig. 2e, f). The TSTs are observed in the range \(\theta \, \lesssim \, 2.13^\circ\) 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 bi-skyrmions and bi-antiskyrmions observed in nonchiral frustrated magnetic films^8,29,30,31, but with zero separation between the antiskyrmions³¹.

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 \(\theta =1.35^\circ\) and \(B=0.25{T}\). 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., \(\theta =2.13^\circ\)) 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⁷⁸.

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 \(\theta \, \lesssim \, 3.15^\circ\).

The MBs in the chiral/nonchiral Heisenberg and Heisenberg-IDMI models are found to disappear at relatively large angles (\(\theta \, \gtrsim \, 4.2^\circ\)), 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 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.⁵³, 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éel-type 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^\circ \le \theta \le 6^\circ\).

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 (\(\theta \, \lesssim \, 2.44^\circ\)), 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 \(\left|Q\right| > 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 \({Ni}{I}_{2}\)³⁸. On the other hand, a sizeable NN DMI \((d=0.15{meV})\) locks the chirality and induces spontaneous Néel-type skyrmions for \(\theta \, \lesssim \, 3.15^\circ\) 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^\circ \le \theta \le 6^\circ\). 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.⁵⁷, 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.⁵⁷ 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 \, < \, \theta \, \lesssim \, 3.15^\circ\). For a short-range above the minimal angles, Akram et al.⁵⁷ 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⁵⁷ 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 could generate merged skyrmions or magnetic strips at any angle \(0.65^\circ \, < \, \theta \, \lesssim \, 2.44^\circ\), 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^\circ \le \theta \le 6^\circ\), 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 (\(\sim \! 15 \, K\)). As a test, we added a sizable easy-axis anisotropy (\({{{{{\mathscr{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⁵⁸.