Topological spin textures in chiral magnets on the honeycomb lattice with magnetic fields

Li-Peng Jin 11930771@mail.sustech.edu.cn Department of Physics, Southern University of Science and Technology, Shenzhen 518005, China Bin Xi xibin@yzu.edu.cn College of Physics Science and Technology, Yangzhou University, Yangzhou 225002, China Yong-Jun Liu College of Physics Science and Technology, Yangzhou University, Yangzhou 225002, China

Abstract

Topological spin textures, like skyrmions, have significant potential for spintronics applications. The main purpose of this work is to study further the topological spin textures on the discrete lattice with magnetic fields. In this work, we study a classical rotated Heisenberg model with Dzyaloshinskii-Moriya interaction, bond-dependent anisotropy and easy-axis anisotropy on the honeycomb lattice via Monte Carlo simulations. We mainly focus on phase diagrams with magnetic fields, especially on the non-trivial phases only existing with fields. The results demonstrate the emergence of field-induced incommensurate skyrmion superlattice including mixed skyrmion-bimeron states, ferromagnetic star, vortex and z-vortex superlattice. We systematically analyze structures of these topological spin textures through spin configurations, structure factors, topological charge and vector chirality. We hope that our results could be useful for the identification of topological phases experimentally.

Keywords: skyrmions, chiral magnets, spin-orbit coupling, Monte Carlo

pacs:

12.39.Dc, 75.70.Kw, 71.70.Ej, 52.65.Pp

I Introduction

Magnetic skyrmions [1, 2] are representative topological spin textures, which have been observed or proposed due to four main mechanisms, such as magnetic dipolar interaction[3, 4], Dzyaloshinskii-Moriya (DM) interaction[5, 6, 7, 8], frustrated exchange interaction[9] and four-spin exchange interactions[10]. The DM interaction arises from the relativistic spin-orbit coupling (SOC), which has been extensively studied in continuous models, such as in chiral magnets[11, 12, 13, 14, 15]. In addition, due to the experimental realization of synthetic SOC for ultra-cold atoms[16], more attention has been paid to discrete lattice models, such as the extended Hubbard model on square lattice[17]. In the strong coupling limit, this model could derive the rotated Heisenberg model with DM interaction and bond-dependent anisotropy[19, 20, 18, 23, 22, 24, 21], which support magnetic textures such as spin spirals and vortex and skyrmion superlattice.

However, the magnetic fields also play an important role in the formation of topological spin textures[25, 26]. Motivated by these, we study further the rotated Heisenberg model on the honeycomb lattice via Monte Carlo simulations. We mainly focus on phase diagrams with out-of-plane magnetic fields, especially on the non-trivial phases only existing with fields. The results of zero fields are consistent with the previous work[17, 27]. On the contrast, the results with fields indicate the emergence of ferromagnetic star (FM star)[28], incommensurate skyrmion superlattice (IC-SkL) including mixed skyrmion-bimeron states[29], and z-vortex superlattice (Z-VL). We demonstrate a distinct type of topological structure, bimeron that consists of two merons (half-skyrmions) separated by a spiral domain with zero topological charge, which could further turn into skyrmion superlattice upon increasing the magnetic fields. The FM star contains a large hexagon with the in-plane spin component rotating $2\pi$ in the core hexagon and six spins in the outer hexagon aligning in a parallel fashion. We characterize all these topological spin textures and give the structure factors for the identification of topological phases experimentally.

The paper is organized as follows. In Sec. II, we start from the bosonic model on the honeycomb lattice with a synthetic SOC. In the strong coupling limit, we obtain the rotated Heisenberg model with DM interaction, bond-dependent anisotropy and easy-axis anisotropy. We give details about the Monte Carlo simulation and how to characterize different spin configurations. In Sec. III.1, we map out phase diagrams with/without fields and mainly focus on phase diagrams with fields, especially on the phases only existing with fields. In Sec. III.2, we describe the structures of topological spin textures and give the structure factors for the identification of these topological phases experimentally. In Sec. IV, we give our discussion and conclusion.

II Model and Method

The low-energy behavior of bosons on the honeycomb lattice with SOC could be described with the following Hamiltonian:

