Magnetic Skyrmion State in Janus Monolayers of Chromium Trihalides Cr(I,X)3

TL;DR

This paper proposes creating Janus monolayers of chromium trihalides to induce Dzyaloshinskii-Moriya interactions, enabling the formation of magnetic skyrmions in two-dimensional materials, which are promising for future nano-scale devices.

Contribution

It introduces a novel fabrication approach for Janus monolayers of Cr(I,X)3 that exhibit skyrmionic states due to engineered DMI, advancing 2D magnetic skyrmion research.

Findings

Cr(I,Br)3 hosts intrinsic skyrmionic states.

Cr(I,Cl)3 exhibits magnetic-field-induced skyrmions.

Janus monolayers enable DMI and skyrmion formation in 2D magnets.

Abstract

Magnetic skyrmions are nano-scale spin structures that are promising for ultra-dense memory and logic devices. Recent progresses in two-dimensional magnets encourage the idea to realize skyrmionic states in freestanding monolayers. However, monolayers such as CrI3 lack Dzyaloshinskii-Moriya interactions (DMI) and thus do not naturally exhibit skyrmions but rather a ferromagnetic state. Here we propose the fabrication of Cr(I,X)3 Janus monolayers, in which the Cr atoms are covalently bonded to the underlying I ions and top-layer Br or Cl atoms. By performing first-principles calculations and Monte-Carlo simulations, we identify strong enough DMI, which leads to not only helical cycloid phases, but also to intrinsic skyrmionic states in Cr(I,Br)3 and magnetic-field-induced skyrmions in Cr(I,Cl)3.

Tables

Table 1. Table 1: Magnetic parameters of CrI 3 , Cr(I,Br) 3 and Cr(I,Cl) 3 monolayers. Note that (i) the J = ( J x ​ x + J y ​ y + J z ​ z ) / 3 𝐽 subscript 𝐽 𝑥 𝑥 subscript 𝐽 𝑦 𝑦 subscript 𝐽 𝑧 𝑧 3 J=(J_{xx}+J_{yy}+J_{zz})/3 coefficient listed below is different from the J ′ = ( J α + J β ) / 2 superscript 𝐽 ′ subscript 𝐽 𝛼 subscript 𝐽 𝛽 2 J^{\prime}=(J_{\alpha}+J_{\beta})/2 parameter of Eq. (1), with the former being more reasonable when comparing strength with D 𝐷 D ; and (ii) S = 3 / 2 𝑆 3 2 S=3/2 is used when extracting the magnetic parameters. The energy unit is in meV.

Mag. Para.	CrI3	Cr(I,Br)3	Cr(I,Cl)3
$J_1$	-2.218	-1.800	-0.983
$K_1$	0.847	0.505	0.422
$D_1$	0	0.270	0.191
$\left|D_1/J_1\right|$	0	0.150	0.194
$J_2$	-0.638	-0.673	-0.754
$K_2$	0.078	0.075	0.094
$D_2$	0.135	0.114	0.138
$\left|D_2/J_2\right|$	0.212	0.169	0.183
$A_{zz}$	-0.262	-0.124	-0.029

Equations

$$\mathcal{H}=\frac{1}{2}\sum _{n=1,2}\sum _{(i,j{)}_n}{\mathrm{S}}_i\cdot {\mathcal{J}}_{n,ij}\cdot {\mathrm{S}}_j+\sum _i{\mathrm{S}}_i\cdot {\mathcal{A}}_{ii}\cdot {\mathrm{S}}_i,$$

$$\mathcal{H}=\frac{1}{2}\sum _{n=1,2}\sum _{(i,j{)}_n}\{J_n^{\prime }{\mathrm{S}}_i\cdot {\mathrm{S}}_j+K_n{S}_i^{\gamma }{S}_j^{\gamma }+{\mathrm{D}}_{n,ij}\cdot ({\mathrm{S}}_i\times {\mathrm{S}}_j)\}+\sum _i{A}_{zz}(S_i^z{)}^2.$$

