Stability and Phase Transition of Skyrmion Crystals Generated by Dzyaloshinskii-Moriya Interaction

TL;DR

This paper investigates the formation and stability of skyrmion crystal phases in a two-dimensional magnetic system influenced by Dzyaloshinskii-Moriya interaction, magnetic field, and temperature, revealing phase transitions and slow relaxation dynamics.

Contribution

It demonstrates the existence of a skyrmion crystal ground state in a 2D Heisenberg model with Dzyaloshinskii-Moriya interaction and studies its finite-temperature stability through Monte Carlo simulations.

Findings

Skyrmion crystal ground state identified in certain (D, H) regions.

Skyrmion phase remains stable at finite temperatures with slow relaxation dynamics.

Transition from skyrmion crystal to paramagnetic state occurs at finite temperature.

Abstract

We generate a crystal of skyrmions in two dimensions using a Heisenberg Hamiltonian including the ferromagnetic interaction J, the Dzyaloshinskii-Moriya interaction D, and an applied magnetic field H. The ground state (GS) is determined by minimizing the interaction energy. We show that the GS is a skyrmion crystal in a region of (D, H). The stability of this skyrmion crystalline phase at finite temperatures is shown by a study of the time-dependence of the order parameter using Monte Carlo simulations. We observe that the relaxation is very slow and follows a stretched exponential law. The skyrmion crystal phase is shown to undergo a transition to the paramagnetic state at a finite temperature.

Equations

$${\mathbf{D}}_{i,j}\cdot {\mathbf{S}}_i\wedge {\mathbf{S}}_j$$

$$\mathcal{H}$$

$${\mathcal{H}}_i=-S_i^x{H}_i^x-S_i^y{H}_i^y-S_i^z{H}_i^{z}y$$

$$H_i^x$$

$$H_i^y$$

$$H_i^{z}z$$

$$M(T)=\frac{1}{N^2(t_a-t_0)}\sum _i|\sum _{t=t_0}^{t_a}{\vec{S}}_i(T,t)\cdot S_i^0(T=0)|$$

$$Q_{+}(T)$$

$$Q_{-}(T)$$

$$M(T,t)=A\exp \left[-(t/\tau {)}^{\alpha }\right]+c,$$

Full text

**PACS numbers: 75.25.-j ; 75.30.Ds ; 75.70.-i **

††preprint: APS/123-QED

I Introduction

Skyrmions have been extensively investigated in condensed matter physics Fert2013 since its theoretical formulation by T. H. R. Skyrme Skyrme in the context of nuclear matter.

There are several mechanisms and interactions leading to the appearance of skyrmions in various kinds of matter. The most popular one is certainly the Dzyaloshinskii-Moriya (DM) interaction which was initially proposed to explain the weak ferromagnetism observed in antiferromagnetic Mn compounds. The phenomenological Landau-Ginzburg model introduced by I. Dzyaloshinskii Dzyaloshinskii was microscopically derived by T. Moriya Moriya . This demonstration shows that the DM interaction comes from the second-order perturbation of the exchange interaction between two spins which is not zero only under some geometrical conditions of non-magnetic atoms found between them. The order of magnitude of DM interaction, $D$, is therefore, perturbation theory obliges, very small. The explicit form of the DM interaction will be given in the next section. However, in many recent papers using the DM interaction, the assumption of small $D$ is not always respected. Therefore, we can think that the demonstration of Moriya Moriya is a special case and the general Hamiltonian may have the same form but different microscopic origin. This is similar the case of the Hubbard model which was initially originated from a second-order perturbation of exchange interaction, but it has been used with liberty for arbitrary ratio $t/U$. This is also the case of the Ising model if we think of it as a limiting case of the Heisenberg model.

The DM interaction has been shown to generate skyrmions in various kinds of crystals. For example, it can generate a crystal of skyrmions in which skyrmions arrange themselves in a periodic structure Rosales ; Buhrandt ; Kwon1 ; Kwon2 . Skyrmions have been shown to exist in crystal liquids Bogdanov2003 ; Leonov ; Ackerman as well as quantum Hall systems. A single skyrmion has also been found. Effects of skyrmions have been investigated in thin films Ezawa ; Yu2 . Artificial skyrmion lattices have been devised for room temperatures Gilbert . Experimental observations of skyrmion lattices have been realized in MnSi in 2009 Muhlbauer ; Bauer and in doped semiconductors in 2010 Yu1 . Needless to say, many potential applications using properties of skyrmions are expected in the years to come. At this stage, it should be noted that skyrmion crystals can also be created by competing exchange interactions without DM interactions Hayami ; Okubo . So, mechanisms for creating skyrmions are multiple.

