 Skyrmion Crystal: Phase Transition

In this section, we show results obtained from MC simulations on a sheet of square lattice of size N x 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 to before averaging physical quantities over the next 106 MC steps per spin. The time to is the "waiting time" during which the system relaxes before we perform averaging during the next

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 where S,(7 t) is the /-th spin at the time t, at temperature T, and Si[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 : We follow each spin during the time evolution and take the spatial average at the end.

We show in Fig. 14.5 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 = 105 + 106 MC steps per spin. As seen, both two curves indicate a phase transition at Tc ~ 0.26] /kB. The fact that M does not vanish above Tc is due to the effect of the applied field. Note 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 S2. In the paramagnetic state, Figure 14.5 Red circles: Order parameter defined in Eq. (14.8) versus T, for H = 0.5 and N = 1800, averaged during f;, = 10s MC steps per spin after an equilibrating time to = 105 MC steps. Blue crosses: the projection of the Sz on 5“ of the ground state as defined in Eq. (14.8) but for the z components only. See text for comments.

the negative Sz will turn to the field direction. We define thus the following parameters using the z spin-components Figure 14.6 shows Q+ and versus T. As seen, at the transition Q+ undergoes a change of curvature and becomes zero. All spins have positive Sz 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 = to + tg. In simple systems, the choices of to and tg can be guided by testing the time-dependence of physical quantities, and the values of to and t„ 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 [261, 273]. We show below how to obtain the value of an order parameter at the infinite time. Figure 14.6 Order parameters defined in Eqs. (14.9)—(14.10) versus T, for H = 0.5 and N = 800, to = 10s, t„ = 106. Red circles: total positive z component ((?+(?'))• Green circles: total negative z component Blue

circles: total z component. See text for comments.

Stability of Skyrmion Crystal at Finite Temperatures

In order to detect the dependence of M{T) on the total MC time t = fo + ta, we calculate the average of M[T) over 106 MC steps per spin, after a waiting time to as said above. We record the values of М(Г) in different runs with fo varying from 104 to 106 MC steps per spin. We plot these results as a function of different total time t in Fig. 14.7 for three temperatures. As said above, to find the value extrapolated at the infinite time, we use the stretched exponential relaxation defined by where t is the total simulation time, a is the stretched exponent, A a temperature-dependent constant, and r the relaxation time. Note that this definition, without the constant c, has been used by many previous authors in the context of spin glasses [9, 32, 94, 225, 335]. We have introduced c which is the infinite-time limit of M[T). We have taken t from 104 to 106 MC steps per spin in the simulation. At the infinite-time limit, c is zero for T » Tc, and с Ф 0 for T < Tc. Figure 14.7 shows M[T, t) as a function of time t in unit of 103 MC Figure 14.7 The order parameter M defined by Eq. (14.8) versus MC time t in unit of 1000 MC steps per spin, for H = 0.5 and N = 800 (a) T = 0.01, values of fitting parameters a = 0.6, A = 0.008 ± 0.00001, r = (364 ± 19)103, c = 0.984505 ± 0.00012; (b) T = 0.094, a = 0.6, A = 0.0653 ± 0.0013, т = (277 ± З6)103, c = 0.822 ± 0.018); (с) T = 0.17, a = 0.8, A = 0.31 ± 0.01, r = 891 ± 100)103, c = 0.43 ± 0.017.

steps per spin, for three temperatures T = 0.01, 0.094 and 0.17. As seen, the fit with Eq. (14.11) 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 a 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, a can only be in the interval [0.8 ± 0.02]. The value of a can vary with temperature as seen here: At low T, a = 0.6, and at a higher T, we have a = 0.8. This variation has been seen in other systems, in particular in spin glasses . Figure 14.8 The order parameter М(Г) versus T for several waiting times t, for H = 0.5 and N = 800: curves 1, 2, 3 and 4 are, respectively, for t = 10s, 2 x 10s, 106, oc (by fitting with Eq. (14.11).

• (iii) The stretched exponential behavior is believed to be due to the non-uniform spin configuration since collinear spin states (ferromagnetic an antiferromagnetic orderings) follow a non- stretched exponential law . The stretched exponential law was found in disordered phases such as spin glasses [9, 32, 94, 225, 335]. Of course if D is small, the configuration tends to a ferromagnetic one, the relaxation tends to a non-stretched law.
• (iv) The relaxation time, within statistical errors, is approximately constant at low T, but it increases rapidly when T tends to Tc as seen in the value of r 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. 14.8 using the results of different run times from t = 104 + 106 MC steps per spin to t = 106 -(- 106. The values at infinite time for each T deduced from Eq.

(14.11) is also shown. We see that while the total time 10s + 106 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. (14.11) to find the value of M(T) at infinite time to be sure that the skyrmion ciystal is stable at finite temperatures.

We have studied finite-size effects on the phase transition at Tc and we have seen that from N = 800, all curves coincide: There is no observable finite-size effects for N > 800. The present 2D skyrmion lattice phase, therefore, remains stable for the infinite dimension, unlike the ferromagnetic Heisenberg model in 2D . Note that the nature of the ordering at low T and the phase transition observed here: In spite of the vortex nature of skyrmions, the transition is not due to the Kosterlitz-Thouless vortex mechanism observed in ferromagnetic XY spin systems [191, 192] since our "vortices" are stable without changing their position up to the transition temperature.

Skyrmion Crystal: Effect of Lattice Elasticity

In this section, we would like to see whether the skyrmion phase observed in the preceding sections is stable under large lattice deformations . For this purpose, it is interesting to assume the lattice degrees of freedom as a dynamical variable. Here we should note that "dynamical" variable refers to the variable that is integrated in the partition function in the context of statistical mechanics. Therefore, a fluctuating lattice is considered to be a direct tool for this purpose, where the lattice fluctuation is allowed not only into the in-plane direction but also into the out-of-plane direction [107, 127, 145, 174, 246, 260, 279, 360]. Skyrmions on fluctuating lattice are interesting also from the viewpoint of liquid crystal skyrmions [39, 207].

It is widely accepted that skyrmion configurations, which are observed on materials such as FeGe, MnSi, or Cu20Se03, are promising candidates for future computer memories because of their dramatic energy-saving property in electric transportation. Due to the topological nature of the spin configuration, its stability against external stimuli such as thermal fluctuations is one of its most interesting properties [37, 40, 99, 240, 296, 325]. For these reasons, many experimental and theoretical studies have been conducted by assuming competing interactions, such as ferromagnetic (FM) interaction, Dzyaloshinskii-Moriya [DM], and Zeeman energy under constant magnetic fields [51, 96, 163, 242,

243, 379]. In particular, skyrmions in 2D systems have attracted a great deal of attention due to their potential applications in memory devices [20, 38, 122, 135]. However, the creation/annihilation of skyrmions and the stability of the skyrmion phase are yet to be studied in the context of such applications.

Along this research direction, several investigations have been conducted on the responses to mechanical stresses, which have a non-trivial effect on the spin configuration [245, 290]. In Ref. , Shi and Wang numerically find that skyrmions are stabilized due to the magnetoelastic coupling even without external magnetic field by introducing an explicit magnetoelastic interaction energy. This confirms that skyrmions may be generated by various mechanisms. We come back to this point later in the next section. Nii et al. recently reported that MnSi is an anisotropic elastic object , and they also reported the experimental results of a stress-induced skyrmion-to-conical phase transition in MnSi, where the applied uniaxial stresses create/annihilate skyrmions under suitable external magnetic fields . It is also reported by Seki et al. that skyrmions in Си2ОБеОз can be stabilized by a uniaxial tensile strain [313, 314], and Levatic et al. reported that mechanical pressure significantly increases the skyrmion pocket in Cu20Se03  and in MnSi . From these experiments, it was clarified that uniaxial mechanical stimuli stabilize/destabilize the skyrmion phase, or, in other words, a small lattice deformation enhances the stability/instability of the skyrmion phase. This indicates that the interplay between two different degrees of freedom spin and lattice, which can be rephrased by spin and orbit, plays an important role in mechanical responses, which has been extensively studied via DM interaction energy. Chen et al. numerically studied skyrmions in MnSi by using a 3D model and suggested that anisotropies in DM and FM interactions play an important role in their stability under uniaxial stress in a wide range of temperatures including low temperature region . These anisotropies are naturally expected when crystalline lattices are considerably deformed. Thus, it is interesting to study the skyrmion stability under lattice deformations at finite temperatures.

The DM interaction on the continuous 2D plane depends on both spins and coordinates, and hence, we are able to calculate mechanical properties such as surface tension under the presence of skyrmions if the coordinate, or equivalently the surface, is treated as a dynamical variable [107, 127, 145, 174, 246, 260, 279, 360]. Moreover, the DM interaction is invariant under an arbitrary rotation around the axis perpendicular to the surface. However, its discrete version, which will be described in the next section, has only a discrete symmetry on 2D rigid lattices because the rigid lattices have only discrete symmetries. In contrast, the fluctuating lattices, where the edge direction of the lattice becomes almost random, are not regarded as a rigid lattice by any rotation. For this reason, the discrete DM interaction is expected to be influenced by the surface fluctuations.

Therefore, it is worthwhile studying skyrmions to see the interplay between spins and surface elasticity on fluctuating surfaces. In this section, we calculate the surface tension r by Monte Carlo (MC) simulations to study the influences of skyrmion configurations on r. Three different spin configurations are expected on the surface: stripe, skyrmion, and FM phases, as in the case of rigid lattices [103, 196, 371, 372]. Phase diagrams are also obtained in the space temperature versus external magnetic field. In this calculation, mechanical stress can be applied to the elastic lattices only in the form of strain, which is given by fixing the boundary to a suitable size or shape. Such an effective stress becomes isotropic in general, because the lattice structure changes relatively freely due to the fact that the vertex position is a dynamical variable. Nevertheless, the uniaxial stress condition is simulated by modifying the boundary lengths of the rectangular lattice, for example, from its isotropic lengths (L1; L2) to (Lx/f, £L2), (£ Ф 1). As a result of this geometry modification, we have to check whether the stripe domain axis aligns along the uniaxial stress direction. We also examine whether this uniaxial stress condition influences the stability of skyrmion phase under an external magnetic field.

Triangulated Lattices and Skyrmion Model

To study skyrmions on the fluctuating lattices, we have to combine the skyrmion model and the elastic surface model [107, 127, 174, 246, 260, 360]. If the two Hamiltonians are merged, the total

Hamiltonian becomes slightly lengthy, but the definition itself is straightforward. Indeed, the Hamiltonian S is given by a linear combination of several terms such as The first and second terms are the discrete Hamiltonians of the surface model of Helfrich and Polyakov [145, 279]. These Si and S2 are called Gaussian bond potential and bending energy, respectively, and j — Г/1) is the length of bond ij, and n, and n, in S2 are the unit normal vectors of the triangles that share the bond ij. The parameter к is the bending rigidity of the surface [107, 127, 174, 246,360].

The third term ASF is the energy for describing the FM interaction between two neighboring spins with the interaction strength Л. The fourth term SSDM describes the DM interaction with the interaction coefficient 8, where е,Д= (r, -r,)/|r, — r,|) is the unit tangential vector from r, to r;. This vector e,y corresponds to 3r in the continuous SDm which varies on the fluctuating lattices. It is also possible to use £,Д= r; — r,) for 3r; however, here we choose to use the unit length vector e,y for 3r. The continuous expression о x d0a can also be reduced to a у x oy in the discrete Sdm, because the differential dao is replaced by the difference cry —ay. Thus, the term SDM plays a role in the interaction between о and r, the spins and the surface shape. The final term Sg is the energy for the external magnetic field B = (0, 0, В).

