Theory of current-driven skyrmions in disordered magnets

An emergent topological particle in magnets, skyrmion, has several unique features distinct from the other magnetic textures such as domain wall, helical structure, and vortex. It is characterized by a topological integer called skyrmion number Nsk, which counts how many times the directions of the magnetic moments wrap the unit sphere. This Nsk gives the chiral nature of the skyrmion dynamics, and leads to the extremely small critical current density jc for the current-driven motion in terms of spin transfer torque effect. The finite jc indicates the pinning effect due to the disorder such as impurities and defects, and the behaviors of skyrmions under disorder have not been explored well theoretically although it is always relevant in real systems. Here we reveal by a numerical simulation of Landau-Lifshitz-Gilbert equation that there are four different skyrmion phases with the strong disorder, i.e., (A) pinned state, (B) depinned state, (C) skyrmion multiplication/annihilation, and (D) segregation of skyrmions, as the current density increases, while only two phases (A) and (B) appear in the weak disorder case. The microscopic mechanisms of the new phases (C) and (D) are analyzed theoretically. These results offer a coherent understanding of the skyrmion dynamics under current with disorder.

Skyrmion in magnets is an emergent topological particle made from many spins swirling and pointing all the directions  . These magnetic skyrmions were first found by neutron scattering experiment as the triangular crystal form 7 , while later as individual particles in gas form as well as the crystal by Lorentz transmission electron microscopy 9 . Their motion is driven by the current 22,23 through the spin transfer torque (STT) effect with the extremely small threshold current density ∼ j 10 c 6 A/m 2 compared with that for the domain wall ~10 10 -10 12 A/m 2 24,25 . It is discussed that this high mobility comes from the chiral dynamics of the skyrmions, i.e., the velocity is perpendicular to the force acting on it, analogous to the charged particle under external magnetic field 26 . This makes a skyrmion to avoid an impurity, and reduces the pinning effect. Recently, room temperature skyrmions due to the interface Dzyalosinskii-Moriya interaction [27][28][29] in multi-layered systems [30][31][32][33][34] have been developed. In such systems, it is reported 32,34 that the current driven skyrmions show pinned behavior and sometimes the skyrmion disappears by the current. The defects/impurities seriously influence the current driven skyrmion dynamics. The dynamics of driven skyrmions with disorder has been studied by numerical simulation based on the Thiele equation for each skyrmion, and three phases, i.e., pinned glass, moving liquid, and moving crystal, have been identified in the strong pinned case [35][36][37][38] . However, as the impurity density becomes high and the inter-impurity distance is shorter than the size of the skyrmion and impurity potential becomes strong, the deformation of skyrmions and spin wave emission occur. These phenomena cannot be captured by the Thiele equation for the center of mass motion which assumes the rigid shape of skyrmions. By comparing the results of Thiele equation [35][36][37][38] with ours presented below, one can identify the effects of the internal deformation and spin wave emission, which lead to the different behaviors in the strong disorder case while both are consistent with each other in the weak disorder case.
In this paper, we investigate the current-driven dynamics of skyrmions with disorder by numerically solving the Landau-Lifshitz-Gilbert (LLG) equation for spins. We found four distinct phases as the current density is increased in the strongly disordered case, i.e., (A) pinnned state, (B) depinned state, (C) skyrmion multiplication/ annihilation, and (D) segregation of skyrmions, while only two phases (A) and (B) appear in the weak disorder case. The microscopic processes in each phases are analyzed as shown in Fig. 1, as well as the spectroscopic information.