$\displaystyle\hat{\mathcal{H}}_{\textrm{boson}}$	$\displaystyle=$	$\displaystyle-t\sum_{\langle ij\rangle}\left(\hat{c}^{\dagger}_{i\tau}\mathcal{R}_{ij}\hat{c}_{j\tau^{\prime}}+\mathrm{H.c.}\right)+\hat{\mathcal{H}}_{U},$
$\displaystyle\mathcal{R}_{ij}$	$\displaystyle=$	$\displaystyle\exp\left[i\bm{A}\cdot\left(\bm{r}_{i}-\bm{r}_{j}\right)\right],$
$\displaystyle\hat{\mathcal{H}}_{U}$	$\displaystyle=$	$\displaystyle\frac{U}{2}\sum_{i\tau}\hat{n}_{i\tau}(\hat{n}_{i\tau}-1)+U^{\prime}\sum_{i}\hat{n}_{i\uparrow}\hat{n}_{i\downarrow},$	(1)

where the first term in $\hat{\mathcal{H}}_{\textrm{boson}}$ describes the hopping term between nearest-neighbor sites, with $t$ representing the hopping amplitude and $\hat{c}^{\dagger}_{i\tau}$ ( $\hat{c}_{i\tau}$ ) is the creation (annihilation) operator at site $i$ , with subscripts $\uparrow$ and $\downarrow$ representing two internal hyperfine states of bosons. The $\bm{A}=(\theta\sigma_{y},-\theta\sigma_{x},0)$ in matrix $\mathcal{R}_{ij}$ is a non-Abelian gauge field by the bosons, which is the well-known Rashba term. The second term $\hat{\mathcal{H}}_{U}$ describes on-site interactions between bosons, where $U$ is intracomponent interaction and $U^{\prime}$ is the intercomponent one, $\hat{n}_{i\tau}=\hat{c}^{\dagger}_{i\tau}\hat{c}_{i\tau}$ is the boson number operator with spin $\tau$ at site $i$ .

At unit filling and in the strong coupling limit $U/t\gg 1$ , the bosonic model can be approximated as a rotated Heisenberg model, the effective low-energy spin Hamiltonian given by:

$\displaystyle\hat{\mathcal{H}}_{\text{spin}}$	$\displaystyle=$	$\displaystyle J\sum_{\langle i,j\rangle}\left\{J_{0}\bm{S}_{i}\cdot\bm{S}_{j}+\bm{D}_{ij}\cdot(\bm{S}_{i}\times\bm{S}_{j})\right.$
		$\displaystyle+A_{c}\left[\bm{S}_{i}\cdot\left(\bm{\hat{e}}_{z}\times\bm{r}_{ij}\right)\right]\left[\bm{S}_{j}\cdot\left(\bm{\hat{e}}_{z}\times\bm{r}_{ij}\right)\right]\}$
		$\displaystyle+J\sum_{i}\Delta({S_{i}}^{z})^{2},$
$\displaystyle J=$	$\displaystyle-$	$\displaystyle 4t^{2}/U,\quad J_{0}=\cos 2\theta,$
$\displaystyle\bm{D}_{ij}$	$\displaystyle=$	$\displaystyle\sin 2\theta(\bm{\hat{e}}_{z}\times\bm{r}_{ij}),A_{c}=2\sin^{2}\theta,$	(2)

we set $|J|=1$ as the energy unit scale. The first term is the nearest-neighbor Heisenberg interaction. The second term is DM interaction, which is anisotropic antisymmetric interaction. The third term is bond-dependent anisotropy, which is symmetry-allowed even in the presence of inversion symmetry. The last term in $\hat{\mathcal{H}}_{\text{spin}}$ is easy-axis anisotropy $(\Delta>0)$ , which is derived from intercomponent interaction $U^{\prime}$ between bosons.

We study the honeycomb lattice with $L\times L$ ( $L=60$ ) unit cells placed in the $xy$ -plane, where each unit cell contains two sites colored with blue and red points shown in Fig. 1, hence the system contains $2\times L\times L=7200$ sites. We set the lattice constant $a=1$ and the magnitude of the spins is fixed by the normalization condition $|\bm{S}_{i}|=1$ . The periodic boundary conditions are imposed. Due to the anisotropy of the honeycomb lattice, there are three types of $\bm{r}_{ij}$ pointing from blue site to red site: $\bm{r}_{ij}^{\alpha}=(1,0)$ , $\bm{r}_{ij}^{\beta}=(-1/2,\sqrt{3}/2)$ and $\bm{r}_{ij}^{\gamma}=(-1/2,-\sqrt{3}/2)$ . The three types of DM vectors ( $\bm{D}_{\alpha},\bm{D}_{\beta},\bm{D}_{\gamma}$ ) are shown with green arrows. We use the convention $\bm{D}_{ij}\cdot(\bm{S}_{i}\times\bm{S}_{j})$ in which the first spin in the cross product $\bm{S}_{i}$ is always on the sublattice color with blue. The lattice is constructed from different sets of primitive vectors $(\bm{a}_{1},\bm{a}_{2})$ shown with yellow arrows.

