Mechanism of Skyrmion Attraction in Chiral Magnets near the Ordering Temperatures

Isolated chiral skyrmions are investigated within the phenomenological Dzyaloshinskii model near the ordering temperatures of quasi-two-dimensional chiral magnets with Cnv symmetry and three-dimensional cubic helimagnets. In the former case, isolated skyrmions (IS) perfectly blend into the homogeneously magnetized state. The interaction between these particle-like states, being repulsive in a broad low-temperature (LT) range, is found to switch into attraction at high temperatures (HT). This leads to a remarkable confinement effect: near the ordering temperature, skyrmions exist only as bound states. This is a consequence of the coupling between the magnitude and the angular part of the order parameter, which becomes pronounced at HT. The nascent conical state in bulk cubic helimagnets, on the contrary, is shown to shape skyrmion internal structure and to substantiate the attraction between them. Although the attracting skyrmion interaction in this case is explained by the reduction of the total pair energy due to the overlap of skyrmion shells, which are circular domain boundaries with the positive energy density formed with respect to the surrounding host phase, additional magnetization “ripples” at the skyrmion outskirt may lead to attraction also at larger length scales. The present work provides fundamental insights into the mechanism for complex mesophase formation near the ordering temperatures and constitutes a first step to explain the phenomenon of multifarious precursor effects in that temperature region.