Model and Simulation.
We focus on the chiral magnet with Dzyaloshinskii-Moriya (DM) antisymmetric spin-orbit interaction [27][28][29] , whose model Hamiltonian on the two-dimensional square lattice is given as i i with x and ŷ being the unit vectors in the x-and y-directions, respectively. The lattice constant is taken as the unit of length. The first term represents the ferromagnetic interaction and the second one is the DM interaction stabilizing the Bloch winding magnetic texture. The easy-axis anisotropy K imp is introduced to the random sites r i ∈ Λ (Λ: set of the random sites), which represents the disorder effect. The Zeeman effect by the external magnetic field h perpendicular to the xy plane is considered by the last term. Let us first discuss the magnetic states without the disorder, i.e., K imp = 0. The competition between the ferromagnetic and the DM interactions results in the single-q helix state with q = D/J, for h = 0. Usually  D J, and therefore the period of the helix ξ = 2π/q = 2π(J/D) is much longer than the lattice constant, which justifies the continuum approximation. The skyrmion crystal state (SkX) appears in the intermediate h between the helix state for low h and the ferromagnetic state for high h 39,40 .
The LLG equation is given by: where α is the Gilbert damping constant (see also Methods). The last two terms in Eq. (2) represent the STT effect due to the spin polarlized electric current density j with the coefficient of the non-adiabatic effect β. We use J for the unit of h, D and K imp , and 1/(γJ) for the unit of time t. (p: polarization of magnet). (See Methods for computational details). Typically ∼ − J 10 3 eV and the unit 1/(γJ) becomes ∼0.7 ps for  γ μ = g / s B (g s : electron spin g-factor, μ B : Bohr magneton). The unit of the electric current density j = |j| is 2eγJ/(pa 2 ) and is typically ∼ . × 1 0 10 13 A/m 2 for the polarization of magnet p = 0.2 and the lattice constant a = 5Å. Figure 1 summarizes the dynamical phase diagrams of the current-driven skyrmions with disorder. In this study, we use D/J = 0.2 and introduce the impurity sites with 5% concentration and investigate the dynamics of skyrmions for various strength of impurity K imp . In this case, a skyrmion has a spatial extent ξ ∼ . The effect of an impurity on the skyrmion discussed here is characterized by an impurity potential with the spatial extent ξ for the center of mass dynamics of the skyrmion. For 5% impurity concentration, because the averaged impurity-impurity distance is about 4 lattice spacings and is much smaller than ξ ∼ 30 lattice spacings, the model discussed here is for the grain defects and not for point defects for the skyrmion dynamics.
We start with the results for h = 0.025 (see Fig. 1(a)). Here, we find the skyrmion pinned state, the skyrmion gas (SkG) state, the skyrmion multiplication behavior and the skyrmion segregation behavior. To investigate these states/behaviors in more detail, we first examine the single skyrmion dynamics in the disordered system. Single skyrmion dynamics: pinned, depinned and multiplication of skyrmion. In Fig. 2(a), the single skyrmion dynamics driven by j = (−j, 0) = (−0.01, 0) in the disordered system for {h = 0.025, K imp = 0.2, α = β = 0.01} is shown. In the left panel of Fig. 2(a), the solid line represents the trajectory of the motion of the skyrmion in the periodic boundary condition. The skyrmion motion is disturbed by impurities and eventually the skyrmion is totally pinned. The final magnetic texture of the pinned skyrmion is shown in the left panel of Fig. 2(a). Under the condition α = β, in the absence of impurities K imp = 0, the Hall angle of the current driven skyrmion motion is zero. In the disordered system, however, as discussed by the previous studies [35][36][37][38]41 , the skyrmion motion shows the Hall behavior during its motion as seen in the trajectory in Fig. 2(a), i.e., on the whole, the skyrmion is moving down in the panel. This is because of the skew dynamics of the skyrmion 1 : The vorticity of the skyrmion drives a dynamics to avoid the impurity similar to the vortices in superconductors 42,43 . Under the successive scattering by the impurities, the Hall behavior is brought about even in the case of α = β in contrast to the clean case 26 . Note that the Hall angle is not ~90 degree with respect the current direction even near the threshold current density different from what is expected from the analytic theory 26,44 .
To evaluate the velocity of the moving skyrmion, we define the quantity, with the emergent e-field 44-46 In the system without impurities K imp = 0, skyrmion has no time-dependent distortion in the current driven motion and it is found that is the skyrmion number and v d = (v d,x , v d,y ) is the velocity of the skyrmion. The total number of skyrmions is given by |N sk | and the single skyrmion magnetic texture {n r } results in N sk = −1 in the present model. The right panels of Fig. 2(a) show the time dependence of the quantities E x /j (blue line) and E y /j (red line). We find a correspondence between the skyrmion trajectory and the plots of E x /j and E y /j in Fig. 2(a): For t ≥ 65000(=t p ), E y /j and E x /j show the damped oscillations around zero reflecting the pinned behavior of the skyrmion. For t ≤ t p , on the other hand, E y /j and E x /j largely fluctuate with large magnitude. From the trajectory of the skyrmion, we calculate the mean velocity of the skyrmion v d for t ≤ t p , and find = . We and it gives (−0.165, −1.079) which is consistent with the relation Eq. (5). Even in the disordered system, E is available to discuss the dynamics of the skyrmion.
We see the pinned behavior for j ≤ 0.02, i.e., with decreasing j, the skyrmion travel-distance decreases and the skyrmion shows a damped oscillation nearby the initial position when j is small enough.
As discussed by previous study 41 , the current driven skyrmion shows the characteristic breathing dynamics. Here, during the traveling, we also see the breathing dynamics and spin wave emission (see also the supplementary movie for Fig. 2). The breathing amplitude becomes larger for larger j. This breathing eventually brings about a substantial distortion of the skyrmion and skyrmion multiplication discussed below.
With increasing the current density j, i.e., for 0.02 < j < 0.2, we find the multiplication of the skyrmion. Figure 2(b) shows a time evolution of the magnetic texture in the skyrmion multiplication for j = 0.1. Under this current density, the moving skyrmion is strongly distorted (see the magnetic texture at t = 75230 in Fig. 2(b)). In the present system, a skyrmion has the skyrmion nunber N sk = −1. In the time-evolution, a part of the skyrmion is squeezed and eventually the skyrmion is separated into two magnetic textures, upper small one and lower elongated one, (see the magnetic texture at t = 75685 in Fig. 2(b)). This occurs within a continuous deformation of the magnetic texture, i.e., it is numerically confirmed that the total skyrmion number N sk does not change along the time-evolution up to t = 75685. At t = 75685, the upper small magnetic texture has N sk = −1 and is a distorted skyrmion (c.f., the skyrmion texture in Fig. 2(a)). The rest of the magnetic texture is a topologically trivial one, i.e., it has N sk = 0: The arrows in the magnetic texture at t = 75685 in Fig. 2(b) indicate the winding texture of the in-plane magnetic moments of the upper end of the lower elongated magnetic texture. This winding texture is not compatible with the DM interaction and is reversed with a spin wave emission along the time-evolution for 75685 Scientific REPORTs | (2018) 8:6328 | DOI:10.1038/s41598-018-24693-5 ≤ t ≤ 75795 (see Fig. 2(b)). At the same time, the skyrmion number changes. The resulting magnetic texture has N sk = −2 in total, i.e., a multiplication of skyrmion occurs.
With further increasing j, i.e., for the current density j ≥ 0.2, the skyrmion multiplication is not seen, and the distortion of the skyrmion becomes small. We find that the skyrmion moves almost along the current direction and the transverse velocity is very small due to the strong STT effect. In addition, we find a substantial spin wave creation/emission by current driven single skyrmion motion in the disordered system with the large current density. Figure 3 shows the time evolution of the emergent e-field for {α = β = 0.01, K imp = 0.2, h = 0.025, j = 0.6} and the black line in Fig. 3(a) represents the trajectory of the skyrmion. (Note that the color code Fig. 3(b) is for e r not for n r ). Because the definition of the e-field Eq. (4) involves the time derivative of the magnetic moments  n r , the spatial distribution and its time evolution give a measure of those of the spin wave(s). As seen in Fig. 3(b), the spatial distribution of the spin wave is not symmetric from the center of the skyrmion. In particular, Fig. 3 clearly shows the spin wave is created and emitted from the moving skyrmion and the large amplitude is concentrated along the trajectory of the skyrmion. It is confirmed that such spin wave creation/emission does not occur in the case that K imp = 0, i.e., in the system without impurities.
Current driven dynamics of many skyrmions. In the absence of impurities K imp = 0, the SkX phase is stabilized for h = 0.025, and the SkX magnetic texture is kept under the current driven motion. At K imp = 0.01, we examine the current driven skyrmion dynamics for j ≥ 0.001 (see Fig. 4). Initial state shown in Fig. 4(a) is the relaxed SkX magnetic texture at K imp = 0.01, and hence the skyrmion arrangement is slightly distorted from the purely triangular one. Figure 4(b) is the magnetic texture under the current j = 0.001 at t = 200000. In this case, the impurity effect is not strong enough to pin all the skyrmions, and the skyrmions show a moving gas like behavior, which is denoted SkG in Fig. 1. For larger j, the impurity effect becomes less important. Figure 4(c) shows a snapshot of the moving skyrmions for j = 1.0 at t = 200000. The (distorted) SkX arrangement seen in the initial state Fig. 4(a) is almost kept during the current driven motion in this case.
For larger K imp , we find the pinned state. The critical current density for K imp = 0.2 is smaller than that for the single skyrmion discussed in Fig. 2. Namely, the skyrmions are moving for j ≥ 0.006 (see Fig. 1(a)) for many skyrmions although the single skyrmion shows the pinned behavior for the current density. This is because the skyrmions kick the other skyrmions before they are eventually pinned during the motion as shown in Fig. 2(a).
To examine the pinned SkX dynamics in more detail, we calculate the Fourier spectra of E x and E y (see Fig. 5): For K imp = 0.2 with h = 0.025, we apply the current j = 0.001 (=j on ). The initial state is the relaxed SkX state in the disordered system. At the first stage, every skyrmions show the damped oscillation nearby the initial positions without large travel-distances, and eventually, the oscillating behavior almost disappears at t = 200000 (=t 1 ). In particular, the magnitudes of E x and E y at t = t 1 are negligibly small. At t = t 1 , we remove the current, i.e., j = j on → 0. Following this, the pinned skyrmions start the relaxation dynamics. In this relaxation dynamics, the skyrmions show again damped oscillation and the oscillating behavior almost disappears at t = 400000 (=t 2 ). For this time interval T = t 2 − t 1 , we calculate the Fourier spectra and the results is shown in Fig. 5. In the same way, the relaxation dynamics for K imp = 0.1 is also examined. Note that this Fourier spectra give the linear response of skyrmions to the external current. Therefore, we expect the resonant behavior around the pinning frequency in analogy to the case of pinned charge density wave [47][48][49][50][51] . It is actually observed in Fig. 5, and the pinning frequency is of the order of ω ∼ 10 3 in Fig. 5(a) with K imp = 0.1, and ω ∼ × − 4 10 3 in Fig. 5

