Skyrmions at vanishingly small Dzyaloshinskii–Moriya interaction or zero magnetic field

Ever since its first realization [1], the magnetic skyrmion has drawn a huge attention because of its promising applications [2] such as memory device, logic gate, writing and deleting information, and emergent properties [3, 4] such as topological Hall effect. The paradigm [3–5] of chiral magnets describes that the nearest neighbor ferromagnetic exchange interaction, J₁, which favors a ferromagnetic ground state competes with the DMI energy, D, arising due to inversion asymmetry. As the DMI tends to make neighboring spins noncollinear, the ground state transforms into a one-dimensional spin-spiral (SS) structure with pitch length λ_D = 2πaJ₁/D where a is the lattice constant. The effective Zeeman energy, H, favors spin orientation along $+\hat{z}$ direction and thus the SS structure is converted into a skyrmion with topological quantum number Q = 1 in the ferromagnetic background. For further increase of H, the size of a skyrmion becomes shorter [3, 4] and eventually diminishes into the background of out-of-plane ferromagnet. The DMI of chiral magnets is quintessential for producing nanoscale SS and skyrmion structures. The out-of-plane magnetic anisotropy, A, tends to orient all magnetic moments along perpendicular to plane and thus the SS becomes unstable for large A in favor of out-of-plane ferromagnet even at H = 0. Finally, an in-plane magnetic anisotropy tends to orient the magnetic moments in a plane and by the influence of the DMI, meron lattice with Q = 1/2 is produced [6–8] before it becomes a planar ferromagnet. The out-of-plane tilting angle of this ferromagnet increases till the structure becomes out-of-plane ferromagnet as H is increased.

1.1. Motivation and model

In contrast to the above standard paradigm, some of the experiments report SS structures with much shorter [9–13] pitch than λ_D, the skyrmions [9, 14–24] at zero magnetic field, and the skyrmions in centrosymmetric [25–27] systems where the DMI is absent. The density functional theory (DFT) data [10, 11] of spin-wave dispersion supports this short-pitched SS structure. Our aim is to introduce a theory in which all these novel magnetic structures along with the usual structures mentioned in the previous section, can be achieved in a single footing. To this end, we consider nearest neighbor positive biquadratic exchange interaction [28, 29] along with the conventional ferromagnetic exchange interaction. On top of the above mentioned experimentally observed phases, our model predicts several other nontrivial magnetic structures that may be experimentally realized.

The presence of the biquadratic exchange interaction in various magnetic systems has been found [30–34] earlier. It generally arises due to the super-exchange mechanism [28] that may also be found in the fourth order [30] of the expansion parameter t/U (where t is the hopping integral and U is the strength of onsite Hubbard interaction) and thus usually smaller than the bilinear exchange interaction which arises in the quadratic order. However, as shown in reference [30], the strength of bilinear term decreases due to fourth order correc-tion. In the antiferromagnetic Fe/Ru(0001) system, the DFT calculation [30] finds the ratio between the strength of biqua-dratic and bilinear terms, J₂/J₁ ∼ 0.65, where J₂, J₁ < 0 (see equation (1) for the sign convention), indicating that the former cannot be ignored. To the best of our knowledge, no ab initio calculation has yet found the opposite signs of the exchange interactions for ferromagnetic systems, i.e., J₁, J₂ > 0 which we on the other hand consider for our model. The microscopic calculation in reference [30] indeed suggests that J₂ can be positive and J₁ is also positive provided t/U ≲ 1/2. We believe that the DFT calculations will reveal such a possibility specially in some of the centrosymmetric systems where skyrmions are observed. We envisage this because as these systems have very weak DMI, there must be some compensating interaction which will frustrate the system from being a ferromagnet. A model with Ruderman, Kittel, Kasuya, and Yosida (RKKY) type frustrating bilinear exchange interaction is already proposed [10, 11] in literature for explaining the DFT data of non-ferromagnetic spin-wave dispersion, but it will not be appropriate for nanoscale magnetic structures as the interaction is long-ranged. Indeed, the importance of four-spin interaction whose two-site contribution is nothing but the biquadratic exchange interaction have been attributed to the atomic-scale skyrmions in Fe/Ir(111) system [13].