The stabilization mechanism of chiral ISs relies on the Dzyaloshinskii-Moriya interaction (DMI) [17][18][19][20], which is generally expressed as the energy terms with the first-order derivatives of the magnetization M with respect to the spatial coordinates, the so called Lifshitz invariants (LI) [21]: These LIs arise in certain combinations depending on the crystal symmetry of an underlying chiral magnet [1]. Moreover, they lead to a unique internal structure of ISs by defining their vorticity and helicity. For cubic helimagnets belonging to 23 (T) (as MnSi [4], FeGe [5], and other B20 compounds) and 432 (O) crystallographic classes, ISs are truly three-dimensional (3D) states. Their two-dimensional (2D) pattern on the plane xy is known to be modified by the LI along z axis, which becomes evident, e.g., in additional surface twists near surfaces of thin-film samples [22][23][24]. Moreover, skyrmions can be embraced by the conical phase with the wave vector along z, which shapes IS internal structure and stipulates attractive nature of inter-skyrmion potentials [25,26]. In, e.g., the C nv symmetry class, however, no Lifshitz invariants are present along the high-symmetry z axis [1]. Thus, ISs retain their 2D structure in both bulk and nanofabricated magnets [27,28].
In frustrated magnets, the limitations of the Derrick-Hobart theorem are overcome by the competing ferromagnetic and antiferromagnetic exchange interactions between spins, e.g., on a triangular lattice, which generate higher-order derivative terms in the fundamental continuum form of the magnetic free energy density [10][11][12][13]. In these frustrated magnets, the spin direction in the skyrmion centre is still opposite to that of the applied magnetic field. However, the topological charge of skyrmions can be both +1 (skyrmions) and −1 (anti-skyrmions) owing to the arbitrary sign of vorticity. In addition, the helicity angle χ is arbitrary and represents a zero mode. χ = 0, π corresponds to skyrmions carrying a monopole moment, A ∝ ∑ j x j · S j , whereas χ = ±π/2 distinguishes skyrmions with a toroidal moment, [29]. These degrees of freedom, however, do not exist for skyrmions in chiral magnets. As a matter of fact, DMI selects only one type of IS according to the symmetry arguments, and other localized solitonic solutions describe metastable or instable particles with higher energy.
In contrast to the fixed and rigid non-collinear spin pattern of skyrmions in chiral magnets, the skyrmions in frustrated magnets, thus, have internal degrees of freedom. It is instructive to compare the properties, shape and mutual interactions of skyrmions in frustrated magnetics to those in chiral helimagnets with DMIs, regarding these internal degrees of freedom. Interestingly, the toroidal and monopole moment densities within these skyrmions exhibit periodic sign reversals along the radial directions with growing distance r from the IS center [11,12]. In these patterns, the spins at the skyrmion outskirt undergo fan oscillations with decaying amplitude, which also give rise to a number of minima in the skyrmion-skyrmion interaction potentials U(r 12 ). On the one hand, mutual skyrmion interaction depends on the topological charges of skyrmions. Two skyrmions or two antiskyrmions attract each other at distances of the order of the skyrmion diameter, the maximal reduction of energy being ∼10% of the skyrmion eigen-energy, which corresponds to the "deepest" minimum of U(r 12 ). Hence frustrated ISs tend to form clusters with a shortrange crystal order, resembling the clustering of vortices in 1.5-type superconductors [30]. On the other hand, interactions between frustrated skyrmions depend on their helicities. This effect couples helicity dynamics to the translational motion of skyrmions during, e.g., current-driven dynamics [11,12]. The oscillating interactions between skyrmions in the frustrated magnets are in sharp contrast with skyrmion-skyrmion interaction in noncentrosymmetric systems, which do not show fan-like oscillations and repel each other at all distances [31]. Hence, one would not expect chiral helimagnets to display similar attractive and complex skyrmion-skyrmion interactions as in frustrated systems at all.
In the present paper, chiral isolated skyrmions are analysed close to the ordering temperature using the phenomenological Dzyaloshinskii model supplemented with the basic Landau expansion for the homogeneous part of the free-energy. This study is motivated by the well-known fact that the longitudinal magnitude of the order-parameter, here the ferromagnetic magnetization M, and the non-collinear twisting are not independent. The non-collinear twisted pattern of skyrmions, or any other localized spin-pattern like kinks, becomes inhomogeneous also with respect to the local magnetization M near the ordering temperature, where the order parameter is soft both in transversal and in longitudinal direction. The inherent frustration by the LIs therefore does not only twist the orientable order-parameter considered as a field distribution of fixed length vector, but it also modulates its magnitude. In a short manner, we may say that for the order-parameter fields in Dzyaloshinskii models the Goldstone modes, i.e., transversal spin distortions, and the longitudinal Higgs mode are intertwined because of the gauge-background encoded in Equation (1). This behavior leads to interesting consequences. The inhomogeneous nature of periodic skyrmion lattice solutions has been described earlier [32,33]. Whereas in a broad temperature range the magnetization modulus |M| remains constant and is independent on the value of the applied magnetic field, its longitudinal stiffness decreases near the ordering temperature, and spatial longitudinal modulations of the magnetization become a sizable effect. We show that isolated skyrmions develop haloes of damped and oscillatory spin twistings in radial directions around their rigid core, accompanied by oscillations of M. This modification of the skyrmion pattern creates attractive inter-skyrmion pair-potentials, and achieves also a longer range of the interaction. The crossover from the repulsive to attractive potential is a consequence of the coupling between the magnitude and the angular part of the magnetization M and occurs at the so-called confinement temperature.
First, from numerical investigations on 2D models of isotropic chiral ferromagnets, we investigate the properties of chiral ISs and find complex combinations of rotational and longitudinal degrees of freedom. In particular, we focus on the internal structure of ISs undergoing a collapse process at high fields. Second, we extend the paradigm to 3D-ISs surrounded by the conical phase in bulk cubic helimagnets and confirm the same attracting interaction due to the oscillatory asymptotic behaviour. Our results identify a fundamental mechanism in chiral magnets and similar systems described by Dzyaloshinskii models, in that particle-like solutions display complex oscillatory transformations of their patterns near phase transition, as their orientational and longitudinal ordering modes are coupled. For chiral helimagnets this point to the long-range entangled character of the skyrmionskyrmion interaction which underlies precursor phenomena in chiral magnets near the ordering temperatures, as observed in experiment [34][35][36].