We note that spin-wave excitations in systems with a DM interaction in the helical phase without skyrmions have been investigated by many authors Diep2017 ; Puszkarski ; Zakeri ; Wang ; Stashkevich ; Moon .

In this paper, we study a skyrmion crystal created by the competition between the nearest-neighbor (NN) ferromagnetic interaction $J$ and the DM interaction of magnitude $D$ under an applied magnetic field $\vec{H}$. We show by Monte Carlo (MC) simulation that the skyrmion crystal is stable at finite temperatures up to a transition temperature $T_{c}$ where the topological structure of each skyrmion and the periodic structure of skyrmions are destroyed.

The paper is organized as follows. Section II is devoted to the description of the model and the method to determine the ground state (GS). It is shown that our model generates a skyrmion crystal with a perfect periodicity at temperature $T=0$. The GS phase diagram in the space $(D,H)$ is presented. Results showing the stability of the skyrmion crystal at finite $T$ obtained from MC simulations are shown in section III. We show in this section that the relaxation of the skyrmions is very slow and follows a stretched exponential law. The stability of the skyrmion phase is destroyed at a phase transition to the paramagnetic state. Concluding remarks are given in section IV.

II Model and Ground State

The DM interaction between two spins $\mathbf{S}_{i}$ and $\mathbf{S}_{j}$ is written as

[TABLE]

where $\mathbf{D}_{i,j}$ is a vector which results from the displacement of non magnetic ions located between $\mathbf{S}_{i}$ and $\mathbf{S}_{j}$, for example in Mn-O-Mn bonds in the historical papers Dzyaloshinskii ; Moriya . The direction of $\mathbf{D}_{i,j}$ depends on the symmetry of the displacement Moriya . For two spins, the DM interaction is antisymmetric with respect to the inversion symmetry.

Theoretical and experimental investigations on the effect of the DM interaction in various materials have been extensively carried out in the context of weak ferromagnetism observed in perovskite compounds (see references cited in Refs. Sergienko, ; Ederer, ). As said in the Introduction, the interest in the DM interaction goes beyond the weak ferromagnetism and the model is used beyond the perturbation limit. It has been shown that the DM interaction is at the origin of topological skyrmions Muhlbauer ; Yu2 ; Yu1 ; Maleyev ; Lin ; Bogdanov ; Rossler ; Seki ; Adams ; Heurich ; Wessely ; Jonietz and new kinds of magnetic domain walls Heide ; Rohart . The increasing interest in skyrmions results from the fact that skyrmions may play an important role in the electronic transport which is at the heart of technological application devices Fert2013 .

In this paper, we consider for simplicity the two-dimensional case where the spins are on a square lattice in the $xy$ plane. We are interested in the stability of the skyrmion crystal generated in a system of spins interacting with each other via a DM interaction and a symmetric isotropic Heisenberg exchange interaction in an applied field perpendicular to the $xy$ plane. All interactions are limited to NN. The full Hamiltonian is given by

[TABLE]

where the DM interaction and the exchange interaction are taken between NN on both $x$ and $y$ directions. Rewriting it it in a convenient form, we have

[TABLE]

For the $i$-th spin, one has

[TABLE]

where the local-field components are given by

[TABLE]

To determine the ground state (GS), we minimize the energy of each spin, one after another. This can be numerically achieved as the following. At each spin, we calculate its local-field components acting on it from its NN using the above equations. Next we align the spin in its local field, i. e. taking $S_{i}^{x}=H_{i}^{x}/\sqrt{H_{i}^{x}**2+H_{i}^{y}**2+H_{i}^{z}**2}$ etc. The denominator is the modulus of the local field. In doing so, the spin modulus is normalized to be 1. As seen from Eq. (II), the energy of the spin $\vec{S}_{i}$ is minimum. We take another spin and repeat the same procedure until all spins are visited. This achieves one iteration. We have to do a sufficient number of iterations until the system energy converges. For the skyrmion case, it takes about one thousand iterations to have the fifth-digit convergence.

An example of GS are shown in Fig. 1: a crystal of skyrmions is seen using $D=1$ and $H=0.5$ (in unit of $J=1$).

In Fig. 2a we show a GS at $H=0$ where domains of long and round islands of up spins separated by labyrinths of down spins are mixed. When $H$ is increased, vortices begin to appear. The GS is a mixing of long islands of up spins and vortices as seen in Fig. 2b obtained with $D=1$ and $H=0.25$. This phase can be called ”labyrinth phase” or ”stripe phase”.