We consider a square lattice, as an illustration, of a magnetic system with Hamiltonian $\mathcal{H}={\mathcal{H}}_{\mathrm{E}\mathrm{X}}+{\mathcal{H}}_{\mathrm{D}}+{\mathcal{H}}_{\mathrm{A}}+{\mathcal{H}}_{\mathrm{Z}}$, where the spin-exchange Hamiltonian

$$\begin{equation}{\mathcal{H}}_{\mathrm{E}\mathrm{X}}=\sum\limits _{\langle ij\rangle }\left[-{J}_{1}\left({\boldsymbol{m}}_{i}\cdot {\boldsymbol{m}}_{j}\right)+{J}_{2}{\left({\boldsymbol{m}}_{i}\cdot {\boldsymbol{m}}_{j}\right)}^{2}\right]\end{equation} \tag{ 1 }$$

with nearest neighbor bilinear and biquadratic exchange couplings (J₁, J₂ > 0) and magnetization unit vector m _i at ı^th site, inversion asymmetry induced Dzyaloshinskii–Moriya interaction ${\mathcal{H}}_{\mathrm{D}}=D{\sum }_{\langle ij\rangle }\left(\hat{z}{\times}{\hat{r}}_{ij}\right)\cdot \left({\boldsymbol{m}}_{i}{\times}{\boldsymbol{m}}_{j}\right)$ with strength D for the systems with C_nv symmetries, energy due to magnetic anisotropy ${\mathcal{H}}_{\mathrm{A}}=A{\sum }_{i}{\left({m}_{i}^{z}\right)}^{2}$ with positive (negative) sign of A representing in-plane (out of plane) anisotropy, and Zeeman energy due to application of magnetic field ${\mathcal{H}}_{\mathrm{Z}}=-H{\sum }_{i}{m}_{i}^{z}$ represented by energy parameter H including both applied and demagnetization fields.

We find that the biquadratic exchange energy, beyond a threshold value, minimizes spin-wave dispersion energy at a nonzero momentum signaling the possibility of noncollinear magnetic structures. The pitch length of the SS structure is found to be much smaller than λ_D. Employing the method of simulated annealing by further inclusion of DMI, magnetic anisotropy and magnetic field in the model, we obtain various noncollinear magnetic structures that are experimentally observed but most of which are not yet theoretically understood for the presence of all or the absence of one or more of these parameters: spin-spiral [9–13, 35]; skyrmion lattice [3, 16, 36]; isolated skyrmions [9, 14, 35, 37, 38]; mixed phase of broken spirals, chiral bubbles, and skyrmions; [21, 39, 40] meron lattice [6, 7, 41].

A naive back of the envelope calculation suggests that ${\mathcal{H}}_{\text{EX}}$ is minimized when the relative orientation of two neighboring magnetic moment is ±cos⁻¹(J₁/2J₂) for J₁ < 2J₂ and zero otherwise. The degeneracy in the mode of orientation (clockwise or anticlockwise) when J₂ > J₁/2 is broken for an infinitesimal D which favors one kind of orientation only, depending on its sign. We thus expect the spin-spiral ground state with pitch length λ_EX ∼ 2πa/cos⁻¹(J₁/2J₂) in the limit of vanishingly small D. λ_EX decreases with the increase of J₂/J₁ and it is about 11 lattice sites for J₂ = 0.6J₁. It can further decrease with the increase of D. We note that even for vanishingly small DMI relevant for centrosymmetric systems, SS structure with pitch-length of a few lattice sites is possible. For the same reason, we expect, contrary to the paradigm, the formation of small size skyrmions even in centrosymmetric systems but with J₂ > J₁/2 when magnetic field is applied. Does this model also favor stabilizing skyrmions in the absence of magnetic field? J₁ supports J₂ for orienting neighboring spins along same direction, but it cannot break symmetry to orient all spins along a particular direction because J₂ will destabilize. However, if large easy-axis magnetic anisotropy spontaneously creates up (down) spin background then the spiral effect generated by the exchange interactions can produce skyrmions with topological number Q = 1 whose center will have down (up) spin moment. We next describe our detailed results obtained using analytical calculation and numerical simulation, that will be in agreement with the above intuitive description.

We find (see appendix A for details) the spin-wave dispersion from the Hamiltonian ${\mathcal{H}}_{\text{EX}}+{\mathcal{H}}_{\mathrm{D}}$ for k_y = ±k_x (along high-symmetry directions) as