Phenomenological Theory
Within the phenomenological theory introduced by Dzyaloshinskii [17][18][19]37] the magnetic energy density of a non-centrosymmetric ferromagnet with spatially dependent magnetization M can be written as where m is the unit magnetization vector; A > 0 and D are coefficients of exchange and Dzyaloshinskii-Moriya interactions; H is a magnetic field applied along the z axis; x i are the Cartesian components of the spatial variable. The functional (2) contains only the basic interactions essential to stabilize skyrmionic states in noncentrosymmetric ferromagnets. In this paper, we omit some less important energy contributions such as magnetic anisotropy, stray-field energy, magneto-elastic coupling etc. First, we restrict all calculations to the two-dimensional case on the xy plane. As an instructive example, we consider the DMI energy density inherent to magnets with the C nv symmetry [38]: where ∂ x = ∂/∂x, ∂ y = ∂/∂y. Near the ordering temperatures the magnetization amplitude varies under the influence of the applied magnetic field and temperature. Commonly this process is described by supplementing the magnetic energy (2) with the additional homogeneous free energy term f (M) [37].
where a 1 and a 2 are corresponding coefficients of Landau expansion. The magnetic moment M enters only in even powers due to the time reversal symmetry. The case a 1 = 0 corresponds to the critical temperature T c where spontaneous magnetization appears. Therefore, one can write a 1 in the form: By rescaling the spatial variable, the magnetic field, and the magnetization where energy density can be written in the following reduced form Coefficient a is expressed as Alongside with three internal variables (components of the magnetization vector m) functional (8) includes only two control parameters, the reduced magnetic field amplitude, h, and the "effective" temperature a(T) (9). By direct minimization of Equation (8) one can derive one-dimensional (spirals) and two-dimensional (skyrmions) phases with the propagation directions perpendicular to the polar axis, as dictated by the form of DMI (3). In the present manuscript, however, we analyze only solutions for localized isolated skyrmions (mentioned 1D and 2D bound states will be considered elsewhere).

Equations
Isolated skyrmions can be thought of as isolated static solitonic textures localized in two spatial directions. The magnetization in the center of skyrmion pointing opposite to an applied magnetic field rotates smoothly in all radial directions and reaches the orientation along the field at the outskirt of skyrmion. We will use the spherical coordinates for the magnetization in ISs, and cylindrical coordinates for the spatial variables [27,28], r = (ρ cos ϕ, ρ sin ϕ). Thus, the structure of isolated skyrmions near the ordering temperature is characterized by the dependence of the polar angle θ(ρ) and modulus m(ρ) on the radial coordinate ρ with ψ = ϕ. The total energy E of such a skyrmion after substituting ψ(φ) (see dependences ψ(φ) corresponding to different crystallographic classes in Ref. [1]) is as follows: where energy density is The Euler equations for the functional (8) with boundary conditions describe the structure of isolated skyrmions. The magnetization of the homogeneous phase m 0 is derived from equation: Equation (13) can be solved numerically. However, before to consider typical solutions θ(ρ), m(ρ) of Equation (13) we focus on the asymptotic behavior of skyrmion solutions and some remarkable results that can be drawn by simple means. For 1D kinks such an analysis was done in Refs. [39,40].

Asymptotic Behaviour of Skyrmion Solutions
The asymptotic behaviour of isolated skyrmions bears exponential character [27,28]: By substituting these asymptotes into the linearized Euler Equation (13) one finds three distinct regions in the magnetic phase diagram on the plane (a, h) with different character of skyrmion-skyrmion interactions ( Figure 1): repulsive interactions between isolated skyrmions occur in a broad temperature range (area (I)) and is characterized by real values of parameter α ∈ , the magnetization in such skyrmions has always "right" rotation sense in accord with DMI; at higher temperatures (area (II)) the skyrmion-skyrmion interaction changes to attractive character with α ∈ C; finally, in area (III) near the ordering temperature, a N = 0.25, only strictly confined skyrmions exist with imaginary values of parameter α ∈ . The equation for parameter α obtained from (17) allows to write the equation for the line separating different regions: We note that this equation is written for a general term f . For the case of Landau expansion (4) the separating line looks like h = 2 ± P(a)(a + 1 ± P(a)/2), P(a) = with turning points p (−0.75, √ 2/4), q (0.06, 0.032 √ 5), and u (−0.5, 0) (solid line in Figure 1). Here, we used first and second derivatives, f m (m 0 ) and f mm (m 0 ), calculated from Equation (4)