We note that the effects of surface fluctuations on spins are taken into consideration only via the DM interaction energy SDm- The tangential vector e,y in SDm is deformed by the surface fluctuations as mentioned above. On the other hand, the FM interaction can also reflect surface fluctuations, because the magnetoelastic coupling is not negligible on the fluctuating surfaces due to the spin-orbit coupling . However, in the case of skyrmion deformations in a strained crystal, it is expected that the deformation of DM interaction plays an essential role, which is precisely checked in Ref. . So, in this section, we simply use the ordinary undeformed Sp in Eq. (14.12). In other words, we neglect the effect of magnetoelastic coupling in the FM interaction and we study only the effects of surface elasticity on spins via the DM interaction deformation.

The partition function is given by where p = l/T,кв = 1). Xa. The symbol p is the inverse temperature, and the simulation unit is defined by the relation kB = 1, and, hence, the temperature T is the reduced temperature. The symbols ./П/dr,X1WBD, /П,А& denote the three-dimensional and two- dimensional multiple integrations of the vertices inside and on the boundaiy, respectively, where Nbo is the total number of vertices on the boundary. The symbol С(Л,,) denotes a constraint on the integration of the boundary vertices so that the projected area remains fixed to Ap. This constraint is a function of r1( • ■ • , However, this function С(г1( • ■ ■ , rWl)D; Л,,) is too complex to write down with Ap, and for this reason we use Cf/l,,) for the short notation .