$$\begin{equation}{E}_{k}=-4\left(1+\frac{16}{{\pi }^{2}}{J}_{2}\right)\mathrm{cos}\left({k}_{x}\right)+4{J}_{2}\enspace \mathrm{cos}\left(2{k}_{x}\right)+4D\enspace \mathrm{sin}\left({k}_{x}\right).\end{equation} \tag{ 2 }$$

Henceforth, we assume all length scales are in the unit of a and energy scale in the unit of J₁. Figure 1(a) shows E_k for D = 0, 0.15, 0.8, −0.8 and J₂ = 0.7. We find two degenerate minima at two nonzero k_x (equal in magnitude but opposite in sign) for D = 0. For any arbitrary nonzero magnitude of D, this degeneracy is broken and one global minimum occurs at k_x = −k_min sgn(D), and k_min increases with |D|. Energy dispersion obtained in ab initio studies [10, 11] fit (figure 1(b)) quite well with the analytical form in equation (2). With the input value of J₁ provided in these ab initio studies, we extract J₂, D and λ = 2π/k_min which are tabulated in table 1. The pitch length λ of SS excellently agrees with experiments [10, 11] on Mn/W and Pd/Fe/Ir systems. Figure 1(c) shows the variation of λ⁻¹ with |D|. We note that while λ⁻¹ is zero at D = 0 for J₂ = 0.4, it is nonzero for J₂ = 0.45. To be precise, SS structure is possible even at D = 0 for 0.42 ≲ J₂; the lower bound of J₂ for SS structure is somewhat less than the naive calculation discussed above. As J₂ is increased, the dependence of λ on D decreases and it becomes shorter. Therefore, the energy scale J₂ provides SS structures with shorter pitch length than the same for the presence of D only. For example, if J₁ = 1 and D = 0.05, the standard paradigm supports spin-spiral pitch length of about 125 lattice sites, while our model with J₂ = 0.5 keeping the same values of J₁ and D supports spin-spiral of pitch length 12 lattice sites. However, if J₂ < 0.42, like the canonical model, large DMI is essential to produce spin spiral of similar pitch-length. Therefore, the present model with J₂ > 0.42 has clear advantage (which we show below) over the canonical model for explaining nanoscale magnetic structures for the physical systems with small DMI, for example, thin films made with centrosymmetric crystals.

Zoom In Zoom Out Reset image size

Alternatively, the long ranged (more than 6th nearest neighbor) exchange couplings [10, 11] which may arise due to long-ranged RKKY type interaction seem to reproduce this short pitch of spin-spiral. The model with bilinear and biquadratic nearest neighbor exchange interaction may be as relevant as the model with long-ranged bilinear exchange interaction to the systems that have been studied by these authors. While both the models describe spin-spiral for vanishingly small DMI, it distinguishes the models as the former (latter) provides equal (unequal) change in spin orientation between two neighbors that may easily be understood through their respective continuum versions. However, such a model with long-ranged and multi-parameter exchange couplings is possibly not relevant when the skyrmions produced due to the superposition of spin spirals along high symmetric directions is much shorter in size compared to this range of the interaction. Further, exchange interaction up to second nearest neighbor (a short-range version of RKKY type interaction) seems to suggest skyrmions with Q > 1, instead of Q = 1 [42].

We next perform numerical simulation (see appendix B for its method) of $\mathcal{H}$ for determining ground state magnetic structures in the parameter space of A, D and H for J₂ = 0.7 which supports SS even for small D as discussed above. Various magnetic phases obtained in the simulation are characterized by the estimation of following parameters: local magnetization ${m}_{i}^{z}$, average magnetization $M=\left(1/N\right){\sum }_{i}\;{m}_{i}^{z}$, local chirality ${\rho }_{i}=\left(1/4\pi \right){\sum }_{\delta ={\pm}1}{\boldsymbol{m}}_{i}\cdot \left({\boldsymbol{m}}_{i+\delta \hat{x}}{\times}{\boldsymbol{m}}_{i+\delta \hat{y}}\right)$, and total chirality $\mathcal{R}=\left(1/N\right){\sum }_{i}\;{\rho }_{i}$; local spin asymmetric parameter $\delta {\theta }_{i}=\vert {\theta }_{i+\hat{y}}-{\theta }_{i-\hat{y}}\vert -\vert {\theta }_{i+\hat{x}}-{\theta }_{i-\hat{x}}\vert$, where N is the number of lattice sites. We find various magnetic phases (I–XI as marked in figures 2 and 3) characterized by different ranges of these parameters shown in table 2.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Table 2. Characterization of magnetic phases in terms of various parameters calculated from the ground states. While all the parameters sometimes may be same for the phases X and XI, they can be distinguished through the Fourier decomposition as the former is an ordered structure and the latter is disordered.