Refer to caption — Figure 1: Sketch of the honeycomb lattice, with the primitive vectors $\bm{a}_{1},\bm{a}_{2}$ , and $\bm{D}_{\alpha},\bm{D}_{\beta},\bm{D}_{\gamma}$ corresponding to the DM vectors in $\alpha$ , $\beta$ and $\gamma$ bonds, respectively.

To reveal the true classical ground states of $\hat{\mathcal{H}}_{\text{spin}}$ in Eq. (2), we use the paralleling-annealing Monte Carlo (MC) simulation[30] on 40 replicas, with temperature $T/|J|$ ranging from $0.001$ to $1.0$ . For each replica, we sample it with a combination of heat-bath and over-relaxation methods[31]. A whole MC step consists of a single heat-bath sweep and subsequent ten over-relaxation sweeps over the entire lattice. We perform $2\times 10^{5}$ MC steps per replica, then we copy out the spin configuration from the lowest-T replica and sample it with a combination of zero-temperature heat-bath and over-relaxation method to obtain the ground states. The zero-temperature heat-bath sampling is simply aligning the spins according to their local fields:

\displaystyle\bm{S}_{i}

\displaystyle=

\displaystyle\frac{\bm{h}_{i}^{loc}}{|\bm{h}_{i}^{loc}|}S,

(3)

with $\bm{h}_{i}^{loc}={\partial\hat{\mathcal{H}}_{\text{spin},i}}/{\partial\bm{S}_{i}}$ .

The important observables that can be used for the identification of magnetic phases are the spin structure factors given by:

\displaystyle\bm{F_{k}}

\displaystyle=

\displaystyle\frac{1}{N}\sum_{ij}\left\langle\bm{S}_{i}\cdot\bm{S}_{j}\right\rangle e^{i\bm{k}\cdot(\bm{r}_{i}-\bm{r}_{j})},

(4)

where $\bm{k}$ is the reciprocal space vector, $\bm{r}_{i}$ is the lattice vector at site $i$ and $N$ the number of sites.

Each skyrmion contributes $\pm 1$ to the topological charge (skyrmion number). To calculate the topological charge, we employ the definition on the lattice version introduced by Berg[32, 33]. The calculation of $Q$ starts by triangulating the entire lattice and then counting the solid angle $\Omega_{l}$ for each triangle $l$ through all three spins $(\bm{S}_{1},\bm{S}_{2},\bm{S}_{3})$ located at its vertices, and the topological charge can be obtained as $Q=\sum_{l}\Omega_{l}/4\pi$

We also calculate the vector chirality $\bm{\kappa}$ for each triangle:

\displaystyle\bm{\kappa}=\frac{2}{3\sqrt{3}}\left(\bm{S}_{1}\times\bm{S}_{2}+\bm{S}_{2}\times\bm{S}_{3}+\bm{S}_{3}\times\bm{S}_{1}\right)

(5)

where three spins $\bm{S}_{1}$ , $\bm{S}_{2}$ and $\bm{S}_{3}$ are chosen anti-clockwise on the triangle.

III Results

III.1 Phase diagrams

We obtain the classical ground-state phase diagram of $\hat{\mathcal{H}}_{\text{spin}}$ in Eq. (2) via the MC simulation, with tuning the parameter $\theta$ and the easy-axis anisotropy $\Delta$ . We map out phase diagrams with/without and get the following phases characterized by the spin structure factors $\bm{F_{k}}$ . We mainly focus on phase diagram with magnetic fields, especially on the non-trivial phases only existing with fields. To identify the different topological phases realized in the phase diagrams, we used the calculated topological charge and vector chirality to visualize magnetic configurations.

III.1.1 phase diagram without magnetic fields