The Internal Structure of Confined Isolated Skyrmions
To obtain numerically rigorous solutions for isolated skyrmions, one has to solve the system of differential Equation (13) with the boundary conditions (14). We use the "pure" shooting method, where the integration proceeds from the initial point ρ = 0 to a distant point ρ → ∞, and we try to match boundary conditions (14) at the end of the integration. We use the same numerical routines as explicitly described in Ref. [41] (Chapter 17) and therefore omit their description here.
The typical solutions as profiles θ(ρ), m(ρ) for isolated skyrmions in each region are plotted in Figure 1b,c. Since in the region II the exponents α are complex numbers, the profiles display antiphase oscillations (Figure 1c). Rotation of the magnetization in such an isolated skyrmion contains two types of rotation sense: if rotation has "right" sense according to the DMI (blue shaded parts in Figure 1c), the modulus increases, and otherwise, modulus decreases in parts of the skyrmion with "wrong" rotation sense (red shaded parts). Such a unique rotational behaviour of the magnetization is a consequence of the strong coupling between two order parameters of Equation (13)-modulus m and angle θ.
The coupling of angular and longitudinal order parameters may be so strong, that oscillations in the asymptotics of isolated skyrmions do not diminish. The purely imaginary exponent α then reflects the region of strict confinement. For −0.5 < a < 0.25 line (20) delimits a small pocket (III) in the vicinity of the ordering temperature. Within this region skyrmions can exist only as bound states and drastically differ from those in the main part of the phase diagram.

Collapse of Skyrmions at High Fields
The solutions of Equation (13) exist only below a critical line h 0 (T) (Figure 1a). As the applied field approaches this line, the magnetization in the skyrmion center m 1 (Equation (14)) gradually shrinks (Figure 2b), and as m 1 becomes zero, the skyrmion collapses. This is in contrast to low-temperature skyrmions which exist without collapse even at very large magnetic fields [27,28] and are protected by the stiffness of the magnetization modulus which maintains topological stability of skyrmions. At high temperatures the softness of the magnetization amplitude allows to destroy the core of the skyrmion by "forcing" through the magnetization vector m 1 along the applied field. The radial profile for the angle θ(ρ) nevertheless displays an increasing localization under fields, as it was also noted for LT-skyrmions ( Figure 2a). As an example, we illustrated the magnetization process for an isolated skyrmion for a = −0.5 ( Figure 2).
The topological protection of the LT-skyrmions, however, is strictly enforced only in the continuum, whereas the atomic lattice of real materials underlies their collapse. Recently, experimental identification of two distinct skyrmion collapse mechanisms, radial symmetric and chimera-type mechanisms, was reported in biatomic Pd/Fe bilayers with face-centred cubic (fcc) stacking on Ir(111) [42]. In this sense, the mechanism considered in the present manuscript is a high-temperature analogue of a radial symmetric skyrmion collapse in Ref. [42].

Inter-Skyrmion Attraction
For the purpose of investigation of skyrmion-skyrmion interaction, we introduce two skyrmions into a square sample and define the interaction energy per skyrmion ε int versus the inter-skyrmion distance L (Figure 3). Due to the strongly oscillatory character of this dependence two isolated skyrmions will tend to locate at some discrete, equilibrium distances from each other and to be placed in minima of inter-skyrmion energy ε int (Figure 3c). On the other hand, a single isolated skyrmion (minimum corresponding to L = 0) cannot elongate into a pair of skyrmions because of the high potential barrier toward the minimum with finite L (the first deepest minimum of ε int ). Adding skyrmions one by one an optimal number of skyrmions in the cluster is found -the isolated skyrmions tend to form a hexagonal lattice with the densest space filling. The deepest minimum of ε int for two interacting isolated skyrmions is very close to the period of hexagonal skyrmion lattice existing for the same control parameters.