Phase	${m}_{i}^{z}$	M	ρi	$\mathcal{R}$	δθi
I	≲1	>0.98	∼0.0	<0.001	0
II	(−1)–(+1)	<0.001	∼0.0	<0.001	≠0
III	(−1)–(+1)	0.1–0.6	0.0–1.0	0.25–0.45	≠0
IV	At the core (−0.98)–(−1.0)	0.45–0.8	At the core 0.95–1.0	0.5–0.9	0
V	At the core (−0.98)–(−1.0)	0.8–0.98	At the core 0.95–1.0	0.04–0.5	0
VI	At the core (−0.98)–(−1.0)	0.3–0.45	0.95–1.0	0.5–0.7	≠0
VII	(−1)–(+1)	0.04–0.16	<0.95	<0.4	≠0
VIII	At the core $\vert {m}_{i}^{z}\vert { >}0.98$	∼0.001	∼0.5	∼0.0	0
IX	∼0.5	≃0.5	≃0.0	<0.001	0
X	(−0.95)–(0.95)	Any value	<0.001	<0.001	≠0
XI	(−0.95)–(0.95)	<0.001	<0.001	<0.001	≠0

3.1. Structures for large DMI

3.2. Structures for vanishingly small DMI

3.3. Structures at zero field

We next analyze magnetic structures at H = 0. Figures 3(a) and (b) respectively show the corresponding phases for A < 0 and A > 0 in D–A plane. Interestingly, apart from expected polarized ferromagnet and SS structures, a mixed phase of broken spin-spirals, chiral bubbles and skyrmions (phase-III) are found for out-of-plane magnetic anisotropy. These skyrmions are likely to be corroborated with the skyrmions observed at zero magnetic field at various experiments [9, 14–22]. Meyer et al [9] have recently reported that skyrmions at zero magnetic field are observed only for the systems with energy dispersion having quartic wavenumber dependence rather than quadratic only. However, the dispersion relation with positive coefficients for quadratic term will not make any qualitative difference. The sign of quadratic and quartic terms ought to be negative and positive respectively. Our dispersion relation (2) for D = 0 at small momenta provides exactly this expectation. The dispersion relation at small momenta reads

$$\begin{align}{E}_{k}& =4\left[-1+{J}_{2}\left(1-\frac{16}{{\pi }^{2}}\right)\right]+\left[2-8{J}_{2}\left(1-\frac{4}{{\pi }^{2}}\right)\right]{k}_{x}^{2}\\ & \quad +\frac{1}{6}\left[-1+16{J}_{2}\left(1-\frac{1}{{\pi }^{2}}\right)\right]{k}_{x}^{4}\end{align} \tag{ 3 }$$

which suggests negative coefficient for quadratic term when J₂ > 0.42, which is consistent with our detailed analysis about the lower bound of J₂ in section 2.

We find that the skyrmions with down as well as up spin at their centers (see zoomed structures of phase-III in figure 3(c)) may simultaneously form in the respective background of up and down spins because both these backgrounds at zero magnetic field are degenerate. In the case of positive anisotropy, apart from the well-known phases of SS and planar ferromagnet, we find a regime of D and A that supports a lattice of merons-a meron with up spin at its center is surrounded by merons with down spin at their centers and vice versa (phase-VIII in figure 3(c)). This lattice has been recently reported in experiments [6, 7, 41], a micromagnetic simulation [7] and also as analytic solution [8] of Euler equation of a single meron followed by physical arguments. For the parameter regime in between planar ferromagnetic and meron lattice phases, we find an SS structure like phase but with the limited range (phase-XI in figure 3(c)) of out-of-plane magnetization, |m_z | < 1. The value of |m_z | increases as we approach closer to meron lattice boundary. However, we find a narrow regime of parameter space towards the planar ferromagnet phase boundary where the structure appears to be the disordered spin island (phase-X in figure 3(c)) phases with |m_z | < 1.