In Fig. 2(a), we plot the phase diagram without magnetic fields, which is consistent with the previous work on the square lattice, including ferromagnetic phase (Z-FM), antiferromagnetic phase (Z-AFM), spiral phases and two topological phases ( $3\times 3$ SkL and $\sqrt{3}\times\sqrt{3}$ VL). The $3\times 3$ SkL and $\sqrt{3}\times\sqrt{3}$ VL phases mean the $3\times 3$ unit cell skyrmion superlattice and $\sqrt{3}\times\sqrt{3}$ unit cell vortex superlattice, respectively. The order parameters in FM, AFM, SkL, VL and spirals are given by average magnetic moment $m$ $(m=N^{-1}\sum_{i}{S_{i}}^{z})$ , topological charge $|Q|$ and vector chirality $\bm{\kappa}$ , respectively. In Fig. 2(b)-(d), we plot these order parameters as a function of $\theta$ at $\Delta=0.2$ .

III.1.2 phase diagram with magnetic fields

With applied out-of-plane external magnetic fields in the Hamiltonian of $\hat{\mathcal{H}}_{\text{spin}}$ , we obtain two different phase diagrams corresponding to different easy-axis anisotropy ( $\Delta=0.2$ and $0.5$ ), marked with dashed lines in Fig. 2(a). In Fig. 3(a) and Fig. 3(e), we map out the two phase diagrams with out-of-plane external magnetic fields. Most of the phases in the phase diagrams of ground states are still stable when applying external magnetic fields. It is wondrous that three field-induced phases have emerged, viz IC-SkL phase, $2\times 2$ FM star phase and $\sqrt{3}\times\sqrt{3}$ Z-VL phase, respectively. The order parameters in these phases are also given by average magnetic moment $m$ , topological charge and chirality, respectively. In Fig. 3, we plot these order parameters as a function of $\theta$ with different parameters. Significantly, the topological charge $|Q|$ shown in Fig. 3(c) and Fig. 3(g) indicate that the number of skyrmions in IC-SkL phase varies with different parameters.

We briefly talk about the non-trivial phases in phase diagrams including ferromagnetic phase, antiferromagnetic phase and spiral phases. Z-FM and Z-AFM phases result from the competition between ferromagnetic or antiferromagnetic exchange interaction and easy-axis anisotropy. Z-FM phase is the ferromagnetic phase where spins orient along the $\pm z$ axis. The spin structure factor $\bm{F_{k}}$ exhibits a peak at $\bm{k}=(0,0)$ , which is the $\bm{\Gamma}$ point in the first Brillouin zone. Analogously, the Z-AFM phase is the antiferromagnetic phase where spins point along with the $\pm z$ axis. The $\bm{F_{k}}$ exhibits peaks at $(\pm 2\pi/3,\pm 2\sqrt{3}\pi/3)$ and $(\pm 4\pi/3,0)$ .

Spiral phases occupy a large region in phase diagram without fields, which mostly result from the DM interaction. Spiral phases are coplanar states, where spins are all parallel at a plane and their direction rotates by a constant angle from one spin to a neighboring spin along the helical axis. When $\theta$ and $\Delta$ are small, the ferromagnetic exchange interaction is much more than the DM interaction, which leads to IC-Spiral phases in zero fields. IC-Spiral phases are incommensurate spiral states. As $\theta$ increases, the DM interaction increasing, the emergence of $N^{*}\times 1$ spiral phases in phase diagram. $N^{*}\times 1$ spiral phases are different sizes of spiral states with $2N^{*}$ sites unit cell, where $N^{*}=2,3,4,5,6$ and $N^{*}$ decreased with $\theta$ increasing. Spins spiral in the $z-\bm{k}$ plane and the $N^{*}\times 1$ spiral phases are triple degeneracy because of the rotation symmetry of lattice. However, with the magnetic fields applied, the spiral phases could turn into IC-SkL phases rapidly, only the $3\times 1$ and $4\times 1$ spiral phases still stable with magnetic fields. In particular, the $2\times 1$ spiral phase turn into the $2\times 2$ FM star phase and the $\sqrt{3}\times\sqrt{3}$ vortex superlattice phase turn into the $\sqrt{3}\times\sqrt{3}$ z-vortex superlattice phase under magnetic fields, which preserve the rotation symmetry.