$$Q=\frac{1}{4\pi }\int \mathrm{m}\cdot (\frac{\partial \mathrm{m}}{\partial x})\times (\frac{\partial \mathrm{m}}{\partial y})dxdy,$$

$$E=E_{spin}+E_{nonspin}=E_0+{\mathrm{S}}_1\cdot {\mathcal{J}}_{12}\cdot {\mathrm{S}}_2+{\mathrm{S}}_1\cdot {\mathrm{K}}_1+{\mathrm{S}}_2\cdot {\mathrm{K}}_2$$

$$J_{xy}=\frac{E_1-E_2-E_3+E_4}{4S^2}$$

$$Q=\frac{1}{4\pi }\sum _l{A}_l$$

$$\cos (\frac{A_l}{2})=\frac{1+{\mathbf{n}}_i\cdot {\mathbf{n}}_j+{\mathbf{n}}_j\cdot {\mathbf{n}}_k+{\mathbf{n}}_k\cdot {\mathbf{n}}_i}{\sqrt{2(1+{\mathbf{n}}_i\cdot {\mathbf{n}}_j)(1+{\mathbf{n}}_j\cdot {\mathbf{n}}_k)(1+{\mathbf{n}}_k\cdot {\mathbf{n}}_i)}}$$

$$\mathrm{sign}(A_l)=\mathrm{sign}[{\mathbf{n}}_i\cdot ({\mathbf{n}}_j\times {\mathbf{n}}_k)]$$

Full text

{affiliations}

Physics Department and Institute for Nanoscience and Engineering, University of Arkansas, Fayetteville, Arkansas 72701, USA

School of Physics and Materials Engineering, Hefei Normal University, Hefei 230601, P. R. China

Key Laboratory of Computational Physical Sciences (Ministry of Education), State Key Laboratory of Surface Physics, and Department of Physics, Fudan University, Shanghai, 200433, China

Collaborative Innovation Center of Advanced Microstructures, Nanjing 210093, P. R. China

Magnetic skyrmions are nano-scale spin clusters with topological stability, and are promising for advanced spintronics [1, 2]. One requirement toward such applications is that the hosting materials should be thin films, so that the nano size of skyrmions can be taken full advantage of. Besides previous studies on bulk MnSi [3, 4, 5, 6], recent works focused on ultrathin films, such as FeGe [7, 8] and rare-earth ion garnet [9, 10], which both take advantage of the Dzyaloshinskii-Moriya interaction (DMI) arising from the heavy metal substrate. However, no skyrmionic state has ever been reported to intrinsically exist in free-standing monolayers, to the best of our knowledge, while two-dimensional (2D) semiconducting magnets, such as monolayer CrI3 [11], are recently attracting much attention due to their novel physics and rich applications [12]. The ferromagnetic monolayer CrI3 crystalizes in honeycomb lattice made of edge-sharing octahedra. Its ferromagnetic order is stabilized by an out-of-plane anisotropy [11], which arises from single ion anisotropy (SIA) and Kitaev-type exchange coupling that both result from the SOC of its heavy ligands [13, 14]. However, the ingredient DMI is absent between the most strongly coupled first nearest neighbor (1st NN) Cr-Cr pairs, because the inversion center between the two Cr atoms prevents its existence [15]. Interestingly, very recent theoretical study proposed the application of electric field to break the inversion center and induce DMI in monolayer CrI3 [16]. Although this clever method leads to CrI3 monolayers becoming closer to adopt a skyrmion phase, the weak effects of electric field in generating DMI, as well as the rather strong out-of-plane anisotropy, hinders the actual creation of skyrmions in this system.