It is interesting to note that skyrmion crystals with texture similar to those shown in Figs. 1 and 2 have been experimentally observed in various materials Yu2 ; Gilbert ; Muhlbauer ; Yu3 , but the most similar skyrmion crystal was observed in two-dimensional Fe(0.5)Co(0.5)Si by Yu et al. using Lorentz transmission electron microscopy Yu1 .

We have performed the GS calculation taking many values in the plane $(D,H)$. The phase diagram is established in Fig. 3. Above the blue line is the field-induced ferromagnetic phase. Below the red line is the labyrinth phase with a mixing of skyrmions and rectangular domains. The skyrmion crystal phase is found in a narrow region between these two lines.

In the following section, we are interested in the stability of the skyrmion crystal phase as the temperature ($T$) is increased from zero.

III Stability of Skyrmion Crystal at Finite Temperatures

In this section, we show results obtained from MC simulations on a sheet of square lattice of size $N\times N$ with periodic boundary conditions. The first step is to determine the GS spin configuration by minimizing the spin energy by iteration as described above. Using this GS configuration, we heat the system from $T=0$ to a temperature $T$ during an equilibrating time $t_{0}$ before averaging physical quantities over the next $10^{6}$ MC steps per spin. The time $t_{0}$ is the ”waiting time” during which the system relaxes before we perform averaging during the next $t_{a}$.

The definition of an order parameter for a skyrmion crystal is not obvious. Taking advantage of the fact that we know the GS, we define the order parameter as the projection of an actual spin configuration at a given $T$ on its GS and we take the time average. This order parameter is thus defined as

[TABLE]

where $\vec{S}_{i}(T,t)$ is the $i$-th spin at the time $t$, at temperature $T$, and $\vec{S}_{i}(T=0)$ is its state in the GS. The order parameter $M(T)$ is close to 1 at very low $T$ where each spin is only weakly deviated from its state in the GS. $M(T)$ is zero when every spin strongly fluctuates in the paramagnetic state. The above definition of $M(T)$ is similar to the Edward-Anderson order parameter used to measure the degree of freezing in spin glasses Mezard : we follow each spin with time evolving and take the spatial average at the end.

We show in Fig. 4 the order parameter $M$ versus $T$ (red data points) as well as the average $z$ spin component (blue data points) calculated by the projection procedure for the total time $t=10^{5}+10^{6}$ MC steps per spin. As seen, both two curves indicate a phase transition at $T_{c}\simeq 0.26J/k_{B}$. The fact that $M$ does not vanish above $T_{c}$ is due to the effect of the applied field. It should be said that each skyrmion has a center with spins of negative $z$ components (the most negative at the center), the spins turn progressively to positive $z$ components while going away from the center.

We can also define another order parameter: since the field acts on the $z$ direction, in the GS and in the skyrmion crystalline phase we have both positive and negative $S_{z}$. In the paramagnetic state, the negative $S_{z}$ will turn to the field direction. We define thus the following parameters using the $z$ spin-components

[TABLE]

Figure 5 shows $Q_{+}$ and $Q_{-}$ versus $T$. As seen, at the transition $Q_{+}$ undergoes a change of curvature and $Q_{-}$ becomes zero. All spins have positive $S_{z}$ after the transition due to spin reversal by the field.

The results obtained at the end of the simulation may depend on the overall time $t=t_{0}+t_{a}$. In simple systems, the choices of $t_{0}$ and $t_{a}$ can be guided by testing the time-dependence of physical quantities, and the values of $t_{0}$ and $t_{a}$ are chosen when physical quantities do not depend on these run times. However, in disordered systems such as spin glasses and in complicated systems such as frustrated systems, the relaxation time is very long and out of the reach of simulation time. In such cases, we have to recourse to some scaling relations in order to deduce the values of physical quantities at equilibrium Phillips ; Ogielski . We show below how to obtain the value of an order parameter at the infinite time.

In order to detect the dependence of $M(T)$ on the total MC time $t=t_{0}+t_{a}$, we calculate the average of $M(T)$ over $10^{6}$ MC steps per spin, after a waiting time $t_{0}$ as said above. We record the values of $M(T)$ in different runs with $t_{0}$ varying from $10^{4}$ to $10^{6}$ MC steps per spin. We plot these results as a function of different total time $t$ in Fig. 6 for three temperatures. As said above, to find the value extrapolated at the infinite time, we use the stretched exponential relaxation defined by

