Triple-meron crystal in high-spin Kitaev magnets

Spin textures with nontrivial topology hold great promise in future spintronics applications since they are robust against local deformations. The meron, as one of such spin textures, is widely believed to appear in pairs due to its topological equivalence to a half skyrmion. Motivated by recent progresses in high-spin Kitaev magnets, here we investigate numerically a classical Kitaev-$\Gamma$ model with a single-ion anisotropy. An exotic spin texture including three merons is discovered. Such a state features a peculiar property with an odd number of merons in one magnetic unit cell and it can induce the topological Hall effect.Therefore, these merons cannot be dissociated from skyrmions as reported in the literature and a general mechanism for such a deconfinement phenomenon calls for further studies. Our work demonstrates that high-spin Kitaev magnets can host robust unconventional spin textures and thus they offer a versatile platform not only for exploring exotic states in spintronics but also for understanding the deconfinement mechanism in the condensed-matter physics and the field theory.

Introduction. -Topological spin textures (TSTs), which have attracted enormous attention in condensed matter physics, can be described by the homotopy theory [1] with a non-trivial mapping π n (S m ), where n and m are the corresponding degrees of freedom in the real and the parameter spaces, respectively. A topological charge Q thus counts the number of times of the real space configuration covering the parameter space. Among these TSTs, the skyrmions (π 2 (S 2 )), originally proposed in the particle physics [2], have become the focus of recent research in magnetism due to their potential applications in spintronics [3][4][5][6][7][8]. Experimentally, they have already been successfully discovered in noncentrosymmetric chiral magnets [9][10][11][12][13], magnetic heterostructures [14][15][16][17] and centrosymmetric magnets [18][19][20]. In a skyrmion, Q is an integer which is given by [1] where S = S(x, y) is the unit spin vector at the position (x, y). In addition, recent studies show that, given a sufficiently large effective easy-plane anisotropy, a skyrmion can be transferred into a meron-antimeron pair [21][22][23][24]. Such meron and antimeron carry one half skyrmion topological charge, and thus are treated as half skyrmions. From the spin configurations, the most distinct difference between skyrmions and merons is that at the perimeter the spins of merons lie within the plane while those of skyrmions point out of the plane. On the other hand, ever since the discovery of the exact solution of the Kitaev honeycomb model [25], great efforts have been made to synthesize materials to realize such a model. These materials, so called Kitaev magnets, are two-dimensional transition-metal compounds with strong spin-orbit coupling [26]. The promising candidates include Na 2 IrO 3 , Li 2 IrO 3 and α-RuCl 3 [27,28]. In these compounds, honeycomb-located cations are surrounded by the edge-sharing octahedral anions. Strong spin-orbit coupling entangles spin and orbital degrees of freedom together, resulting in an effective spin-1/2 model with the Heisenberg (J) and Kitaev (KS γ i S γ j ) interactions. Further research suggests that a bond-dependent symmetric off-diagonal Γ and/or Γ interactions should also be taken into account [29][30][31]. This JKΓΓ model and its relevance have been proposed [32] to explain the possible quantum spin liquids observed in experiments.
However, the Kitaev interaction is not limited to spin-1/2 compounds. Recent numerical and theoretical analysis [33][34][35] suggest that it may also exist in the van der Waals materials CrI 3 , CrGeTe 3 and CrSiTe 3 . Interestingly, the crystal structures of these compounds bear a strong resemblance to those of the honeycomb iridates but Cr owns an effective S = 3/2 spin. Similar to its spin-1/2 counterparts these compounds can also be described by the JKΓΓ model but possibly with an additional single-ion anisotropy (A). Although the sought-after quantum spin liquids are not preferable in these high-spin Kitaev magnets [36], stable TSTs originating from frustration [37][38][39][40][41] are highly desirable. This calls for extensive and immediate studies to search for the TSTs in these high-spin Kitaev magnets. However, so far, most works focus simply on the magnetic The two-dimensional honeycomb lattice lies in the ab plane and c is perpendicular to the plane. The xyz coordinate systems are along the three Cr-Te bonds. In the abc coordinate system, x, y and z are given by respectively. The Kitaev interactions on the bonds in blue, green and red correspond to x, y and z Ising type, respectively. orders [42][43][44] but relevant topological phenomena are rarely touched. We fill in this gap by investigating the classical KΓA model numerically. A triple-meron crystal (TmX) with three (anti)merons in one magnetic unit cell (MUC) is discovered, suggesting the mechanism for generating (anti)merons is beyond our present knowledge. Furthermore, we show that such a TmX can cause the topological Hall effect (THE) [4] which is absent in the meron-antimeron crystals [3].
Model. -The crystal structures [33,[47][48][49] of the Crbased compounds mentioned above are sketched in Fig. 1. Due to their similarity, here we exemplify them by the CrGeTe 3 . In CrGeTe 3 six Te atoms are at the corners of an octahedron and the Cr atom sits at the center. These edge-shared octahedrons form honeycomb layers stacking along the c axis, which are stabilized by a weak interlayer van der Waals coupling. Due to the symmetry of the crystal, two coordinate systems [33], abc and xyz are involved in our later discussion (Fig. 1). In principle, J, K, Γ, Γ and A terms are all allowed by the lattice symmetry. However, for a given compound, some terms may be zero. For example, the density-functional theory calculation shows that in CrGeTe 3 the Heisenberg term J becomes zero at some compressive strain [50]. Here, for simplicity, we consider the KΓA model, which is given by where i, j means that the summation runs over all nearest neighbors. S γ i = S i · γ( γ = x, y, z), e.g., the γcomponent of the spin vector at site i associated with the γ bonds on the honeycomb lattice and α, β are other two components. This Hamiltonian can be parameterized as K = sinθ cos φ, Γ = sin θ sin φ and A = cos θ, with θ ∈ [0, π] and φ ∈ [0, 2π). Now the model (2) depends on only two parameters θ and φ. Furthermore, one can see that the transformation φ → φ + π is equivalent to flipping all the spins in one sublattice. Therefore, we only need to consider the parameters in the region φ ∈ [0, π).
In our Monte Carlo simulations, the maximum size is 2 × 36 × 36 and the temperature ranges from 0.005 to 0.2. The ground state is then obtained by optimizing the Monte Carlo data. Various sizes are used to confirm we obtain the correct MUC (SM [55] Sec. II). The groundstate energy, the spin configurations and the static spin structure factor are calculated to map out the phase diagram. In the following, all the figures are plotted in the abc coordinate system.
Results. -In Fig. 2(a), we present the ground-state phase diagram of the KΓA model in the region θ ∈ [0.45π, 0.60π], φ ∈ [0.5π, 1.0π] (see SM [55] Sec. VI for the full phase diagram). Of all these phases, the one in red marked by TmX is the most exciting discovery in our work. In the following, we will focus on this phase and present the details. For simplicity, we choose one representative point in the TmX phase, i.e., θ = 0.515π, φ = 0.68π to demonstrate its unusual properties. In Fig. 2(b), we plot the ground-state spin configuration at the representative point. One can see there are 18 spins in one MUC. Moreover, these 18 spins are divided into three hexagons with their edges shared, which are marked by A, B and C, respectively. After a simple inspection one finds the spin at the core of each hexagon points out of the plane, while the edge spins lie almost in the plane. This is the characteristic feature of (anti)merons, indicating that there are three (anti)merons in one MUC. This finding is very distinct from our previous knowledge that the meron and antimeron should emerge in pairs [21][22][23][24]. In a common sense, the meron and antimeron with a topological charge ∓1/2 can be mapped [21] onto the southern and northern hemispheres, respectively, and wrap them once. Hence, considering a large enough easy-plane interaction, such as an in-plane anisotropy or the dipoledipole interaction, a skyrmion could be dissociated into two halves from the equator, transforming into a meronantimeron pair [21]. However, as shown in Fig. 2(b), there are three Bloch-type (anti)merons. Each one is composed of ten spins, including one polarized core spin, its three nearest neighbours (NNs) and six next-nearest neighbours (NNNs). All the core spins locate on the same sublattice, and are completely polarized along the c axis (down for A and C, up for B). The three NNs of each (anti)meron locate on the other sublattice and have the same polar angle. In particular, in B and C they are on the opposite hemisphere to the core spin. At the edges, the equator is formed by six NNNs locating on the same sublattice as the core spin. All NNNs lie almost in the ab plane. These characters suggest that B and C are of the AFM type. We want to mention that in our model FM merons are readily available by the transformation φ → φ + π, equivalent to flipping all spins on one sublattice.
One could notice that, if swirling around one meron in any direction, its NN and NNN spins rotate in the same direction while the core spins in B and C are antiparallel. In addition, for one loop around the meron, the spins in A rotate 4π, wrapping the hemisphere twice. By contrast, the spins wrap the hemisphere once in B and C, suggesting the topological charge of A should be twice of those of B and C. In order to confirm our analysis, we calculate the topological charges of the three merons at the representative parameter point following the definition given by Berg and Lüscher [2] (SM [55] Sec. III), which turn out to be Q A = 1.055231, Q B = 0.472850 and Q C = −0.528081, indicating the existence of a meron (Q C ) and antimeron (Q B ) pair and a high-Q antimeron (Q A ). Particularly, the high-Q antimeron is topologically equivalent to a half high-Q antiskyrmion [57]. It is known theoretically that the topological charge for a ideal (high-Q) (anti)meron should be exactly n/2 with n an integer. However, in the lattice system, the three (anti)merons share the outmost spins as the equator. When frustrations disturb the outmost edge slightly off the ab plane, we can expect that the topological charge of each (anti)meron may deviate from its theoretical value accordingly. Since the increase or decrease of solid angles should occur synchronously and one can expect that the summation of Q A , Q B and Q C should be an integer. Indeed, this is what we have observed from our results, i.e., it is 1.0. Particularly, in the TmX phase there are points with their Q A , Q B and Q C agreeing well with their theoretical values and the summation giving 1.0 as well.
The crystals formed by magnetic spin textures with a nonzero Q often exhibit nontrivial transport properties. The prototypical example is the THE [4]. However, it has been shown [3] that there is no THE in the meronantimeron crystals due to the vanishing net Q. Thus, the TmX provides an opportunity to explore the THE in the meron crystals. To show this, we study the doubleexchange model, which describes itinerant electrons coupled with the spin texture [4], where c † iσ (c iσ ) is the creation (annihilation) operator of electrons with the spin σ at the site i. t is the hopping integral between two nearest neighbors i and j. J is the coupling strength between the electron and on-site magnetization S i , and σ denotes the Pauli matrix. Since each unit cell of the TmX carries a finite Q, it provides a finite magnetic flux and leads to the THE.
The nth band structure E n (q) and the corresponding eigenvector φ n (q) can be easily obtained by diagonalizing the Hamiltonian (3) in the strong coupling limit [4] (see the SM [55] Sec. IV for details). We plot the band structure along the high-symmetry lines in Fig 3(a). One may notice that there is a Dirac-cone structure with a certain large gap in the vicinity of Γ point between the bands n = 6 and n = 7. With the E n (q) and φ n (q) available, some important quantities such as the Berry curvature, the Chern number C n and the topological Hall conductance σ xy are readily obtained. In Fig 3(b), we plot σ xy given by the Kubo formula [4] (see SM [55] Sec. IV for other quantities), where Ξ is the size of the system and f is the Fermi distribution function. One may notice that at zero temperature the Hall conductance is quantized to be −1 in the band gap between n = 6 and n = 7 where E f ∈ [−0.73, −0.61], which indicates a finite total Chern number, i.e., En<E f C n = σ xy /(e 2 /h).
As shown in Fig. 2, the TmX phase and the 18 B phase are separated by the line θ = π/2, where A = 0. In both phases, there are 18 spins in one MUC. On the line θ = π/2 (A = 0), these two states are degenerate (SM [55] Sec. I). One may notice that in the TmX phase A < 0, which is thus of the easy-axis type. In the literature, several mechanisms have been proposed to stabilize merons, such as the dipole-dipole interaction [58], the easy-plane anisotropy [5,21] as well as their interplay [59], the relative phase shifts among multiple helical spin density waves [3], and the interplay between the biquadratic interaction and the Dzyaloshinskii-Moriya interaction [60]. One question naturally arises as to why the TmX phase is stable under the easy-axis anisotropy. Actually, the TmX state is rooted in the Γ model, which has highly degenerate ground states [1]. Its degeneracy is partly lifted by the Kitaev term, resulting in the degenerate TmX phase and 18 B phase (see SM [55] Sec. I and III). Moreover, as shown in Fig. 4 2 /N is different in these two states. Therefore, in the presence of a finite A, the degeneracy is lifted. Particularly, S 2 c is larger in the TmX state, and hence when A < 0 the TmX state is energetically favored.
Summary and Outlook. -In this work, we report our discovery of the TmX state in the high-spin Kitaev magnets and the topological effect led by such a state is demonstrated. Although the Kitaev interaction was first proposed in the spin-1/2 Kitaev honeycomb model [25], recent studies show that it may also exist in highspin magnets. This offers great opportunities to study TSTs in such high-spin Kitaev magnets. Our discovery of the TmX state represents a fantastic progress in this direction and laid the foundation for the forthcoming Kitaev spintronics. On the other hand, the merons, as the classical solutions of Yang-Mills equation, were proposed as the mechanism for the quark confinement. Owing to their one half topological charge, they can only exist in pairs [62][63][64]. This assessment is also supported by numerous studies in magnets [21][22][23][24] and photonic systems [65]. Our findings are obviously beyond this paradigm and thus extend the scope of the deconfinement [66], which leads to our conjecture that merons with a half-integer Q should be in pairs while those with an integer Q can appear solely. Further studies, with Kitaev magnets as a feasible starting point, are necessary to explore such a deconfinement phenomenon.
We are grateful to Shi-Zeng Lin, Yan Zhou, Meng Xiao and Yalei Lu for helpful discussions. This work is supported by the National Natural Science Foundation Supplementary Material on "Triple-meron crystal in high-spin Kitaev magnets" In this Supplementary Material, we present more details of the numerical methods, the corresponding results and the complete phase diagram. The energy optimization method and how to use it to determine the phase transition points are shown in Sec. I. In Sec. II we show how to determine the magnetic unit cell numerically. In Sec. III, we explain the origin of the TmX phase. In Sec. IV, the details for calculating the topological charges and the topological Hall effect caused by the TmX are provided. In Sec. V, we show results from the atomistic dynamics simulations and in Sec. VI we provide the complete phase diagram. The phase diagram is determined by the ground-state energy, the spin configurations and the static spin structure factors.
Sec. I: The ground-state energy and phase-transition points by the Energy Optimization Method Based on the spin configuration obtained by the Monte Carlo simulation, we could further use an energy optimization method to target the ground-state energy with a controllable precision. In some cases, we can derive an analytic expression for the ground-state energy. Here, we take the TmX as an example to illustrate this method.
The paradigmatic spin configuration in the TmX state, which includes three kinds of merons, is shown in Fig. 2 of the main text. As presented in Fig. sm-1, there are eighteen spins in one magnetic unit cell. For each meron, the core spin and its three NNs can be written as (S c |S i S j S k ) for short. In the abc coordinate system, the classical spin at site i can be parameterized as where ϑ i ∈ [0, π] and ϕ i ∈ [0, 2π). Pertaining to the core spin and its NNs in each meron, the only variable parameter is the polar angle ϑ i . For merons centered at site 1, 9, and 17, the auxiliary angles are found to be meron A : and meron C : (sm-4) Here, elements in the second and third rows of Eqs. (sm-2)-(sm-4) are the polar angles ϑ i and azimuthal angles ϕ i , respectively. The other six spins rely on two angles (ϑ 0 , ϕ 0 ), subjecting to the following relation Hence, the Hamiltonian H could be recast as a multivariable function F(ϑ 0 , ϕ 0 ; ϑ 1 , ϑ 2 , ϑ 3 ) which depends on five variational parameters, and the ground-state energy is determined by minimizing F in the allowed parameter spaces. We emphasize that proper trial angles ({ϑ i }, ϕ 0 ) guided by the Monte Carlo simulation are essential to reach the global minima of F.
As shown in the Fig. sm-12, the classical KΓA model has a rich phase diagram which contains scores of magnetically ordered phases. Here, five of them are conventional magnetic orders which also frequently appear in other spin models. In the limit of ferromagnetic (FM) Kitaev point, there are two FM phases termed FM I phase and FM V phase. In the former the spins lie in the ab plane while in the latter the spins point out of the plane. Meanwhile, in the vicinity of antiferromagnetic (AFM) Γ limit, there are an AFM V phase for a negative A, and a zigzag (ZZ) phase and an in-plane 120 • I phase in the case of positive A. We note that the ZZ phase lies in the ac plane with a tilted angle α τ away from the a axis. The ground-state energy E g and the polarization directions in these phases can be given analytically, which are shown in Tab. sm-1.  In the KΓA model, when K < 0 and Γ > 0 the system is highly frustrated and many magnetic orders with a large unit cell could be stabilized. To illustrate it, we start by performing the Monte Carlo simulations along the line of φ = 0.68π and tune the polar angle θ from 0.25π (ZZ phase) to 0.7π (12 C phase). In between, we find six distinct phases, 18 B , TmX, 48 A , 6 B , 20 and 8 B , with the increase of θ. While closed forms of their ground-state energy are not available, the expression F that only relies on a few variational parameters could be obtained by the energy optimization method. The ground-state energy is thus targeted by minimizing F w.r.t. the variational parameters in the neighboring of optimal trial values. The energy of these phases is shown in Fig. sm-2, and the transition points are identified by the level-crossing points between neighboring phases. Since the magnetic unit cell of each phase is unknown aforehand, we need to determine their magnetic unit cell numerically. In Monte carlo simulations, this can be done by varying the cluster sizes and compare their ground-state energy. Here we exemplify this procedure by applying it to the TmX phase. Again we choose θ = 0.515π, φ = 0.68π as the representative point. If the cluster does not match the magnetic unit cell, we will end with a higher energy than that of the true ground state in the Monte Carlo simulations. In such a case, its ground state energy may approach the true one as the size increases. Therefore we need to adjust the cluster size and compare the energy to determine the magnetic unit cell, which is shown in Fig. sm-4 with primitive lattice vetors r 1 , r 2 given in Fig. sm-3. In Fig. sm-4(a), we fix L r1 = 6 and vary L r2 from 3 to 12. We find that when L r2 = 3n with n an integer we have the same lowest energy E g = −1.0065447701 while the energy is higher for L r2 = 3n + 1 and 3n + 2. Moreover, the difference decreases as n increases. This is well consistent with our expectation. In Fig. sm-4(b), we fix L r2 = 6 and vary L r1 from 3 to 12. We can find the same results as that from Fig. sm-4(a). In Fig. sm-4(c), we fix the ratio L r1 /L r2 = 1 and vary L r1 from 3 to 36. The lowest energy E g is at L r1 = 3n and it is independent of n. Moreover, it is the same as that in Fig. sm-4(a) and Fig. sm-4(b). These results obviously tell us the magnetic unit cell of the TmX phase has 2 × 3 × 3 = 18 spins.
We want to mention that to determine the unit cell of a magnetic ordered phase, this procedure is necessary. We also apply the same method to determine other phases in the phase diagram.

Sec. III: Some formal analysis on the TmX phase
Once we know the magnetic unit cell and the corresponding magnetic order, we can do further analysis. As demonstrated in Sec. I, the variational function F(ϑ 0 , ϕ 0 ; ϑ 1 , ϑ 2 , ϑ 3 ) of E g in the TmX phase can be written as The ground state energy E g (= min(F/18)) is readily available by numerically optimizing such a variational function. The equations ∂F/∂ξ = 0(ξ = ϑ 0 , ϕ 0 , ϑ 1 , ϑ 2 , ϑ 3 ) are safisfied at the minimum and they are used to check the results. In Fig. sm-5, we compare the ground-state energy obtained by the energy optimization method and the Monte Carlo method. We find that they are well consistent, demonstrating the validity of the function F(ϑ 0 , ϕ 0 ; ϑ 1 , ϑ 2 , ϑ 3 ). From the global phase diagram Fig. sm-12, one can see that ZZ, 18 B , TmX, 20, 8 B , AFM V and 120 • I phases intersect at the point Γ = 1 (i.e., A = 0, K = 0). The ground state of such pure Γ model is exactly known to be classical spin liquid [1] with its E g = −Γ. Moreover, E g = −Γ for the Γ model can also be obtained exactly from F(ϑ 0 , ϕ 0 ; ϑ 1 , ϑ 2 , ϑ 3 ). This further suggests that the TmX state is rooted in the Γ model and it is one of the highly degenerate ground states of the Γ model. To obtain a stable TmX state, a reasonable strategy is to add proper interactions to lift the degeneracy so that the TmX state becomes energetically favored. The Kitaev interaction can partly lift the degeneracy and the degenerate TmX and 18 B become the ground states. As we have shown in the main text, the degeneracy of the TmX and 18 B is further lifted by the single-ion anisotropy A. For a negative A, the TmX state is energetically favored. First, we discuss briefly how to calculate the topological charges Q of the merons in our lattice model. The definition was given by Berg and Lüscher [2], which is illustrated in Fig. sm-6. Actually, in a lattice the topological charge can be obtained by summing up the solid angle on every elementary triangle. For this purpose, the hexagon of each meron is splitted into six equilateral triangles, each is formed by a core spin and its two next-nearest neighbors. Among these six triangles, three of them contain an additional spin (the nearest neighbor of the core), which further divides the equilateral triangle into three isosceles triangles. The solid angle for each triangle is thus calculated following the definition given by Berg and Lüscher [2]. In our work we present the values of the topological charges up to six digits for simplicity although the numerical errors are of the order 10 −14 or smaller. In Fig. sm-7 we plot the topological charges Q of the three merons and their summation as a function of θ. Although the outmost spins might be driven off ab-plane, we can see that Q B and Q C are still close to ±1/2 and Q A is roughly 1. Importantly, their summation is always 1.0. We want to mention that in the TmX phase we can find parameters with the topological charges of Q A , Q B , Q C much closer to their theoretical values. For example, at θ = 0.515π, φ = 0.656π, Q A = 0.999881, Q B = 0.500058 and Q C = −0.499939, respectively. Moreover, due to the quadratic form of our model, a degenerate ground state is obtained by flipping all the spins. Under this transformation, a straightforward conclusion is that the topological charges of corresponding merons are sign-reversed.
Next, we discuss the details about calculating topological Hall effect (THE). Since there is no THE in the meronantimeron crystals [3], the THE caused by the TmX is highly nontrivial for the merons crystals. The double-exchange model used in main text means the coupling between itinerant electrons and the local magnetization is much larger than the hoping amplitude, so the spin of itinerant electron should be parallel to the direction of local magnetization S i . Hence, the energy of Kondo coupling term turns to be a constant. By defining a rotation matrix U i on each site which rotates the eigenvector of the itinerant electron from σ z to σ · S i , the above Hamiltonian can be transformed to a spinless free-fermion model [4]: The band structure E n (q) as well as wave functions φ n (q) can be directly obtained through the above effective Hamiltonian. The topological property of a band can be determined through its Chern number C n , which by definition is the integration of the Berry curvature Σ over the first Brillouin zone: with A(q) = −i φ n (q)|∇ q |φ n (q) the Berry connection. The Chern number of the first 9 bands presented in main text are 0, 0, 0, 1, 0, -2, 2, 1, 2. Since there is no band gap between 7-8 band and between the 8-9 bands, we focus on the 4-5 and 6-7 bands. Fig. sm-8 and Fig. sm-9 shows the full band structure and Berry curvature of selected bands. There is a clear Dirac cone structure with certain large gap in the vicinity of Γ point between 6 and 7 bands. As a results, one can see that the Berry curvature mainly concentrates around Γ point (Fig. sm-8). As a contrast, there are two Dirac cones at Γ and K 2 points between 4 and 5 bands resulting in very narrow global gap, and the Berry curvature dominates around Γ and K 2 (Fig. sm-9).
One may notice that the quantized Hall conductance shown in the main text is just the total Chern number when the chemical potential is inside the gap between the bands 6 and 7, i.e., σ xy = e 2 /h En<E f C n , since the Berry curvature could be equivalently obtained through Kubo formula as well. There is also a very narrow conductance platform between the bands 4 and 5, due to the tiny global gap there. Since the TmX state is found in the ground-state phase diagram, an important question to be answered is whether we can obtain such state experimentally. While Monte Carlo simulation is an elaborate and powerful metheod for equilibrium state problems, the dynamics of magnetic properties like: domain-wall motion, Hysteresis loop, piningdepinning problem, etc., are usually treated with Landau-Lifshitz-Gilbert equation [5,6]. To this end, we perform atomistic spin-dynamics simulations on field-cooling process of TmX state. Here, we choose again the representative point θ = 0.515π, φ = 0.68π. The Landau-Lifshitz-Gilbert (LLG) equation is numerically solved in the xyz coordinate system: the effective field felt by S i , α a dimensionless damping factor and h c ext an external magnetic field along c direction. For finite temperatures, a stochastic fluctuation field h f l is then added to the effective field to act as the thermal noise, which should have a zero average value and be uncorrelated in spin component, space and time: with 2 = 2α k B T and T the temperature determining the strength of thermal fluctuations. In our simulation, we set k B = = 1 (a dimensionless system), α = 0.01 (different values of α have also been tested) and solve the LLG equation using the fourth-order Runge-Kutta method on a 288 × 288 × 2 honeycomb lattice. 6 × 10 6 iterations are performed between every two successive temperatures. As we have already noticed, the TmX state is degenerate but with an opposite Q under the transformation S → −S. These two states might appear simultaneously and annihilate each other when starting from a random initial state. Such a degeneracy can be easily lifted by applying an external magnetic field h c ext along the c-axis. As shown in Fig. sm-10(a)-(c), one can already find the TmX state when the temperature T = 0.01. As T decreases, the TmX state becomes more clear, suggesting the field cooling is indeed effective to obtain a stable TmX state (Supplementary Movie Sabc). However, one may notice that there are two kinds of TmX states with their core spins on different sublattices, say, a and b, which are separated by narrow domains (Fig. sm-10(c) and Fig. sm-10(d)). Such a phenomenon arises from the bipartite nature of the honeycomb lattice. To eliminate such domains, it is necessary to lift the degeneracy further. In our numerical simulations, a simple strategy is to change A on one sublattice, say a, to αA with α equal to, for example, 0.9. We want to mention that in experiments one can break the symmetry of the lattice by applying strain. The corresponding results are shown in Fig. sm-10 can find a nearly perfect TmX state at T = 0.001 (Supplementary Movie Sefg). We want to stress that lifting such degeneracy is definitely helpful to eliminate domains but the effect depends on α as well as the initial configuration. As shown in Fig. sm-11, the parameters and simulation details in Fig. sm-11 and Fig. sm-10 are the same except that the random initial configuration is different. However, in Fig. sm-11(g) there are still domains. Even though, the effect of lifting the degeneracy can be demonstrated by counting the number of merons with their cores on the different sublattice, say N a or N b . In Fig. sm-11(c) where T = 0.001, N a : N b = 7524 : 7328, which is very close to 1 : 1. This ratio again evidences the degeneracy of the TmX state. As we change the coefficient of the single-ion anisotropy A on the sublattice a to 0.9A in Fig. sm-11(e)-(h), at T = 0.001 ( Fig. sm-11(g)), the region of the energetically favored Tmx state is obviously increased, i.e., N a : N b = 12337 : 3937. One may naively expect that core spins on the b sublattice is energetically favored because α < 1.0 and A < 0. Actually, this is not true. Changing α simultaneously changes the spin directions of NNs and NNNs and one can not just compare the energy of the core spins. We confirm the core spins prefer staying on the a sublattice by the Monte Carlo simulations.
Sec. VI: Complete phase diagram, Spin configurations and static spin structure factors In this section, we show the complete phase diagram, the spin configurations and the corresponding static spin structure factors F ξξ q = 1 N ij e iq·(R i −R j ) S ξ i S ξ j (ξ = a, b, c). The complete phase diagram of the KΓA model is shown in Fig. sm-12 as a semicircle plot. The radial distance represents θ and the polar angle represents φ. In the region π/2 < φ < π, the model is highly frustrated and a large number of phases are found. The number in the phase diagram marks the number of spins in one magnetic unit cell and the subscripts A, B, C and D are used to distinguish different phases marked by the same number. The phase transitions of all the ordered phases are of the first order except that the transitions to the wave phase is unclear. Note that on the line θ = π, where K = 0 and Γ = 0, the Hamiltonian is decoupled and then become trivial. From the symmetry, it is ready to know that the line φ = π, or equivalently φ = 0, is the phase transition line.
In the following, we will show the spin configuration and the static spin structure factors of each phase. We will not show the five phases discussed in Sec. I since they are common and their properties are well known in magnets. For simplicity, we choose one representive point in each phase to illustrate the spin structure, which are shown in Fig. sm-13-Fig. sm-26 with the phase name and corresponding parameters given in the caption. In each figure, the spin configurations are shown on the left panels and the corresponding static spin structure factors are on the right panels. On the right panel, the a, b or c in the subfigures marks the spin component and the right bottom subfigures show the summation of all three components. Since there two regions belonging to the 6 A phase, we choose one point in each region and plot them in Fig. sm-15 and Fig. sm-16. The magnetic unit cell of the wave phase is rather complex. It may vary as a function of θ and φ. At present we can not draw a reliable conclusion. The two figures sm-35 and sm-36 are plotted just for reference.