Here we propose a more effective approach that consists in fabricating Janus monolayers of chromium trihalides Cr(I,X)3 (X = Br, Cl). One example of Janus monolayers is the transition metal dichalcogenides MoSSe, which originates from the well-known monolayer MoS2 but with one layer of S being substituted by Se. Such Janus monolayer is experimentally achievable, since MoSSe has been reproducibly obtained using different methods [17, 18], which strongly suggests the feasibility of also creating Cr(I,X)3 Janus monolayers. In this manuscript, we apply density functional theory (DFT) and parallel tempering Monte Carlo (PTMC) simulations to study the magnetic interactions and to investigate skyrmionic states of two Janus monolayers, namely Cr(I,Br)3 and Cr(I,Cl)3, as well as the prototype CrI3 (see structural schematics in Figs. 1a and 1b). As we will show, these Janus monolayers exhibit not only strong enough DMI, but also a decrease in magnetic anisotropy energy (MAE), which both contribute to stabilizing skyrmionic states. We will also demonstrate that the presence of Kitaev interaction and the application of magnetic field both benefit to the emergence of skyrmions in these Cr(I,X)3 Janus monolayers.

Results

Hamiltonian and magnetic parameters

We consider the following Hamiltonian for describing the magnetic interactions of the Cr(I,Br)3 and Cr(I,Cl)3 Janus monolayers, as well as the prototype CrI3 monolayer, using the generalized matrix form up to the second nearest neighbors (2nd NN) for the spins, $\bm{{\rm S}}$:

[TABLE]

where the first term is the exchange coupling that runs over all first $(i,j)_{1}$ and second $(i,j)_{2}$ NN Cr pairs, respectively, and the second term represents SIA that runs over all Cr sites. $\mathcal{J}$ and $\mathcal{A}$ are 3$\times$3 matrices, for which the elements can be extracted using the four-state energy mapping method and density-functional theory (DFT) (see Method section and Ref. [13, 19, 20, 21] for details). The $\mathcal{J}$ matrices can be further decomposed into an isotropic parameter $J^{\prime}$ (to be defined later), the symmetric Kitaev term $K$ [13] and the antisymmetric DMI vector $\bm{{\rm D}}$ [19]; and the $\mathcal{A}$ matrix simply reduces to its $A_{zz}$ components by symmetry. The Hamiltonian can then be rewritten as

[TABLE]

where $S^{\gamma}$ is the $\gamma$ component of S in a local basis {$\alpha,\beta,\gamma$} that is related to the Kitaev interaction. Such local basis, as well as the direction of DMI vector, are schematized in Figs. 1c and 1d, and will be discussed in detail in the following paragraphs.

Note that the Kitaev interaction is a bond-wised anisotropic exchange coupling that was first proposed by A. Kitaev [22] and then found in Na2IrO3 by G. Jackeli and G. Khaliullin [23]. A non-negligible Kitaev coefficient has also been previously identified in CrI3 and CrGeTe3 and was found to be crucial for determining and understanding the different magnetic anisotropies of these two latter compounds [13]. A Kitaev term was also needed in order to reproduce the temperature evolution of magnetization in CrBr3 few layers [24] and is the key ingredient to create quantum spin liquid states (when the $J^{\prime}$ parameter is vanishing [25]).

In the present case, for the 1st NN of Cr(I,X)3, as well as that of CrI3, diagonalization of the symmetric part of $\mathcal{J}_{1}$ matrix leads to (i) nearly degenerate lower-in-energy $J^{\prime}_{1,\alpha}$ and $J^{\prime}_{1,\beta}$, where the eigenvectors $\alpha_{1}$ and $\beta_{1}$ lie in the Cr2L2 (L for ligand) plane; and (ii) higher-in-energy $J^{\prime}_{1,\gamma}$ being associated with the perpendicular $\gamma_{1}$ axis (see Figs. 1c and 1d). Similarly for the 2nd NN, the low-energy easy plane is spanned by the perpendicular 3rd NN Cr pair ($\alpha_{2}$ axis) and the two most related ligands ($\beta_{2}$ axis); the hard axis ($\gamma_{2}$ axis) is thus perpendicular to the easy plane (see schematics in Figs. 1c and 1d). For both the 1st and the 2nd NN, the aforementioned isotropic exchange parameter is thus defined as $J^{\prime}=(J^{\prime}_{\alpha}+J^{\prime}_{\beta})/2$ and the Kitaev coefficient as $K=J^{\prime}_{\gamma}-J^{\prime}$. Note that we further define another isotropic parameter as $J=(J_{xx}+J_{yy}+J_{zz})/3$, which is equivalent to the $J$ coefficient involved in the commonly used $J\bm{{\rm S}}_{i}{\cdot}\bm{{\rm S}}_{j}$ term. Such isotropic $J$ is thus different from $J^{\prime}$ and its ratio with the DMI parameter ($|D/J|$) can provide information about the existence of magnetic skyrmions or not.