Phenomenon of Skyrmion Confinement
The confinement temperature a L = −0.75 subdivides the temperature interval in low-and high-temperature parts: (i) in the main part (a < a L = −0.75) the rotation of the local magnetization vector determines the chiral modulation, while the magnetization amplitude remains constant; (ii) at high temperatures (a L = −0.75 < a < a N = 0.25) spatial variation of the magnetization modulus becomes a sizeable effect, and strong interplay between longitudinal and angular variables is the main factor in the formation and peculiar behaviour of chiral modulations in this region.
The characteristic temperature a L is of fundamental importance for chiral magnets. It is of the same order of magnitude as the temperature interval ( Figure 1a), where chiral couplings cause inhomogeneous precursor states around the magnetic order temperature (for details see Ref. [32]). Here, a C is the conventional Curie temperature for centrosymmetric systems. When the temperature drops below a C , the energy density of a ferromagnetically spin-aligned state is the lowest one. Dzyaloshinskii-Moriya interactions with negative energy density favour the rotation of the moments. Therefore, the transition to the ferromagnetic state is preceded by a transition to modulated states at the temperature a N . Due to the relativistic origin and corresponding weakness of the DM exchange the shift as well as are small. For MnSi ∆a 1 is estimated to be 0.9 K [32]. The shift ∆a 2 is estimated to be three times as large as ∆a 1 [31,35]. The crossover and confinement effects arise as a generic property of the asymptotic behavior of chiral solitons at large distances from the core. These effects also apply to kinks [39,40] and hopfions [43].