Besides FM, AFM and spiral phases, we find rich varieties of topological spin textures in phase diagrams, such as SkL, VL, IC-SkL, FM star and Z-VL phases, which have significant potential for spintronics applications. Remarkably, these topological spin textures remain a large region in phase diagrams even with magnetic fields, which is important for spintronics applications. We describe the structures of these topological spin textures in details in the following and especially focus on these non-trivial phases only existing with magnetic fields.

III.2 Structures of topological spin textures

We describe both the real-space spin configurations and structure factors of all topological spin textures, where the structure factor is an important observation that can be used for the identification of magnetic phases experimentally. Most of the phases in the phase diagram without fields are still stable when applying external magnetic fields. Besides SkL and VL topological phases, three field-induced topological phases have emerged, viz IC-SkL phase, FM star phase and Z-VL phase. We investigate all these topological spin textures in details as follows.

III.2.1 $3\times 3$ skyrmion superlattice and Incommensurate skyrmion superlattice

The extremely strong competition between DM interaction and $A_{c}$ term leads to the $3\times 3$ SkL phase. The SkL phase is a non-coplanar state, where the spins form an 18-site skyrmion structure. In Fig. 4(a), we give the real-space spin configuration of the SkL phase and the dashed line denotes the 18-site skyrmion magnetic unit cell. The elementary skyrmion contains a large hexagon with $3\times 3$ unit cell of the honeycomb lattice. In the hexagon of $3\times 3$ unit cell, a spin with negative $S^{z}$ magnetization at the core, spins from the core lie in the $xy$ -plane rotating $2\pi$ and six outermost spins with positive $S^{z}$ magnetization. The structure factor of SkL phase has peaks at $\bm{k}=\frac{2}{3}\mathbf{M}$ points shown in Fig. 4(b). We calculated the topological charge for the $3\times 3$ SkL phase by the definition of topological charge shown in Fig. 2(c). The elementary skyrmion was marked with the dashed line in Fig. 4(a), where the topological charge $|Q|=1$ for each 18-site skyrmion.

There are similar IC-SKL phases exist in the phase diagrams with magnetic fields. The IC-Spiral states and larger size spiral states are not stable when applying the magnetic fields. With the magnetic fields increasing, the IC-Spiral states and large-size spiral states are superseded by IC-SkL phases rapidly. The IC-SkL phases are complicated states including the different sizes of skyrmion superlattice and mixed skyrmion-bimeron states, which result from the extremely strong competition among the ferromagnetic exchange interaction, DM interaction and the external magnetic fields.

IC-SkL phases are non-coplanar topological states, including non-coplanar skyrmion states and mixed skyrmion-bimeron states. In Fig. 5(a)-(d), we give four spin configurations of IC-SkL states corresponding to different parameters. In Fig. 5(a)-(b), we give the spin configurations of IC-SkL states by tuning the parameter $\theta$ . From Fig. 5(a) to Fig. 5(b), we fixed the value of $\Delta=0.2$ and $H/|J|=0.2$ , then set up the $\theta$ equal to $0.10\pi$ and $0.13\pi$ . The increasing of the DM interaction with tuning $\theta$ leads to the size of skyrmions decreasing and the number of skyrmions increasing, corresponding to the order parameter topological charge in Fig. 3(c). As shown in Fig. 6(a)-(b), the region occupied by the peaks in the structure factors also indicate the size of skyrmions decreasing.

In Fig. 5(c)-(d), we give the spin configurations and corresponding topological charge density of IC-SkL states by tuning the external magnetic fields $H/|J|$ . From Fig. 5(c) to Fig. 5(d), we fixed the value of $\Delta=0.5$ and $\theta=0.15\pi$ , then set up the $H/|J|$ equal to $0.2$ and $0.5$ . As the increasing of the magnetic fields with parameter $H/|J|$ , bimerons states contracting and the mixed skyrmion-bimeron states are superseded by the skyrmion superlattice gradually. In the spin configuration of mixed states (Fig. 5(c)), we could find skyrmions and bimerons, which consist of two merons carrying the topological charge $Q=-1/2$ and connected by a spiral domain with zero topological charge and are relatively stable in low fields. In Fig. 6(c)-(d), we give the structure factors of mixed skyrmion-bimeron states and skyrmion lattice, where the structure factor of mixed states has spread peaks because of the existence of bimerons.

III.2.2 $\sqrt{3}\times\sqrt{3}$ vortex superlattice and $\sqrt{3}\times\sqrt{3}$ z-vortex superlattice