Furthermore, the directions of the DMI vectors are also distinct between 1st and 2nd NN. Specifically and as further illustrated in Figs. 1c and 1d, the $\bm{{\rm D}}_{1}$ vectors of all Cr(I,X)3 are basically parallel to $\gamma_{1}$, which satisfy Moriya’s rule [15] and thus testify the accuracy of our calculations. On the other hand, the $\bm{{\rm D}}_{2}$ vectors all stay in the ($\alpha_{2},\beta_{2}$) easy plane and are close to the $\beta_{2}$ axis.

Let us now look at the precise values of the magnetic parameters obtained from DFT and four-state method [13, 19, 20, 21]. As shown in Table I, both $J_{1}$ and $J_{2}$ yield negative values, which imply ferromagnetism (FM) for all investigated systems, which is consistent with the measured ferromagnetism of CrX3, with X = Cl, Br and I [26, 24, 11]. The $J_{1}$, $K_{1}$ and $A_{zz}$ parameters all decrease in magnitude when the X ion of Cr(I,X)3 varies from I to Cl, via Br. This decrease in $J_{1}$ results from the shrinking of lattice constants (leading to shortened Cr-Cr distance and enhanced direct antiferromagnetic exchange), while that of $K_{1}$ and $A_{zz}$ root in the weakening of the SOC strength (as consistent with the location of I, Br and Cl in the periodic Table). On the other hand, the $J_{2}$ and $K_{2}$ coefficients show no significant changes with the X ion. Moreover, the Janus monolayers Cr(I,Br)3 and Cr(I,Cl)3 exhibit remarkable $D_{1}$ values of 0.270 meV and 0.194 meV, respectively, while the prototype CrI3 has no finite $D_{1}$, as aforementioned. Such values, altogether with $J_{1}$, yield large $|D_{1}/J_{1}|$ ratios of 0.150 for Cr(I,Br)3 and 0.194 for Cr(I,Cl)3, which are within the typical range of 0.1-0.2 known to generate skyrmionic phases [2]. Such $|D_{1}/J_{1}|$ ratios are much larger than the reported value 0.071, obtained when applying an extremely large 2 V/nm electric field to CrI3 monolayer [16]. Notably, our three investigated systems also show non-negligible $D_{2}$ values, which leads to $|D_{2}/J_{2}|$ ratios comparable to their 1st NN counterparts, as detailed in Table I. The significant $|D_{1}/J_{1}|$ and $|D_{2}/J_{2}|$ ratios of Cr(I,Br)3 and Cr(I,Cl)3 demonstrate that the fabrication of Janus monolayers is an effective way to induce considerable DMI and is thus promising to create magnetic skyrmions.