Three-Dimensional Attracting Skyrmions within the Conical Phase of Cubic Helimagnets
For cubic helimagnets mentioned in the introduction Dzyaloshinskii-Moriya interactions are reduced to the following form: In these material systems, isolated skyrmions retain their two-dimensional structure as long as they are surrounded by the homogeneous state for relatively high magnetic fields. In this case, the 2D structure in the plane xy is replicated along the third z direction. At the same time, at some moderate field values, the homogeneous state is replaced by the conical spiral with the propagation vector along the field. The structure of the cone is characterized by the constant polar angle of the magnetization, θ cone = 2h/m, linearly varying azimuthal angle, ψ = z  Color plots of the out-of-plane magnetic moment, m z (x, y) (a), and the modulus, m(x, y) (b), plotted for one cross-section with the fixed coordinate z. In (c), an alternative way of representing the internal structure of a nonaxisymmetric skyrmion is used. After the magnetization components corresponding to the conical phase have been extracted, the 3D-model of the modulus m represents a cylinder-like core centered around the magnetization opposite to the field and exhibiting some reduced value and two coils with the magnetization along the field and enhanced modulus value. (d) two-dimensional cross-cuts across the skyrmion center exhibiting distributions m z (y) and m(y). After averaging over all energy profiles within one conical period, the energy distribution (e) exhibits the following composite parts: the core with the positive energy density, the ring with the negative energy density, and a number of additional rings with alternating energy density.
The internal spin pattern of IS with its axis along the wave vector of the conical phase is depicted in Figure 4. In the field range of cone stability, the distribution of the m z -component in a cross-section xy splits into the central core region that nearly preserves the axial symmetry ( Figure 4a) and the domain-wall region, connecting the core with the embedding conical state. This part of the skyrmion cross-section is asymmetric and accommodates a crescent-like shape, which undergoes additional screw-like modulation along z axis, matching the rotating magnetization of the conical phase. Remarkably, the distribution of the modulus m is more involved (Figure 4b): the modulus exhibits elongation in the crescent part and reduction on the opposite side from the skyrmion center.
A convenient way to depict these non-axisymmetric skyrmions, which has been proven to be particularly illustrative in addressing the character of skyrmion-skyrmion interaction [26], is as follows (Figure 4c): we extract the spins corresponding to the conical phase and then plot the remaining spins as spheres colored according to their magnitude. In this way, all intricate details of the internal structure are explicitly exposed, which makes skyrmion crosscuts along different directions superfluous. In particular, in Figure 4c we depict the 3D-model of the modulus m. From the point of view of the modulus, ISs are composed of a cylinder-like (blue) core centered around the magnetization opposite to the field, a thick (light-red) coil with the magnetization along the field, and somewhat "thinner" (dark-red) coil, which indicates the "ripples" of the magnetization considered before for 2D ISs. The 2D-cross-section along the skyrmion center in Figure 4d shows oscillating m z -component and the modulus m. The oscillatory behavior of the magnetization manifests itself also in the color plot of the energy density averaged over the cone period along z. In Figure 4e, the energy pattern of a LT-skyrmion with the "positive" core in the center and the "negative" ring around it, is complemented by smaller rings with alternating energy density.
The internal structure of LT-ISs embraced by the conical phase has been extensively investigated in our previous papers [25,26]. In the LT-region, only one, the most pronounced coil, persists in 3D-models of the m z -component, which originates from the domain boundary between an isolated skyrmion and the surrounding conical phase. The 3D model of the modulus (Figure 4c) becomes redundant since |m| = const. By theoretical modeling [26] and from the direct experimental observation [25], such skyrmions were shown to attract each other. Then, the cluster formation of such LT-ISs was envisioned as a process of zipping loops: a coil of one skyrmion penetrates the voids between the coils of another one. By this, the magnetization on the way from the center of one skyrmion to the center of another rotates as in ordinary axisymmetric skyrmions. We emphasize that the attracting interaction originates owing to the incompatibility of isolated skyrmions with the surrounding conical phase and from the necessity to adjust their internal structure at each elevation along z axis.
In the "soft" HT version of such isolated skyrmions, in addition to the global energy minimum of skyrmion-skyrmion interaction potential, there should be a number of smaller minima presumably formed owing to the oscillatory asymptotics, as was described before for 2D ISs. Thus, isolated skyrmions near the ordering temperatures become "aware" of each other's presence via oscillatory magnetic ripples. They may start moving towards each other and gather into clusters thus occupying the deepest minimum of the mutual interaction potential.

Conclusions
Using the basic Dzyaloshinskii theory for isotropic chiral magnets, we thoroughly scrutinized the internal structure of isolated 2D skyrmions within the homogeneous state as well as of 3D skyrmions surrounded by the conical spiral in chiral magnets near the ordering temperature. In both cases, the magnetization exhibits coupling of longitudinal and orientational degrees of freedom, which becomes apparent in the fine-tuning of the skyrmion core structures and/or in the oscillatory "ripples" of the magnetization encircling the skyrmion. As a result, the skyrmionic particles develop attractive soliton-soliton interactions and ultimately become confined in extended clusters or textures. The theory and results from numerical simulations demonstrate why a multitude of different small pockets of different phases at the experimental phase diagrams is generically expected in a distinct temperature interval, interleaved between paramagnetic and helix magnetic state. Generically, the attractive couplings between the localized chiral patterns may yield first-order phase transition, between a paramagnetic-like skyrmionic precursor and the helical ground state, as reported for FeGe [34], similarly to the behavior of flux-line lattices in the type-1.5 superconducting systems. Furthermore, the zero-field fully chiral precursor state in MnSi [36] hints at the existence of an amorphous phase of interacting skyrmions, which corresponds to the fully confined region III of the basic phase diagram elucidated here. Hence, the confinement effects of chiral skyrmions strongly change the picture of the formation and evolution of chiral modulated textures and may shed new light on the problem of precursor states observed in chiral magnets.