[TABLE]

where $t$ is the total simulation time, $\alpha$ is the stretched exponent, $A$ a temperature-dependent constant, and $\tau$ the relaxation time. Note that this definition, without the constant $c$, has been used by many previous authors in the context of spin glasses DeDominicis ; Almeida ; Campbell2011 ; Suzuki ; Malinowski2011 . We have introduced $c$ which is the infinite-time limit of $M(T)$. We have taken $t$ from $10^{4}$ to $10^{6}$ MC steps per spin in the simulation. At the infinite-time limit, $c$ is zero for $T\gg T_{c}$, and $c\neq 0$ for $T<T_{c}$. Figure 6 shows $M(T,t)$ as a function of time $t$ in unit of $10^{3}$ MC steps per spin, for three temperatures $T=0.01$, 0.094 and 0.17. As seen, the fit with Eq. (7) presented by the continuous line is very good for the whole range of $t$.

Several remarks are in order:

(i) the precision of all parameters are between 1% to 5% depending on the parameter,

(ii) the value of $\alpha$ can vary a little bit according to the choice and the precision of the other parameters in the fitting but this variation is within a very small window of values around the value given above. For example, at $T=0.17$, $\alpha$ can only be in the interval $[0.8\pm 0.02]$. The value of $\alpha$ can vary with temperature as seen here: at low $T$, $\alpha=0.6$, and at a higher $T$, we have $\alpha=0.8$. This variation has been seen in other systems, in particular in spin glasses Ngo2014 .

(iii) the relaxation time, within statistical errors, is approximatively constant at low $T$, but it increases rapidly when $T$ tends to $T_{c}$ as seen in the value of $\tau$ at $T=0.17$. This increase is a consequence of the so-called ”critical slowing-down” when the system enters the critical region.

Let us show $M(T)$ as a function of $T$ in Fig. 7 using the results of different run times from $t=10^{4}+10^{6}$ MC steps per spin to $t=10^{6}+10^{6}$. The values of at infinite time for each $T$ deduced from Eq. (7) is also shown. We see that while the total time $10^{5}+10^{6}$ MC steps per spin is sufficient at low $T$, it is not enough at higher $T$. That was the reason why we should use Eq. (7) to find the value of $M(T)$ at infinite time to be sure that the skyrmion crystal is stable at finite temperatures.

We have studied finite-size effects on the phase transition at $T_{c}$ and we have seen that from $N=800$, all curves coincide: there is no observable finite size effects for $N\geq 800$. The present 2D skyrmion lattice phase therefore remains stable for the infinite dimension, unlike the ferromagnetic Heisenberg model in 2D Mermin . Note that the nature of the ordering at low $T$ and the phase transition observed here, in spite of the vortex nature of skyrmions, are not those of the Kosterlitz-Thouless vortex mechanism observed in ferromagnetic XY spin systems KT since our ”vortices” are stable without changing their position up to the transition temperature.

IV Concluding Remarks

In this paper, we have shown that the competition between a ferromagnetic interaction $J$ and a Dzyaloshinskii-Moriya interaction $D$ under an applied magnetic field $H$ in two dimensions generate a skyrmion crystal in a region of the phase space $(D,H)$. The spin model is the classical Heisenberg model. We have numerically determined the ground state by minimizing the local energy, spin by spin, using an iteration procedure. The skyrmion lattice is then heated to a finite temperature by the use of Monte Carlo simulations. We have shown that the skyrmion lattice is stable up to a finite temperature $T_{c}$ beyond which the system becomes disordered. We have also shown that the relaxation follows a stretched exponential law. We believe that such a stability can be experimentally observed in real systems.

Perhaps, our model is the simplest model able to generate a skyrmion crystal without the need to include more complicated interactions such as long-range dipolar interactions, easy and uniaxial anisotropies Kwon2 ; Lin ; Yu3 . Note that other authors have shown that skyrmion crystals can also be generated only with frustrated short-range interactions without even the DM interaction Hayami ; Okubo . So, we have to keep in mind that there are many interaction mechanisms of different nature which can generate skyrmion crystals. The origin of a skyrmion crystal determines its structure, its stability, its dynamics and in short its properties (see a discussion on this point in Yu3 ). Experiments can be used therefore to determine the interaction mechanism which is at the origin of the formation of a skyrmion crystal.

Acknowledgements.

SEH acknowledges a financial support from Agence Universitaire de la Francophonie (AUF).