Another important factor that affects the formation of skyrmions is the magnetic anisotropy energy (MAE), which is defined as the energy difference between out-of-plane FM (zFM) and in-plane FM (xFM) states, $\epsilon=E_{zFM}-E_{xFM}$, where $\epsilon<0$ ($>0$, respectively) favors out-of-plane (in-plane, respectively) FM. Here, for Cr(I,X)3 and CrI3, the total anisotropy is the result of the interplay between Kitaev interaction and SIA [13] and is calculated using the model energy of Eq. (2) with the parameters of Table 1. As shown in Fig. 2a, CrI3, which has an $\epsilon$ = -0.693 meV/Cr, favors strong out-of-plane anisotropy and Cr(I,Br)3, with $\epsilon$ = -0.218 meV/Cr, exhibits mild out-of-plane anisotropy. In contrast, Cr(I,Cl)3 shows weak in-plane anisotropy since $\epsilon$ = 0.111 meV/Cr (note that Fig. 2a also displays the $\epsilon$ directly calculated from DFT, that agree rather well with those obtained by the model – which thus further testifies the accuracy of our model Hamiltonian). As we will see latter, such differences in anisotropy between Cr(I,Br)3 and Cr(I,Cl)3 plays a role on the morphology of the skyrmions.

Intrinsic skyrmionic states in Cr(I,Br)3 Janus monolayer

As detailed in the Method section, parallel tempering Monte Carlo (PTMC) simulations using the Hamiltonian of Eq. (1) are performed over a 50$\times$50$\times$1 supercell to find spin structures with low energies. These latter spin structures are then further relaxed with the conjugate gradient (CG) method, until arriving at energy minima. Such optimization scheme guarantees the converged spin structures to be either metastable states or the ground state at the temperature of zero Kelvin.

Practically, the ground state of Cr(I,Br)3 is determined to be out-of-plane ferromagnetism (zFM), of which the energy is set to be zero. The first metastable state (the convention of terminology here is that the $n$th metastable state has the $(n+1)$th lowest energy) of Cr(I,Br)3 is an out-of-plane cycloid structure, as shown in Fig. 3a, with an energy of 0.056 meV/Cr. Its propagation direction is symmetrically equivalent along $\bm{{\rm(a+2b)}}$, $\bm{{\rm(2a+b)}}$ and $\bm{{\rm(b-a)}}$. Such cycloid structure yields no finite topological charge $Q$ that is defined as

[TABLE]

where $\bm{{\rm m}}$ is the unit vector lying along the local magnetic moment’s direction.

On the other hand, skyrmionic states with finite $Q$ are found as the second and even higher metastable states. As shown in Fig. 3b, anomalous spin patterns occur at the boundary of dark and bright zones, which are actually the domains with spins having “in” and “out” out-of-plane components , respectively, i.e., z-components that are negative and positive, respectively, along the z-axis (note that, on the other hand, the arrows displayed in Fig.3 represent the in-plane components of the magnetic dipoles). We name such boundary as domain wall. The novel spin structure of Fig. 3b renders $Q=1$ and is thus topologically identical to the common bubble-like skyrmion [3, 4, 5, 6]. Such skyrmionic state has an energy of 0.069 meV/Cr (with respect to the ground state) and is also 0.013 meV/Cr higher than the cycloid state in our used 50$\times$50$\times$1 supercell. Moreover, it is found that either the upper or the lower domain wall can host multiple skyrmions. For example, state possessing two skyrmions ($Q=2$) in the supercell can have (i) one skyrmion at each domain wall (Fig. 3c) or (ii) two skyrmions at the same domain wall (Fig. 3d). Interestingly, such two cases possess degenerate energy of 0.082 meV/Cr. In fact, the energy is a linearly increasing function with simply the number of skyrmions, as shown in Fig. 2c, with the maximum number of skyrmions existing within the 50$\times$50$\times$1 supercell being four (see Fig. 3e for that four-skyrmion state).