The pitch length of SS structure obtained in this simulation method for various values of J₂ agrees very well with the analytical estimate for an infinite system (figure 1(c)). For increasing accuracy in estimating λ from the simulation for lower |D|, we have considered 128 × 128 lattice. The minimum value of |D| for which SS structure is obtained for a given lattice size is at the threshold value of J₂ ≈ 0.42 which agrees with its analytical estimation.

The spin configuration of some of the magnetic phases may be represented in terms of the truncated Fourier decomposition [42, 48, 49]:

$$\begin{equation}{\boldsymbol{S}}_{i}={\boldsymbol{B}}_{0}+\frac{1}{2}\sum\limits _{\alpha =1}^{2}\left({\boldsymbol{B}}_{\alpha }\enspace {\text{e}}^{\text{i}{\boldsymbol{q}}_{\alpha }\cdot {\boldsymbol{r}}_{i}}+{\boldsymbol{B}}_{\alpha }^{{\ast}}\enspace {\text{e}}^{-\text{i}{\boldsymbol{q}}_{\alpha }\cdot {\boldsymbol{r}}_{i}}\right)\end{equation} \tag{ 4 }$$

where ${\boldsymbol{q}}_{1,2}=\left(q/\sqrt{2}\right)\left(1,\enspace {\pm}1\right)$ are two orthogonal wave vectors representing two high symmetric directions in a square lattice. The uniform ferromagnetic state is represented by ${B}_{0}^{z}=1$ and ${B}_{0}^{x}={B}_{0}^{y}={\boldsymbol{B}}_{\alpha }=0$. We find that the SS structure corresponds to B ₀ = B ₂ = 0 and ${\boldsymbol{B}}_{1}=\left(-\mathrm{i}/\sqrt{2},\enspace -\mathrm{i}/\sqrt{2},\enspace 1\right)$. The direction of propagation of the spin-spiral is along q ₁ which switches over to q ₂ when the sign of DMI changes. Superposition of modulations along both q _α directions with equal amplitude creates the skyrmion structures with circular symmetry. In this case, B ₀ = 0, ${\boldsymbol{B}}_{1}=B\enspace {\text{e}}^{\text{i}{\phi }_{1}}\left(-\mathrm{i}\enspace \mathrm{sin}\enspace \chi ,\mathrm{i}\enspace \mathrm{cos}\enspace \chi ,1\right)$ and ${\boldsymbol{B}}_{2}=B\enspace {\text{e}}^{\text{i}{\phi }_{2}}\left(-\mathrm{i}\enspace \mathrm{cos}\enspace \chi ,-\mathrm{i}\enspace \mathrm{sin}\enspace \chi ,1\right)$ where B = −1/2 satisfying ${S}_{i}^{z}=-1$, ${S}_{i}^{x}={S}_{i}^{y}=0$ at the center of the skyrmions. The magnetic vortices are also formed due to the superposition of both q ₁ and q ₂ but with unequal amplitudes.

By introducing biquadratic nearest neighbor exchange interaction, we comprehensively show numerous exotic magnetic structures such as spin-spiral with unusually short pitch-length, skyrmions at zero magnetic field and/or vanishingly small Dzyaloshinskii–Moriya interaction, and meron lattice at zero magnetic field. We thus have made a compelling case of explaining a large number of recent experiments observing these unusual magnetic structures in a single model. Our phase diagrams will encourage further investigations of several other magnetic structures including these over their respective ranges of parameter space. Our theory should also be relevant for recent observation [50] of spin-spiral in YBaCuFeO₅ with extremely small Dzyaloshinskii–Moriya interaction [51]. We predict that these systems should also host skyrmions and other magnetic structures shown in the phase diagrams.

Our study is based on the model consisting of nearest neighbor ferromagnetic bilinear and a positive biquadratic exchange interactions. Ab initio studies are indeed necessary for investigating such interactions in some of the systems, especially the centrosymmetric systems, where skyrmions are observed.

All data that support the findings of this study are included within the article (and any supplementary files).

In this appendix, we show detailed calculation of the spin-wave dispersion relation for the exchange Hamiltonian:

$$\begin{equation}{\mathcal{H}}_{\mathrm{E}\mathrm{X}}=\sum\limits _{\langle ij\rangle }\left[-{J}_{1}\left({\boldsymbol{m}}_{i}\cdot {\boldsymbol{m}}_{j}\right)+{J}_{2}{\left({\boldsymbol{m}}_{i}\cdot {\boldsymbol{m}}_{j}\right)}^{2}\right].\end{equation} \tag{ A1 }$$

