Phase diagram of the two-dimensional dipolar Heisenberg model with the Dzyaloshinskii-Moriya interaction and the Ising anisotropy

Ultrathin magnetic films have attracted much attention for applications toward magnetic recording and their unique magnetic properties have drawn scientific interest Heinrich ; Bader . The competition between magnetic anisotropies, short-range exchange and long-range dipolar interactions causes complex magnetic orderings such as spin-reorientation transitions between in-plane and out-of-plane magnetic phases and a variety of stripe patterns Pappas ; Allenspach ; Ramchal ; Won ; Qiu .

The theoretical aspects of these phenomena have been often investigated by using the two-dimensional (2D) dipolar Ising Booth ; MacIsaac-Ising ; Toloza ; Rastelli06 ; Cannas ; Pighin-Ising ; Rastelli ; Vindigni ; Rizzi ; Fonseca ; Ruger ; Horowitz ; Bab ; Leib ; Komatsu1 or Heisenberg Pescia ; Moschel ; Hucht ; MacIsaac1 ; MacIsaac2 ; Bell ; Santamaria ; Rapini ; Whitehead ; Carubelli ; Mol ; Pighin2 ; Pighin ; Mol2 ; Komatsu2 model with a magnetic anisotropy, and the phase diagrams with multiple stripe-ordered phases have been shown in several parameter regions. Reentrant transitions associated with planar ferro, stripe, and paramagnetic phases have also been presented Komatsu1 ; Komatsu2 .

The Dzyaloshinskii-Moriya (DM) interaction plays an important role in systems whose spatial reversal symmetry is broken, and it causes weak ferromagnetism or helimagnetism. Recently a topologically protected magnetic structure called (magnetic) skyrmion has been observed experimentally Muhlbauer ; Yu ; Jonietz11 ; Yu11 ; Schulz12 ; Sampaio13 ; Jiang15 ; Boulle16 , and because of its potential applications to spintronics devices, physical properties of skyrmions have drawn much attention. For skyrmion-lattice phases, the competition between the DM and exchange interactions is important.

Skyrmion-lattice phases are more stabilized in 2D systems (thin films) than in 3D ones Yu . In 2D systems, the phase diagrams associated with skyrmion-lattice phases have been actively investigated by theoretical and computational methods with the use of the Heisenberg model with the DM interaction and with and without several anisotropies, and their phase diagrams in specific parameter regions have been shown Yi ; Kwon ; Banerjee ; Lin ; Rowland ; Nishikawa ; Bernard . However, studies on the effect of the dipolar interaction on the models Kwon ; Bernard have hardly been performed and have not been well understood because of numerical difficulty for treating long-range interactions.

In the present paper, we study the 2D classical Heisenberg model with the DM interaction, Ising anisotropy, and dipolar interaction under zero and finite fields. We show the field vs. Ising anisotropy ($H$-$\eta$) phase diagrams for typical three cases of the strength of the dipolar interaction, i.e., zero, weak, and strong interactions.

We perform a direct simulation of order parameters in the model by using a Monte Carlo method. The most serious difficulty in numerical computation is $O(N^{2})$ ($N$: the total number of spins) computational time originating from the fully-connected $O(N^{2})$ long-range interactions, and we adopt the stochastic cut-off (SCO) $O(N)$ method Sasaki to reduce the computational cost.

The $H$-$\eta$ phase diagram in the ground state without the dipolar interaction was studied using variational approaches for equivalent continuum models. It was shown that in the case of small DM interaction, for small $|\eta|$, triangular skyrmion-lattice and helical phases exist at high and low (including zero) fields, respectively, which are located between the canted ferromagnetic (in-plane magnetic) phase at small $\eta$ and polarized ferromagnetic phase at large $\eta$ Kwon ; Banerjee .

Afterward, Lin et al. presented that a square skyrmion-lattice phase can exist at finite fields and at a specific small region of $\eta$ Lin . This phase is located at higher fields than the helical phase and between the canted ferromagnetic phase at smaller $\eta$ and the triangular skyrmion-lattice phase at larger $\eta$.

In the present paper, we report the following new findings. Without the dipolar interaction, when the DM interaction is the same order of the exchange interaction, the square skyrmion-lattice phase exists in a wider region of the phase diagram including zero field. The dipolar interaction expands the regions of the ordered states whose boundaries shift to the positive $H$ and $\eta$ directions, increases the density of skyrmions, and shortens the pitch length (stripe width) of helical structures. The dipolar interaction induces another type of square skyrmion-lattice phase with a different alignment, and in the weak strength, a reentrance of a vertical (or horizontal) helical phase through a diagonal helical phase takes place with increasing $\eta$. These complex situations cause five kinds of triple points in the phase diagrams.