Furthermore, it is found that the domain wall can host not only skyrmions, but also anti-skyrmions. As a matter of fact and as shown in Fig. 3f, (i) the spins rotate clockwise from left to right along the lower domain wall, which renders $Q=2$ and characterizes two skyrmions, as similar to the case of Fig. 3d; (ii) but, in contrast, the spins rotate anticlockwise from right to left along the upper domain wall, which yields $Q=-2$ and indicates two antiskyrmions. The total topological charge of Fig. 3f is thus $Q=0$. Comparing Fig. 3e and 3f, we find that either domain wall (upper or lower) can host either skyrmions or anti-skyrmions. The energies of the two spin configurations in Figs. 3e and 3f are further found to be exactly the same, indicating that skyrmion and anti-skyrmion are energetically degenerate, which has never been reported, to the best of our knowledge. Besides, the existence of antisyrmions is actually rather rare in surface or interface systems [27]. It is also found here that skyrmions and anti-skyrmions can not exist at the same domain wall, since they tend to annihilate with each other which would in fact lead to the transformation to the cycloidal state (resulting in $Q=0$).

Magnetic field induced skyrmion states in Cr(I,Cl)3 Janus monolayer.

Let us now turn our attention to the Cr(I,Cl)3 system. As consistent with its positive MAE, the ground state of Cr(I,Cl)3 is determined to be in-plane zigzag-canted FM, as shown in Fig. 4a. The spin vectors at the A and B sites of the honeycomb lattice make an angle of 6.7∘. Such spin canting is further evidenced by DFT calculations on a unit cell, as a moment of 0.3 $\mu_{B}$ along the $y$-direction emerges from the initial FM state fully lying along the x-direction (xFM), with an energy lowering of 0.024 meV/Cr with respect to that xFM state (note that the resulting angle between spins is equal to 6.2∘ in DFT calculations, which compares rather well with the aforementioned 6.7∘ given by numerical simulations using the magnetic Hamiltonian of Eq. (2)). Such zigzag-canted FM state is degenerate when the main component of the magnetization is along the $\bm{{\rm a}}$ direction (which is precisely our x-direction), or along $\bm{{\rm b}}$ and $\bm{{\rm(a+b)}}$ directions.

The first metastable state is a mostly in-plane cycloidal structure superimposed with small out-of-plane components forming a spin wave, as shown in Fig. 4b. Note that such spin arrangement is reminiscent of the well-known complex spin organization of BiFeO3 that consists in the coexistence of an in-plane cycloidal structure and an out-of-plane spin density wave [28]. Such cycloidal structure in Cr(I,Cl)3 has an energy higher by 0.026 meV/Cr from the ground state and is degenerate for propagation directions being along $\bm{{\rm(a+2b)}}$, $\bm{{\rm(2a+b)}}$ or $\bm{{\rm(b-a)}}$ directions. Interestingly, skyrmionic states are not “naturally” stable in Cr(I,Cl)3 within our Hamiltonian using the original DFT-derived magnetic parameters. On the other hand, one can stabilize skyrmionic states there by applying an out-of-plane magnetic field $B$ in order to compensate the in-plane anisotropy. It is numerically found that a critical field of 0.65 T tunes the MAE from positive to negative. However, the magnetic field ($<$ 2 T) does not polarize the spins to the fully out-of-plane direction, but rather renders the in-plane zigzag-canted FM now possessing an out-of-plane component, as shown in Fig. 4c. The previously metastable cycloid state can not survive and transform to such FM state when $B>$ 0.2 T. Practically, skyrmions are numerically found to be stabilized as metastable states when the magnetic field ranges between 0.5 T and 1.3 T. For instance and as shown in Fig. 4d for a magnetic field of 0.8 T, one skyrmion, with a very low energy of 0.015 meV/Cr, is created as a metastable state, rendering $Q=1$. Such skyrmion is quite similar to the common bubble-like skyrmions [3, 4, 5, 6], but also possesses some unique characteristics. For instance, it is not isotropic in the plane and there is always a bright (spin-out) zone being at one side of the dark (spin-in) skyrmion center, as shown in Figs. 4d-f. In fact, up to three skymions can be stabilized within the 50$\times$50$\times$1 supercell and the energy also linearly increases with the skymion numbers in Cr(I,Cl)3, as shown in Fig. 2d. Interestingly, we further find that, wherever the skyrmions are initialized in the supercell, after optimization, their core tend to be aligned equivalently along $\bm{{\rm a}}$, $\bm{{\rm b}}$ and $\bm{{\rm(a+b)}}$ directions, which are the degenerate directions for the in-plane zigzag-canted FM state, as shown in Figs. 4e and 4f. Such well organized structure may be of benefit for potential application in memory devices.