By employing the method discussed in reference [52], we write m_j 's in terms of Bosonic creation and annihilation operators $\left({a}_{j}^{{\dagger}},{a}_{j}\right)$ with the commutation relation $\left[{a}_{j},{a}_{l}^{{\dagger}}\right]={\delta }_{jl}$ as ${m}_{j}^{+}={m}_{j}^{x}+\mathrm{i}{m}_{j}^{y}={\left(2S-{a}_{j}^{{\dagger}}{a}_{j}\right)}^{1/2}{a}_{j}$, ${m}_{j}^{-}={m}_{j}^{x}-\mathrm{i}{m}_{j}^{y}={a}_{j}^{{\dagger}}{\left(2S-{a}_{j}^{{\dagger}}{a}_{j}\right)}^{1/2}$, and ${m}_{j}^{z}=S-{a}_{j}^{{\dagger}}{a}_{j}$ where the magnitude of spin S is assumed to be large for perturbative expansion and finally we take the limit S → 1. As we substitute m _i in equation (A1) by the above Bosonic operators, m _i ⋅ m _j should be replaced by m _i ⋅ m _j − S².

Expressing ${a}_{j}=\frac{1}{\sqrt{N}}{\sum }_{\boldsymbol{k}\;}{\text{e}}^{-\text{i}\boldsymbol{k}\cdot {\boldsymbol{r}}_{j}}{b}_{\boldsymbol{k}}$ and ${a}_{j}^{{\dagger}}=\frac{1}{\sqrt{N}}{\sum }_{\boldsymbol{k}\;}{\text{e}}^{\text{i}\boldsymbol{k}\cdot {\boldsymbol{r}}_{j}}{b}_{\boldsymbol{k}}^{{\dagger}}$ in Fourier representation with $\left[{b}_{\boldsymbol{k}},\enspace {b}_{{\boldsymbol{k}}^{\prime }}^{{\dagger}}\right]={\delta }_{\boldsymbol{k}{\boldsymbol{k}}^{\prime }}$, we find ${\mathcal{H}}_{\mathrm{E}\mathrm{X}}={H}_{0}+{H}_{1}$ (ignoring terms without Bosonic operators) where

$$\begin{align}{H}_{0}& =-2{J}_{1}S\sum\limits _{\boldsymbol{k}}\left(\mathrm{cos}\left({k}_{x}a\right)+\mathrm{cos}\left({k}_{y}a\right)-2\right){b}_{\boldsymbol{k}}^{{\dagger}}{b}_{\boldsymbol{k}}\\ & =-2{J}_{1}S\sum\limits _{\boldsymbol{k}}\zeta \left(k\right){b}_{\boldsymbol{k}}^{{\dagger}}{b}_{\boldsymbol{k}}\end{align} \tag{ A2 }$$

with ζ(k) = cos(k_x a) + cos(k_y a) − 2 and

$$\begin{equation}{H}_{1}=\frac{{S}^{2}}{N}\sum\limits _{\boldsymbol{k},{\boldsymbol{k}}_{1},{\boldsymbol{k}}^{\prime }}\mathcal{J}\left(\boldsymbol{k},{\boldsymbol{k}}_{1},{\boldsymbol{k}}^{\prime }\right){b}_{\boldsymbol{k}}^{{\dagger}}{b}_{{\boldsymbol{k}}_{1}}{b}_{{\boldsymbol{k}}^{\prime }}^{{\dagger}}{b}_{\boldsymbol{k}-{\boldsymbol{k}}_{1}+{\boldsymbol{k}}^{\prime }}.\end{equation} \tag{ A3 }$$

For the purpose of determining dispersion relation, we set k = k ₁ and thus

$$\begin{equation}{\tilde {H}}_{1}={S}^{2}\sum\limits _{\boldsymbol{k},{\boldsymbol{k}}^{\prime }}\mathcal{J}\left(\boldsymbol{k},{\boldsymbol{k}}^{\prime }\right){b}_{\boldsymbol{k}}^{{\dagger}}{b}_{\boldsymbol{k}}{b}_{{\boldsymbol{k}}^{\prime }}^{{\dagger}}{b}_{{\boldsymbol{k}}^{\prime }}\end{equation} \tag{ A4 }$$

with