The outline of the present paper is organized as follows. The model and method are explained in Sec. II. In Sec. III, the order parameters and magnetic structures are analyzed and the phase diagrams are presented. After the overview of this long section in Sec. III.1, discussions about the characteristics of the model with no, weak, and strong dipolar interactions are given in Secs. III.2, III.3, and III.4, respectively. A variational analysis is performed in Sec. III.5 for understanding the mechanism of the reentrant transition of the helical phase. The summary is given in Sec. IV.

In this study, we consider the system on a square lattice in the $xy$ plane composed of classical Heisenberg spins $\mbox{\boldmath$S$}_{i}$ with $|\mbox{\boldmath$S$}_{i}|=1$ at the position $\mbox{\boldmath$r$}_{i}$ represented by the following Hamiltonian:

with the coupling constants $J$ for the nearest-neighbor ferromagnetic interaction, $D$ for the DM one, $D^{\prime}$ for the dipolar one, $\eta$ for the Ising anisotropy, and $H$ for the magnetic field along the $z$ axis. Vector $\mbox{\boldmath$r$}_{ij}$ is defined as $\mbox{\boldmath$r$}_{ij}\equiv\mbox{\boldmath$r$}_{j}-\mbox{\boldmath$r$}_{i}$, $\mbox{\boldmath$e$}_{x}$ and $\mbox{\boldmath$e$}_{y}$ are unit vectors in the $x$ and $y$ directions, respectively, and site $i+\hat{x}$ ($i+\hat{y}$) is the nearest neighbor one of site $i$ in the positive $x$ ($y$) direction. We consider both positive and negative ($XY$ anisotropy) values for $\eta$. Here we fix the coupling constants of short-range interactions as $J=D=1$.

A previous study of the model for $J=D=1$ and $\eta=0$ without the dipolar interaction by Nishikawa et al. Nishikawa showed that a triangular skyrmion-lattice phase exists and the skyrmion lattice does not have the long-range positional order but has the long-range orientational order. In the present study, to identify the triangular and square skyrmion-lattice phases, we introduce the orientational order parameter which detects the $n$th rotational symmetry ($C_{n}$) as follows. First we detect a domain in which $S_{i}^{z}<S_{\mathrm{th}}\leq 0$ for $H\geq 0$ are fulfilled, and regard this domain as one skyrmion. Then, we define its position as the mean value of the positions of the spins inside the domain:

where ${\cal S}_{k}$ is the $k$th domain. For weak fields, down-spin regions are accidentally generated by thermal fluctuations between skyrmions. These down spins connect the skyrmions momentarily (Fig. 5(b)). To detect skyrmions precisely, we tune the cutoff value $S_{\mathrm{th}}$ depending on $H$. We take $S_{\mathrm{th}}$ smaller for weak fields as $S_{\mathrm{th}}=-0.5$ for $H\leq 0.2$, while $S_{\mathrm{th}}=0$ for $H>0.2$.

The orientational order parameter is given as

where $N_{s}$ denotes the number of skyrmions, $\theta_{jk}$ is the angle between vector $\mbox{\boldmath$x$}_{k}-\mbox{\boldmath$x$}_{j}$ and the $x$ axis, and $\partial j$ is the set of skyrmions adjacent to the $j$th one (Fig. 1). For the triangular and square skyrmion lattices, $n=6$ and $8$ are taken, respectively. We give below the reason why we take $n=8$ instead of $n=4$ to detect the square lattice.

We judge the adjacent skyrmions by the Delaunay triangulation Shewchuk . When the lattice is slightly distorted from the perfect square lattice, about half of the pairs of the next-nearest-neighbor skyrmions are regarded as adjacent ones, and this contribution to $\Psi_{4}$ partially cancels that from the pairs of nearest-neighbor skyrmions. On the other hand, when $\Psi_{8}$ is considered, the above two contributions have the same sign and such cancellation does not take place. We confirmed that $\Psi_{8}$ works as the order parameter of the square skyrmion-lattice phase. For skyrmion-lattice phases, we also apply the local chirality Yi commonly used for the detection of them:

In the model (1), $\chi$ takes a negative value for skyrmion-lattice phases. It should be noted that $\chi=0$ by definition for the configuration with skyrmion and anti-skyrmion pairs at $H=0$ (Fig. 5(b)). Furthermore, $\chi$ detects skyrmions without orientational order, i.e., a liquid state. Therefore, we use $\chi$ to estimate the upper limit of skyrmion-lattice phases at high fields.

In the helical phase, a stripe pattern of up and down spins is formed. To detect the helical phase, we measure another type of orientational order parameters Whitehead defined as

with

Both of these order parameters correspond to the breaking of the symmetry under the $\frac{\pi}{2}$-rotation. The value of $O_{\rm hv}$ is the maximum when [1,0] or [0,1] helical structure is formed, whereas that of $O_{\rm d}$ is the maximum when [1,1] or [1,$-$1] helical structure is formed.