Discussion and Conclusion

Note that we also studied the effects of Kitaev interactions in our investigated systems. It is found that, if the Kitaev terms are turned off in our considered Hamiltonian, the skyrmionic state can no longer be stabilized in Cr(I,Br)3 Janus monolayers. As shown in Fig. 2b, $K_{1}$ (respectively, $K_{2}$) results in negative (respectively, positive) MAE and thus contributes to out-of-plane (respectively, in-plane) anisotropy. The interplay between Kitaev terms and SIA is responsible for the total anisotropy, which plays an important role in the morphology of skyrmions. For instance, we numerically found that arbitrarily modifying the MAE gives the following results: (i) a mild out-of-plane anisotropy benefits the formation of skyrmions, while (ii) a too small anisotropy results in skyrmions with large diameter, as shown in Fig. S1. Moreover, the frustration arising from Kitaev interaction [13] adds to the disorder of the system, which also facilitate the formation of skyrmions. We also studied an analogous system, that is Cr(Br,Cl)3. It exhibits in-plane zigzag-canted FM ground state and a metastable in-plane cycloid, as well as magnetic field induced skyrmion states, which makes it very much similar to Cr(I,Cl)3.

In summary, we proposed the fabrication of Cr(I,X)3 (X = Br, Cl) Janus monolayers to induce large DMI and subsequent magnetic skyrmion states. By combining DFT and MC simulations, we find that Cr(I,Br)3 can intrinsically host metastable skyrmionic phases, while a skyrmionic state of Cr(I,Cl)3 can be stabilized by applying an out-of-plane magnetic field. Our study thus suggests a feasible approach to create skyrmions in semiconducting magnets consisting of chromium trihalides Janus monolayers. Such presently predicted skyrmionic phases are not only useful for memory and logic devices, but can also be promising for energy storage using topological spin textures [29].

{methods}

DFT calculations are performed using the Vienna ab-initio simulation package (VASP) [30]. The projector augmented wave (PAW) method [31] is employed with the following electrons being treated as valence electrons: Cr 3$p$, 4$s$ and 3$d$, I 5$s$ and 5$p$, Br 4$s$ and 4$p$, and Cl 3$s$ and 3$p$. The plane wave energy cut off is chosen to be 350 eV. The length of the $c$ axis of all investigated systems is set to be 19.807 Å , which yields thick enough vacuum layers of more than 16.5 Å. The local density approximation (LDA) [32] is used, with an effective Hubbard $U$ parameter chosen to be 0.5 eV for the localized 3$d$ electrons of Cr ions [13]. $k$-point meshes are chosen such as they are commensurate with the choice of 4$\times$4$\times$1 for the unit cell (that contains 10 atoms). The choices on LDA, Hubbard $U$ and $k$-point meshes were previously tested to be accurate [13]. The Hellman-Feynman forces are taken to be converged when they become smaller than 0.001 eV/Å on each ion.

The four-state energy mapping method is applied to obtain the elements of the $\mathcal{J}$ and $A_{zz}$ matrices [20, 21, 13, 19]. The four-state method considers one specific magnetic pair (or site) at a time. Without loss of generality, the total energy can be written as:

[TABLE]