. Clearly it increases as K imp increases.
During the relaxation dynamics in the time interval T, the time-dependent distortion of the each skyrmion is negligibly small. Therefore, the Thiele equation (see Methods) is available to understand the skyrmion dynamics. For > t t 2 , j = 0, so that the Thiele equation gives the velocity of a skyrmion v d,x (t) ≈ −F y /(4π) and v d,y (t) ≈ F x /(4π) by the pinning force F = (F x , F y ) for small α (=0.01 in the present case). Note that E y (t) = −v d,x (t) and E x (t) = v d,y (t) (Eq. (5) with N sk = −1 for a skyrmion), and hence the relation x y x y is seen in Fig. 5 for the relaxation dynamics with tiny oscillation amplitude of the skyrmions at around the pinned skyrmion positions. At large ω limit, the impurity effect becomes less important. In the absence of pinning force, the Thiele equation gives (u(t) is the unit step function) for α = β being the condition discussed here, and hence v d,x (ω) = j on /(−iω). This consequence Eq. (8) is seen as −iωTE y (ω)/(N sk j on ) → 1 and E x (ω) → 0 at the large ω limit in Fig. 5.
We also examine the power spectra |E y (ω)| 2 as shown in Fig. 6 for the depinned skyrmion dynamics. For the current j = 0.008 at K imp = 0.1 and 0.2, skyrmions flow under the impurity effect. This corresponds to the conduction noise of the moving skyrmions. The analogy to the charge density wave reminds us the sharp peak at the "washboard frequency" ω wf = 2πv/λ with v = |v d | being the velocity of skyrmions and λ the period of the SkX state. In the present case, ω wf will be given by λ = | | L N / sk (L: system size) and v = |E(ω = 0)/N sk |, and it is indicated by green arrows in Fig. 6. However, not only at ω = ω wf , but also at whole ω region, we do not find characteristic peak structure in Fig. 6. This is because the particle nature of the skyrmions is more appropriate compared with the density wave picture. There are enough spaces between skyrmions, which behave more like individual particles, and hence there is no well-defined periodicity of the motion. In the log-log plot, |E y (ω)| 2 shows weak ω dependence at small ω and at large ω it drops as ω ∼ − . 6 6 . This behavior is basically consistent with those observed by previous studies [35][36][37][38] .
For K imp = 0.2 and 0.02 < j < 0.4, we find the skyrmion multiplication (i.e, the dark blue region in Fig. 1(a)), similar to the single skyrmion dynamics shown in Fig. 2(b). Figure 7 shows an example of the skyrmion multiplication at the parameter set {h = 0.025, K imp = 0.2, α = β = 0.01, j = 0.04}. Initial state (Fig. 7(a)) is the relaxed SkX state without current, j = 0, in this disordered system. Under the current j = 0.04, skyrmions flow by the STT effect, and due to the impurity effect, strong distortion and multiplication of the skyrmions occur and total number of skyrmions |N sk | increases by almost twice at t = 200000 (see Fig. 7(b) and (d)).
With further increasing j, we find the segregation of the skyrmions. Figure 8 shows the skyrmion dynamics for K imp = 0.2 and j = 0.6. The initial state (see Fig. 8(a)) is the relaxed SkX magnetic texture in this disordered system. Along the current driven skyrmion dynamics, the skyrmions show the stripe formation and finally, the skyrmion stripe along horizontal direction is formed. The important ingredient for the skyrmion segragation is the spin wave emission of the skyrmion discussed in Fig. 3. In the previous study 52 , it is shown that the spin wave attracts the skyrmion. Because the spatial distribution of the spin wave is not symmetric from the center of the skyrmion, the spin wave mediated attractive interaction between the skyrmions brings about the spatially anisotropic skyrmion stripe formation as seen in Fig. 8.
The skyrmion density inside the stripe becomes denser along the time-evolution, and eventually, some of skyrmions are crushed and the skyrmion number decreases. Figures 8(d-f) are the snapshots of the magnetic texture during the typical skyrmion annihilation process. The skyrmion segregation occurs along the current driven motion, so that the skyrmion stripe is relaxed and the skyrmion density becomes sparse after the current is turned off, j = 0. It is well known that the distortion of magnetic structures associated with the motion leads to the inertia, i.e., the mass of the magnetic texture. Therefore, it is an interesting issue weather the mass for the center of mass motion is generated by the disorder for the skyrmions. Figure 9 shows the skyrmion relaxation dynamics after the current is turned off, i.e., j = 0.6(=j on ) for 0 < t ≤ 200000 and j = 0 for 200000 < t. The magnetic texture shown in Fig. 9(a) is identical with Fig. 8(c), but the center of mass in y direction is shifted. (Note that periodic boundary condition is applied here). Figure 9(b) and (c) are the snapshots at t = 204000, and t = 400000, respectively. To examine the skyrmion velocity, the time dependences of E x and E y are also shown in Fig. 9(e) and (f). As discussed in Fig. 3, for the large j, due to the strong STT effect, the skyrmion trajectory shows the small deviation from the skyrmion motion in the absence of the impurities. Therefore, we find tiny E x and E y /j on ≈ N sk for 0 < t ≤ 200000 and this large current density j on = 0.6. At t = 200000, the current is turned off and at the same time, E y is suddenly suppressed. The skyrmion dynamics in this disordered system appears as the damped oscillation of E x and E y with tiny amplitude, and finally all the skyrmions become silent at the trapped positions by the impurities. Therefore, there is no inertia effect for the center of mass motion of skyrmions although the relaxation occurs for the internal structure.
Case of h = 0.04. We also examine the skyrmion dynamics for larger magnetic field h. For the parameter set {J = 1.0, D = 0.2, h = 0.04} corresponding to Fig. 1(b), the perfect ferromagnetic state is the ground state and the SkX state is the metastable state for K imp = 0. As seen in Fig. 1(b), the current driven dynamics is essentially the same as the case in Fig. 1(a), i.e., for h = 0.025, except for absence of the skyrmion multiplication and moving SkX phase with larger j and (very) smaller K imp . The large magnetic field h makes the skyrmion size small. At the same time, the range of the repulsive interaction between skyrmions is also reduced. The SkX with such small repulsive interaction is easily distorted by the current driven motion in the disordered system. For K imp = 0.2 and j ≥ 0.4, we find the clear skyrmion segregation behavior similar to that discussed above. In the skyrmion segregation, the annihilation of the skyrmion also occurs. However, the skyrmion annihilation occurs even for the smaller current density. The example for j = 0.01 is shown in Fig. 10. The clear contrast between the multiplying in Fig. 1(a) and decreasing in Fig. 1(b) of skyrmions in the intermediate current region corresponds to the stability of the skyrmions in the ground state without the current, i.e., skyrmions are stable at h = 0.025 ( Fig. 1(a)) while they are metastable at h = 0.04 ( Fig. 1(b)).