At $\theta$ is close to $0.5\pi$ , the DM interaction is almost to zero and the $A_{c}$ term reaches its maximum. When $\Delta$ is close to $0$ , the spins coplanar in $xy$ plane because of the $A_{c}$ term, which lead to another topological spin texture, $\sqrt{3}\times\sqrt{3}$ VL phase. The VL phase is a coplanar state and spins in $xy$ -plane form a $\sqrt{3}\times\sqrt{3}$ unit cell vortex, where each vortex contains 6 sites. In Fig. 7(a), We give the real-space spin configuration of the VL phase and the dashed line denotes the 6-site vortex magnetic unit cell. The elementary vortex contains a hexagon with a 6-site vortex and in each vortex, the spins in $xy$ -plane wind anti-clockwise $2\pi$ around each hexagon. The structure factor of VL phase has peaks at $(\pm 2\pi/3,\pm 2\sqrt{3}\pi/9)$ and $(0,\pm 4\sqrt{3}\pi/9)$ shown in Fig. 7(c).

We get the $\sqrt{3}\times\sqrt{3}$ Z-VL phase in the region between VL phase and Z-AFM phase with magnetic fields, which mostly result from the competition between the $A_{c}$ term and the external out-of-plane magnetic fields. The $A_{c}$ term is favorable to form the in-plane vortex superlattice state. Whereas, when applying the magnetic fields, the out-of-plane components of spins show ferromagnetic order. The Z-VL phase is also a non-coplanar state, where the in-plane components of spins form a $\sqrt{3}\times\sqrt{3}$ unit cell vortex. In Fig. 7(b), We give the real-space spin configuration of the Z-VL phase and the dashed line denotes a 6-site magnetic unit cell. However, the out-of-plane components of spins show ferromagnetic order. The structure factor of Z-VL phase is similar to VL phase. There has an extra peak at $(0,0)$ , which is $\bm{\Gamma}$ point, corresponding to the ferromagnetic order of out-of-plane components shown in Fig. 7(d).

III.2.3 $2\times 2$ FM star

We get the $2\times 2$ FM star state in the region of $2\times 1$ spiral phase with applied external magnetic fields. The FM star state is another non-coplanar state, where the spins form an 8-site magnetic unit cell with vortex-like structure. The magnetic fields lead to the transformation from the $2\times 1$ spiral phase to the $2\times 2$ FM star. As in zero fields, the $2\times 1$ spiral phase breaks the rotation symmetry, whereas the FM star phase preserves the rotation symmetry. The FM star phase is favored by the competition between $A_{c}$ term, DM interaction, and magnetic fields. In Fig. 8(a), we give the real-space spin configuration of the FM star phase and the dashed line denotes an 8-site magnetic unit cell. The whole FM star contains a large hexagon with $2\times 2$ unit cell of the honeycomb lattice. In the hexagon of an 8-site magnetic unit cell, the in-plane spin component rotating $2\pi$ in the core hexagon and six spins in the vertices of the outer hexagon are aligned in a parallel fashion.

To characterize and identify the $2\times 2$ FM star phase, we also calculate the structure factor. In Fig. 8(b), we give the structure factor of FM star phase, which has peaks at $(\pm\pi/3,\pm\sqrt{3}\pi/3)$ , $(\pm 2\pi/3,0)$ and an extra peak at $(0,0)$ . These peaks are $\mathbf{M}$ points of the first Brillouin zone and the peak at $\bm{\Gamma}$ point corresponding to the ferromagnetic order of six outermost spins. The $2\times 2$ FM star phase remains stable in a large parameter region with magnetic fields, which could be expected to be found in future experimental work.

IV Discussion and Conclusion

In this paper, we find various topological spin textures. Remarkably, these topological spin textures remain a large region in phase diagrams even with magnetic fields, which have significant potential for spintronics applications and could further understand the effect of magnetic fields on topological spin textures. These novel topological spin textures could be reproduced in both real chiral magnets and ultra-cold atoms systems. We hope that our results could be achieved in future experimental work.