where $\bm{{\rm S}}_{1}{\cdot}\mathcal{J}_{12}{\cdot}\bm{{\rm S}}_{2}$ is the exchange coupling between the magnetic site 1 and site 2; $\bm{{\rm S}}_{1}{\cdot}\bm{{\rm K}}_{1}$ ($\bm{{\rm S}}_{2}{\cdot}\bm{{\rm K}}_{2}$, respectively) represents the coupling between site 1 (site 2, respectively) and all the other magnetic sites that are different from sites 1 and 2; and the energy from interactions among all these other magnetic sites, together with non-magnetic energy $E_{nonspin}$, is gathered in $E_{0}$. To calculate, e.g., the $J_{xy}$ component of $\mathcal{J}_{12}$, four spin states are considered for the {Cr1,Cr2} pair: {$(S~{}0~{}0),(0~{}S~{}0)$}, {$(S~{}0~{}0),(0~{}-S~{}0)$}, {$(-S~{}0~{}0),(0~{}S~{}0)$} and {$(-S~{}0~{}0),(0~{}-S~{}0)$}; while other spins are perpendicular to those of Cr1 and Cr2. The $J_{xy}$ can then be constructed with the energies obtained from DFT calculations using the following equation:

[TABLE]

Such four-state method is accurate when using a large enough supercell (to prevent the coupling between some sites present in the supercell and periodically repeated sites from happening). For example, SIA and the 1st NN parameters are calculated using a 2$\times$2$\times$1 supercell and the 2nd NN coefficients using a 3$\times$3$\times$1 supercell. See Ref. [20, 21, 13, 19] for more details.

Parallel tempering Monte Carlo (PTMC) simulations and conjugate gradient (CG) method. PTMC simulations with heat bath algorithm [33] are performed using the Hamiltonian of Eq. (1) to update the spin structures. The results shown in the manuscript are based on a 50$\times$50$\times$1 supercell, while a 54$\times$54$\times$1 supercell is also adopted, which leads to similar skyrmion states. 200,000 MC sweeps are performed at each temperature. After the MC simulations, a CG method [34] is applied to further optimize the spin configuration. Specifically, the directions of spins are described with $(\theta_{i},\phi_{i})$, which are independent variables. Such $(\theta_{i},\phi_{i})$ can thus be optimized locally by the CG method to minimize the force on each spin. The energy convergence criteria is set to be 1E-6 eV. The combination of PTMC simulations and the CG method guarantees that our reported magnetic configurations are all located at global/local energy minima.

Calculation of topological charge $Q$. Although Eq. (3) is commonly used to define topological charge $Q$, it is not convenient to employ it to calculate $Q$ for discrete lattice of spins $\bm{n}_{i}$ $(n_{i,x},n_{i,y},n_{i,z})$. Consequently, we alternatively adopt the definition of Berg and L$\rm\ddot{u}$stcher [35],

[TABLE]

[TABLE]

[TABLE]

where $l$ runs over all elementary triangles made of neighboring spin sites and $A_{l}$ is the solid angle formed by the three spin vectors $\bm{n}_{i}$, $\bm{n}_{j}$ and $\bm{n}_{k}$ of the $l$th triangle. Note that the triangles should cover all lattice area with no overlap and $\bm{n}_{i}$, $\bm{n}_{j}$ and $\bm{n}_{k}$ should be anticlockwisely ordered. See details for exceptional cases in Refs. [35, 36].

{addendum}

All data generated or analysed during this study are included in this published article (and its supplementary information files).

This work is supported by the Office of Basic Energy Sciences under contract ER-46612. H.X. is supported by NSFC (11374056), the Special Funds for Major State Basic Research (2015CB921700), Program for Professor of Special Appointment (Eastern Scholar), Qing Nian Ba Jian Program, and Fok Ying Tung Education Foundation. J.F. acknowledges the support from Anhui Provincial Natural Science Foundation (1908085MA10). The Arkansas High Performance Computing Center (AHPCC) is also acknowledged.

C.X. and J.F. contribute equally to this work. The idea of this paper is conceived during the discussion of all authors. C.X. performed the calculations. H.X. and L.B. supervised this work. C.X. wrote the first draft of the manuscript. All authors contributed to the discussion of the results and comments on the first draft of the manuscript.

See supplementary information for details about our predictions.

The authors declare that they have no competing financial interests.

Correspondence should be addressed to H.X. (email: hxiang@fudan.edu.cn) or L.B. (email: laurent@uark.edu).