$$\begin{align}\mathcal{J}\left(\boldsymbol{k},{\boldsymbol{k}}^{\prime }\right)& =4{J}_{2}\left[\mathrm{cos}\left({k}_{x}a\right)\mathrm{cos}\left({k}_{x}^{\prime }a\right)+\mathrm{cos}\left({k}_{y}a\right)\mathrm{cos}\left({k}_{y}^{\prime }a\right)\right]\\ & \quad -\left(4{J}_{2}-{J}_{1}/2\right)\left[\mathrm{cos}\left({k}_{x}a\right)+\mathrm{cos}\left({k}_{y}a\right)\right.\\ & \quad \left.+\mathrm{cos}\left({k}_{x}^{\prime }a\right)+\mathrm{cos}\left({k}_{y}^{\prime }a\right)\right]+\left(8{J}_{2}-2{J}_{1}\right)\end{align} \tag{ A5 }$$

We thus obtain an effective Hamiltonian (quadratic in Bosonic field) as

$$\begin{equation}{H}_{\text{EX}}^{\mathrm{e}\mathrm{ff}}=\sum\limits _{\boldsymbol{k}}\tilde {\mathcal{J}}\left(k\right){b}_{\boldsymbol{k}}^{{\dagger}}{b}_{\boldsymbol{k}}\end{equation} \tag{ A6 }$$

where

$$\begin{equation}\tilde {\mathcal{J}}\left(k\right)=-2{J}_{1}\zeta \left(k\right)+\sum\limits _{{\boldsymbol{k}}^{\prime }}\mathcal{J}\left(\boldsymbol{k},{\boldsymbol{k}}^{\prime }\right)\langle {b}_{{\boldsymbol{k}}^{\prime }}^{{\dagger}}{b}_{{\boldsymbol{k}}^{\prime }}\rangle \end{equation} \tag{ A7 }$$

by putting S = 1. The ensemble average $\langle {b}_{{\boldsymbol{k}}^{\prime }}^{{\dagger}}{b}_{{\boldsymbol{k}}^{\prime }}\rangle$ (see figure 4 schematically) may be obtained using the expression of H₀. To that end, ζ(k) has four-fold symmetry in the Brillouin zone, with one of such independent regime is 0 ⩽ (k_x , k_y ) ⩽ π in which ζ(k) > (<)0, which we denote as ζ⁺(ζ⁻), for k_x + k_y ⩽ (>)π. Therefore,

$$\begin{align}\tilde {\mathcal{J}}\left(k\right)& =-2{J}_{1}\zeta \left(k\right)+\frac{1}{{\pi }^{2}}\left[{\int }_{0}^{\pi }\mathrm{d}{k}_{x}^{\prime }{\int }_{0}^{\pi -{k}_{x}^{\prime }}\enspace \mathrm{d}{k}_{y}^{\prime }\sum\limits _{n}\frac{\mathcal{J}\left(\boldsymbol{k},{\boldsymbol{k}}^{\prime }\right)}{\mathrm{i}{{\Omega}}_{n}-{\zeta }^{-}}\right.\\ & \quad +\left.{\int }_{0}^{\pi }\mathrm{d}{k}_{x}^{\prime }{\int }_{\pi -{k}_{x}^{\prime }}^{\pi }\enspace \mathrm{d}{k}_{y}^{\prime }\sum\limits _{n}\frac{\mathcal{J}\left(\boldsymbol{k},{\boldsymbol{k}}^{\prime }\right)}{\mathrm{i}{{\Omega}}_{n}-{\zeta }^{+}}\right]\end{align} \tag{ A8 }$$

where Ω_n is bosonic Matsubara frequency. We hence find energy dispersion

$$\begin{align}\tilde {\mathcal{J}}\left(\boldsymbol{k}\right)& =-2\left({J}_{1}+\left(16/{\pi }^{2}\right){J}_{2}\right)\left(\mathrm{cos}\left({k}_{x}a\right)+\mathrm{cos}\left({k}_{y}a\right)\right)\\ & \quad +2{J}_{2}\left(\mathrm{cos}\left(2{k}_{x}a\right)+\mathrm{cos}\left(2{k}_{y}a\right)\right)\end{align} \tag{ A9 }$$