Results

We studied a 2D skyrmion system on fluctuating surfaces with periodic boundary conditions using a Monte Carlo simulation technique. Note that in this model, not only spins but also lattice vertices are integrated into the partition function of the model. From the obtained results, we conclude that skyrmions are stable even on fluctuating surfaces, contrary to initial expectations. These results have been recently published in Ref. . The reader is referred to this paper for details.

The Dzyaloshinskii-Moriya (DM) interaction is assumed together with the ferromagnetic interaction and Zeeman energy under an external magnetic field. We assume that the lattice fluctuations are only reflected in the DM interaction, and the magnetoelastic coupling in the ferromagnetic interaction is not taken into consideration. Note that in our model the competition of Sf and SDM suffices to create a skyrmion crystal. To keep the model as simple as possible in order to analyze the elasticity effect, we did not include long-range magnetic dipole-dipole interaction and anisotropy in the model. These ingredients will yield another kind of skyrmions which give rise to another behavior under elasticity and uniaxial strains. In addition to those skyrmion Hamiltonians, the Helfrich-Polyakov Hamiltonian for membranes is assumed. A surface governed by these Hamiltonians spans a rectangular frame of a fixed projected area, and the surface (or frame) tension is calculated. In this lattice calculation, the effect of stress is introduced by strains, which are imposed on the lattices by fixing the edge lengths of the boundary constant with periodic boundary conditions. In this sense, this effective stress is automatically introduced by the boundary frame. No mechanical stress can be imposed on the lattice as an input, but it is obtained as an output of the simulations.