Discussion
We have studied numerically the dynamics of skyrmions with disorder under the current. There are several unique features of this system compared with the other related systems. Vortices in type-II superconductors 42,43 are the most similar system, which also have the chiral dynamics and are driven by the current. However, there are three essential differences. First, the magnetic flux density B cannot terminate or originate because of the absence of the magnetic monopoles and anti-monopoles. Therefore, vortices and B should be considered together as a three-dimensional system even when one treats the two-dimensional superconductor. In contrast, the emergent magnetic field = ⋅ ∂ × ∂ n n n b e ( ) r r r r i i jk j k , 1 2 associated with the skyrmions can have monopoles and antimonopoles, which exist on surfaces in two-dimensional thin film. Furthermore, the magnetic charge is zero for a Bloch type skyrmion, and hence the dipolar interaction can be neglected there. Secondly, the Goldstone boson in superconductor is gapped due to the Higgs phenomenon, i.e., it is shifted up to be the conventional plasmon and cannot mediate the interaction between the vortices. In the case of chiral magnets, the spin wave excitation is in the low energy region even though gapped by the external magnetic field, which mediates the interaction between the skyrmions. Lastly, the dissipation for the vortex motion is usually large, and the inertia is irrelevant, i.e., overdamped motion. On the other hand, the Gilbert damping constant α is usually much smaller than unity, and the system is in the underdamped regime. To compare with the previous work on the dynamics of skyrmions with disorder, which regards each skyrmion with point particle and simulates the Thiele equation, we mention that the internal deformation of skyrmions, multiplication/annihilation of skyrmions, and also the exchange of spin wave between skyrmions, are essential for the strongly disordered case as discussed above. For the future studies, the extension of the present study to the three-dimensional case is an important issue. In that case, the skyrmion becomes a string, the end points of which are monopole and anti-monopole. The dynamics of this string and (anti) monopole in the presence of disorder will offer a rich physics including the manipulation of metastable states 53 .
For the impurity sites r i ∈ Λ (Λ: set of the random sites), we use the random-number generator developed by M. Matsumoto and T. Nishimura (http://www.math.sci.hiroshima-u.ac.jp/~m-mat/MT/emt.html). We have confirmed that the current-driven dynamics is semi-quantitatively similar for different impurity configurations, and hence the impurity average does not modify our conclusions in this paper.
For the single skyrmion, the Thiele equation is expressed to be, sk d d where ẑ is the unit vector perpendicular to the plane forming the two dimensional system studied in this paper, j = (−j, 0) and the tensor κ is defined by κ = ∫d 2 r∂ i n r ∂ j n r (i, j = x, y).