We also measure the uniform in-plane magnetization

with the total number of spins $N$, and it has a nonzero value in the low-temperature phase when $\eta$ is sufficiently small.

In Monte Carlo (MC) simulations of the present study, we use the stochastic cut-off (SCO) $O(N)$ method Sasaki , which reduces the computational cost for dipolar systems. Each Monte Carlo run for a set of $H$ and $\eta$ starts from a uniformly random state at a high temperature. For $\eta\geq 3.0$, the system is cooled from $T=0.34$ to $T=0.16$ at intervals of $T=0.02$ and further cooled to $T=0.1$, and for $\eta<3.0$, it is cooled from $T=0.6$ to $T=0.1$ at intervals of $T=0.05$. At each temperature, 400,000 MC steps are used for the measurement after 100,000 MC steps for the equilibration. The order parameters are estimated by averaging over four independent simulations for the system size $L=84$ ($N=L^{2}$), and the error bar is estimated by $\pm\sigma/\sqrt{4}$, where $\sigma$ is the sample standard deviation. In Figs. 4 and 5, snapshots of spin structures for $L=36$ are shown, which are qualitatively the same as those for $L=84$.

We show a benchmark to compare the SCO method and a naive MC method in Fig. 2. The computational time [s] for 100,000 MCS is plotted as a function of $N$ for the parameters $D^{\prime}=0.6$, $H=1.5$, and $\eta=2.0$. We confirm that the SCO method costs $O(N)$ computational time and is more efficient than the naive MC method which costs $O(N^{2})$ computational time.

We investigated the two-dimensional classical Heisenberg model with the ferromagnetic exchange ($J$), Dzyaloshinskii-Moriya ($D$) and dipolar ($D^{\prime}$) interactions, and the Ising anisotropy ($\eta$) with $J=D=1$. We estimated the order parameters for square and triangular skyrmion-lattice phases, diagonal and vertical helical phases, and in-plane magnetization with the stochastic cut-off $O(N)$ Monte Carlo method. We presented the field ($H$) vs. $\eta$ phase diagrams at a low temperature $T=0.1$ for three typical values of the dipolar interaction, $D^{\prime}=0.0$, $0.3$, and $0.6$ in Fig. 3.

When there is no dipolar interaction ($D^{\prime}=0$), for small $D$ compared to $J$, any skyrmion-lattice phases do not exist at zero and low fields, and the triangular skyrmion-lattice phase appears at high fields. However, for $D$ close to $J$ ($D=J=1$ in the present study), the diagonal square skyrmion-lattice phase exists at zero and low fields between the canted ferromagnetic phase at smaller $\eta$ and the diagonal helical phases at larger $\eta$. At high fields the triangular skyrmion-lattice phase appears at larger $\eta$ than the diagonal square skyrmion-lattice phase. There exists a triple point at which the diagonal square skyrmion-lattice, triangular skyrmion-lattice, and diagonal helical phases coexist.

The dipolar interaction leads to other types of skyrmion lattice and helical phases, which yield four other types of triple points. The effect of the dipolar interaction shifts the phase boundaries to the positive $\eta$ and $H$ directions and stabilizes the ordered phases, whose regions in the $H-\eta$ phase diagram are expanded. The dipolar interaction also increases the skyrmion density and reduces the skyrmion size in the skyrmion-lattice phase and the pitch length (stripe width) in the helical phase.

In both cases of the weak ($D^{\prime}=0.3$) and strong ($D^{\prime}=0.6$) dipolar interactions, any skyrmion-lattice phases do not exist at zero and low fields, which is different from the case without the dipolar interaction ($D^{\prime}=0$). For the weak dipolar interaction ($D^{\prime}=0.3$) at high fields, instead of the diagonal square skyrmion-lattice phase, the upright square skyrmion-lattice phase appears between the canted ferromagnetic phase at smaller $\eta$ and the triangular skyrmion-lattice phase at larger $\eta$. With increasing $\eta$ at zero and low fields, the canted ferromagnetic phase changes to the vertical helical phase, to the diagonal helical phase, and to the vertical helical phase again. This reentrant transition is qualitatively explained by using the variational functions for spins. For the strong dipolar interaction ($D^{\prime}=0.6$), at zero and low fields, only the vertical helical structure exists at larger $\eta$ than the canted ferromagnetic phase.

In the present paper, we focused on the model for mterials with non-centrosymmetric lattice structures, in which the DM interaction plays an important role. Very recently, a square skyrmion-lattice phase was discovered in GdRu₂Si₂ Khanh , which has a centrosymmetric structure. There the DM interaction is not essential, and four-spin interactions are considered to be important. The present paper shows a possible scenario of square skyrmion-lattice phases driven by the DM interaction with or without the dipolar interaction.