In conclusion, we study the rotated Heisenberg model on the honeycomb lattice via Monte Carlo simulations, which could describe the low-energy behavior of a bosonic model with SOC in the deep Mott insulating region or two-dimensional chiral magnets with strong SOC. We obtain the classical phase diagrams with/without fields and mainly focus on phase diagrams with fields, especially on the non-trivial phases only existing with fields. We found $N^{*}\times 1$ spiral phases, IC-Spiral phase, SkL phase and VL phase in zero fields. As the magnetic fields were applied, we found IC-SkL phases, where the number of skyrmions changes with $\theta$ and mixed skyrmion-bimeron states appear at $\Delta=0.5$ , $\theta=0.15\pi$ , $H/|J|=0.2$ , FM star phase and Z-VL phase. We describe both spin configurations and spin structure factors of topological spin textures in details, where the structure factor of mixed skyrmion-bimeron states has spread peaks compared with SkL states. We also found the structure factor of FM star phase has peaks at $\mathbf{M}$ points and the structure factor of VL phase has peaks at $\mathbf{K}$ points. Overall, our results of these topological spin textures could have significant potential for the identification of topological phases experimentally.

V Acknowledgments

We thank Ji-Ze Zhao for the useful discussion. This work was supported by the National Natural Science Foundation of China (Grants No. 11774300).

References

[1] Skyrme T H R 1962 Nucl. Phys. 31 556.
[2] Bogdanov A and Huber A 1994 J. Magn. Magn. Mater. 138 255.
[3] Lin Y S, Grundy P J, and Giess E A 1973 Appl. Phys. Lett. 23 485.
[4] Yu X Z et al 2012 Proc. Natl Acad. Sci. 109 8856.
[5] Dzyaloshinskii I 1958 J. Phys. Chem. Solids. 4 241.
[6] Moriya T 1960 Phys. Rev. 120 91.
[7] Mühlbauer S et al 2009 Science 323 915.
[8] Yu X Z et al 2010 Nature 465 901.
[9] Okubo T, Chung S, and Kawamura H 2012 Phys. Rev. Lett. 108 017206.
[10] Heinze S et al 2011 Nat. Phys. 7 713.
[11] Nagaosa N and Tokura Y 2013 Nat. Nanotechnol. 8 899.
[12] Iwasaki J, Mochizuki M and Nagaosa N 2013 Nat. Commun. 4 1463.
[13] Lin S Z, Batista C D, Reichhardt C and Saxena A 2014 Phys. Rev.Lett. 112 187203.
[14] Tokunaga Y et al 2015 Nat. Commun. 6 7638.
[15] Fert A, Reyren N and Cros V 2017 Nat. Rev. Mater. 2 17031.
[16] Lin Y J, Jiménez-García K and Spielman I B 2011 Nature (London) 471 83.
[17] Cole W S, Zhang S Z, Paramekanti A and Trivedi N 2012 Phys. Rev. Lett. 109 085302.
[18] Banerjee S, Rowland J, Erten O and Randeria M 2014 Phys. Rev. X 4 031045.
[19] Shekhtman L, Entin-Wohlman O and Aharony A 1992 Phys. Rev. Lett. 69 836.
[20] Yildirim T, Harris A B, Entin-Wohlman O and Aharony A 1994 Phys. Rev. Lett. 73 2919.
[21] Cai Z, Zhou X F and Wu C J 2012 Phys. Rev. A 85 061605.
[22] Zhao J Z, Hu S J and Zhang P 2015 Phys. Rev. Lett. 115 195302.
[23] Sun F D, Ye J W and Liu W M 2015 Phys. Rev. A 92 043609.
[24] Xi B, Hu S J, Luo Q, Zhao J Z and Wang X Q 2017 Phys. Rev. B 95 014405.
[25] Yu X Z et al 2018 Nature 564 95.
[26] Wang Z T, Su Y, Lin S Z and Batista C D 2021 Phys. Rev. B 103 104408.
[27] Sun F D and Ye J W 2016 arXiv:1603.00451 [cond-mat.quant-gas].
[28] Janssen L, Andrade E C and Vojta M 2016 Phys. Rev. Lett. 117 277202.
[29] Iakovlev I A, Sotnikov O M, and Mazurenko V V 2018 Phys. Rev. B 97 184415.
[30] Hukushima K, and Nemoto K 1996 J. Phys. Soc. Jpn. 65 1604.
[31] Miyatake Y et al 1986 J. Phys. C: Solid State Phys. 19 2539.
[32] Berg B and Lüscher M 1981 Nucl. Phys. B 190 412.
[33] Yin G et al 2016 Phys. Rev. B 93 174403.