We obtain numerical results that are consistent with the experimentally observed fact that a lateral pressure imposed on a 3D material creates skyrmions under a longitudinal magnetic field. Our result indicates that skyrmions are stable on 2D fluctuating surfaces, where there is no crystalline structure in contrast to the rigid lattices including real crystalline materials.

We further examine how uniaxial strains influence spin configurations. To impose a uniaxial strain on the lattice, we assume an anisotropic condition for the edge length of the rectangular boundaiy. As a result of this calculation, we find that the effective pressure or stress increment makes the stripe domain align to the stress direction under the same magnetic field and temperature as in the case without the uniaxial strain. As mentioned in Ref. , the responses of stripe domains under uniaxial strains depend on the physical mechanism which generate skyrmions. Our model is not consistent with the reported experimental results of Ref. , but the fact that our stripes are parallel to uniaxial strains may be found in other materials and may have interesting transport applications. Of course, a more general magnetoelastic coupling included not only in SDM but also in Sf should be considered in our future work [322, 355]. In contrast, almost no difference is observed in the skyrmion phase, and this also confirms the stability of the skyrmion phase on fluctuating surfaces.

We would like to comment on the studies of skyrmion stability on fluctuating lattices from the viewpoint of Thiele equation . In Ref. , an interaction between skyrmion pairs on two separated surfaces is studied by Thiele equation. The authors in Ref.  assume that the bilayer lattice is flat, and hence, the effects of surface fluctuations are not taken into consideration. Therefore, to clarify the influence of the surface fluctuations on the skyrmion stability, it is interesting to study their model on a bilayer which consists of fluctuating and triangulated lattices. One possible technique is simply to assume that the fluctuating lattices are frozen and only spin variables are time dependent.

Finally, we also comment on a possibility of experimental studies. Skyrmions can only be observed on the materials such as FeGe, MnSi, or Cu20Se03 as mentioned in Ref. , and we have no possibility at present to observe and study skyrmions on elastic surfaces such as ferroelectric polymers under the condition similar to those described in this. However, it seems possible to measure the effects of local stresses even on the above-mentioned materials FeGe, etc. by using the technique developed for individual manipulation of superconducting magnetic vortex . It is interesting to study the stability of skyrmions under local stress or strains. This may be the first step toward the understanding of the role of surface fluctuations in the skyrmion stability.

Conclusion

In this chapter, 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 Tc 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 ciystal without the need to include more complicated interactions such as long-range dipolar interactions, easy and uniaxial anisotropies [196, 213, 372]. Note that other authors have shown that skyrmion crystals can also be generated only with frustrated short-range interactions without even the DM interaction [141, 263]. 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 ). Experiments can be used, therefore, to determine the interaction mechanism which is at the origin of the formation of a skyrmion ciystal in the material under investigation.

We have also studied the effect of the lattice elasticity on the skyrmion stability in two dimensions. The Helfrich-Polyakov Hamiltonian for membranes has been assumed in addition to other skyrmion terms in the previous sections. A surface governed by these terms spans a rectangular frame of a fixed projected area, and the surface (or frame) tension is calculated. We found that the fluctuating surface does not destroy the skyrmion structure. The effect of stress has been calculated. We obtained numerical results consistent with the experimentally observed fact. Our result indicates that skyrmions are stable on 2D fluctuating surfaces where there is no crystalline structure. We further examined how uniaxial strains influence the spin configuration. We find that the effective pressure or stress increment makes the stripe domain align to the stress direction under the same magnetic field and temperature as in the case without the uniaxial strain. Details are given in Ref. .