at zero temperature, by dropping the constant terms. We note that $\tilde {\mathcal{J}}\left(\boldsymbol{k}\right)$ has four-fold symmetry in the Brillouin zone as the energy is same at (±k_x , ±k_y ). We thus find effective exchange interaction: ${\mathcal{H}}_{\mathrm{E}\mathrm{X}}^{\mathrm{e}\mathrm{ff}}={\sum }_{\boldsymbol{k}}\tilde {\mathcal{J}}\left(\boldsymbol{k}\right){b}_{\boldsymbol{k}}^{{\dagger}}{b}_{\boldsymbol{k}}\equiv {\sum }_{\boldsymbol{k}}\tilde {\mathcal{J}}\left(\boldsymbol{k}\right){\boldsymbol{m}}_{\boldsymbol{k}}\cdot {\boldsymbol{m}}_{-\boldsymbol{k}}$, and DM interaction energy

$$\begin{align}{\mathcal{H}}_{\mathrm{D}\mathrm{M}}& =\mathrm{i}\sum\limits _{\boldsymbol{k}}\left[{\mathcal{D}}_{x}\left(\boldsymbol{k}\right)\left({m}_{\boldsymbol{k}}^{z}{m}_{-\boldsymbol{k}}^{x}-{m}_{\boldsymbol{k}}^{x}{m}_{-\boldsymbol{k}}^{z}\right)\right.\\ & \qquad \left.-{\mathcal{D}}_{y}\left(\boldsymbol{k}\right)\left({m}_{\boldsymbol{k}}^{z}{m}_{-\boldsymbol{k}}^{y}-{m}_{\boldsymbol{k}}^{y}{m}_{-\boldsymbol{k}}^{z}\right)\right]\end{align} \tag{ A10 }$$

Zoom In Zoom Out Reset image size

where ${\mathcal{D}}_{x}\left(\boldsymbol{k}\right)=2D\enspace \mathrm{sin}\left({k}_{x}a\right)$, ${\mathcal{D}}_{y}\left(\boldsymbol{k}\right)=2D\enspace \mathrm{sin}\left({k}_{y}a\right)$, and m _k is Fourier transform of m_i .

Without loss of generality, we assume $\boldsymbol{k}={k}_{x}{\hat{e}}_{x}$ and considering the variation of magnetization in the z–x plane, we find the effective Hamiltonian as

$$\begin{equation}\tilde {\mathcal{H}}=\sum\limits _{{k}_{x}}\left(\begin{matrix}\hfill {m}_{\boldsymbol{k}}^{z},\hfill & \hfill {m}_{\boldsymbol{k}}^{x}\hfill \end{matrix}\right)\left[\begin{matrix}\hfill \tilde {\mathcal{J}}\left({k}_{x}\right)\hfill & \hfill \mathrm{i}{\mathcal{D}}_{x}\left({k}_{x}\right)\hfill \\ \hfill -\mathrm{i}{\mathcal{D}}_{x}\left({k}_{x}\right)\hfill & \hfill \tilde {\mathcal{J}}\left({k}_{x}\right)\hfill \end{matrix}\right]\left(\begin{matrix}\hfill {m}_{-\boldsymbol{k}}^{z}\hfill \\ \hfill {m}_{-\boldsymbol{k}}^{x}\hfill \end{matrix}\right)\end{equation} \tag{ A11 }$$

leading to the energy dispersion

$$\begin{equation}{E}^{{\pm}}\left({k}_{x}\right)=\tilde {\mathcal{J}}\left({k}_{x}\right){\pm}{\mathcal{D}}_{x}\left({k}_{x}\right).\end{equation} \tag{ A12 }$$

Clearly, E⁺(D, k_x ) = E⁻(−D, k_x ) = E⁺(−D, −k_x ). It is thus sufficient that we consider E⁺(D, k_x ) ≡ E(k_x ) only.

We obtain magnetic phase diagrams by performing simulated annealing from a large temperature for the spin model described by $\mathcal{H}$ in a square lattice 32 × 32 size with periodic boundary condition. Some of the key results have been reconfirmed for 128 × 128 size also. The simulation is carried out by adopting standard Metropolis algorithm for updating local spins up to 4 × 10⁶ steps for each temperature. We gradually reduce the temperature in each step of the simulation [53] following the relation T_n+1 = αT_n where T_n is the temperature in the nth step. We set the initial temperature T₀ = 11 and α = 0.99 to reach final temperature T ∼ 5 × 10⁻⁴ in 10³ steps.