# Trochoidal motion and pair generation in skyrmion and antiskyrmion dynamics under spin-orbit torques

###### Abstract

Skyrmions and antiskyrmions in magnetic ultrathin films are characterised by a topological charge describing how the spins wind around their core. This topology governs their response to forces in the rigid core limit. However, when internal core excitations are relevant, the dynamics become far richer. We show that current-induced spin-orbit torques can lead to phenomena such as trochoidal motion and skyrmion-antiskyrmion pair generation that only occurs for either the skyrmion or antiskyrmion, depending on the symmetry of the underlying Dzyaloshinskii-Moriya interaction. Such dynamics are induced by core deformations, leading to a time-dependent helicity that governs the motion of the skyrmion and antiskyrmion core. We compute the dynamical phase diagram through a combination of atomistic spin simulations, reduced-variable modelling, and machine learning algorithms. It predicts how spin-orbit torques can control the type of motion and the possibility to generate skyrmion lattices by antiskyrmion seeding.

Two-dimensional spin structures such as vortices and skyrmions Bogdanov & Yablonskii (1989); Bogdanov & Hubert (1999) possess a nontrivial topology that affords them a degree of stability Hagemeister et al. (2015); Rohart et al. (2016); Stosic et al. (2017). These structures are characterised by a topological winding number or ‘charge’,

(1) |

where is a unit vector representing the orientation of the magnetic moments in time and space. Skyrmions () and antiskyrmions , for example, possess opposite charges and can appear in pairs through the continuous deformation of the uniform state () Koshibae & Nagaosa (2016); Everschor-Sitte et al. (2017); Stier et al. (2017). The description of the dynamics of skyrmions and antiskyrmions can be approximated by assuming a rigid core, which leads to a reduced set of variables describing their motion. This dynamics is captured by the Thiele equation Thiele (1973),

(2) |

which describes the damped gyrotropic motion of the (anti-)skyrmion core position, , in response to a force . Here is the gyrovector, is a damping constant, and is a structure factor related to the damping (see Methods). While the dynamics in Eq. 2 is non-Newtonian, the gyrotropic response depends on , i.e., its topology, and dictates the direction in which the core moves. This conceptual framework has been useful to understand vortex dynamics Guslienko et al. (2002); Choe et al. (2004), spin-torque vortex oscillators Ivanov & Zaspel (2007); Mistral et al. (2008), and the current-driven motion of skyrmions Sampaio et al. (2013); Nagaosa & Tokura (2013); Lin et al. (2013); Koshibae & Nagaosa (2016); Lin & Hayami (2016); Everschor-Sitte et al. (2017); Leonov & Mostovoy (2017).

In most studies to date, however, the robustness of the symmetry between opposite topological charges, as expressed in Eq. (2), has not been examined in detail. In particular, the roles of core deformation beyond inertial effects Büttner et al. (2015), the internal degrees of freedom, and the underlying symmetry of the magnetic interactions that stabilise the skyrmions remain an open question. This issue is of particular importance since nanometre-scale skyrmions are desirable for possible device applications Fert et al. (2013, 2017) and antiskyrmions have been observed in Heusler compounds Nayak et al. (2017) and predicted to occur at transition metal interfaces Hoffmann et al. (2017). We show here that the symmetries of the magnetic interactions, combined with spin-orbit torques (SOTs), play an important role in determining how the (anti)skyrmion core moves. In particular, the choice of the Dzyaloshinskii-Moriya interaction (DMI) can lead to qualitatively different motion for opposite charges. Namely, deviations from rectilinear motion and skyrmion-antiskyrmion pair generation can occur above certain SOT thresholds for the skyrmion or the antiskyrmion depending on the choice of DMI.

To explore this issue in greater depth, we studied theoretically the spin dynamics of skyrmions and antiskyrmions in an ultrathin 3d transition metal ferromagnet on a 5d normal metal substrate as shown in Figure 1(a).

A prominent example of this material combination is PdFe/Ir(111), where a large DMI is induced in the Fe monolayer through interfacial coupling to the strong spin-orbit interaction in the Ir substrate Fert & Levy (1980); Dupé et al. (2014). This allows individual skyrmions to exist as metastable states, which has been brought to light in recent experiments Romming et al. (2013). Moreover, it has been shown that a variety of antiskyrmion states () are also metastable when frustrated exchange interactions are taken into account in such ultrathin films Dupé et al. (2016); Böttcher et al. (2017); von Malottki et al. (2017); Leonov & Mostovoy (2015) and more generally in bulk chiral magnets Zhang et al. (2017); Hu et al. (2017), which lead to an attractive interaction between skyrmions Rózsa et al. (2016). We employed density-functional theory calculations to obtain estimates of the exchange, anisotropy, and DMI energies for PdFe/Ir(111), which were then used to parametrise an atomistic spin model for studying the dynamics (see Methods). Minimising this energy allows us to determine the equilibrium spin configuration of the static skyrmion and antiskyrmion profiles, as shown in Figure 1(c) and 1(d), respectively. Note that the exchange and DMI possess a six-fold symmetry that is consistent with the Ir(111) surface [Fig. 1(b)].

The spin dynamics are computed by time integrating the Landau-Lifshitz-Gilbert equation with additional SOT terms due to the applied current,

(3) |

where is the Planck constant, is the effective field, is the Gilbert damping constant, is the strength of the field-like torque, and is the strength of the damping-like torque. is the orientation of the effective spin polarisation, which models an applied electric current along the direction in the film plane [Fig. 1(a)]. While in-plane currents should in principle flow through both the ferromagnet and normal metal substrate, we assume that the majority of the current flows only through the substrate since the layer resistivity is significantly larger for the ultrathin ferromagnet (one or two monolayers thick), given the importance of interfacial scattering Sondheimer (2001) and its relative thickness in comparison to the substrate. We can therefore neglect spin transfer torques generated within the 3d ferromagnet and assume only field-like and damping-like contributions from the spin-orbit coupling in the 5d substrate. An example of the ensuing current-driven motion of a skyrmion is shown in Fig. 2a, where the average velocity is plotted as a function of the SOT for the case of . This behaviour is consistent with the Thiele equation, which predicts a linear variation of the skyrmion velocity as a function of the SOT Sampaio et al. (2013).

In contrast to skyrmions, the average velocity of antiskyrmions does not increase monotonically with the SOT [Fig. 2(a)].

A linear regime is found at low currents up to a first threshold, , where a discontinuity in the velocity curve can be seen. Above this threshold, the velocity continues to increase linearly as a function of but with a different slope. A second threshold is found as the strength of the SOT is increased, where the velocity decreases with the applied current. The calculated trajectories for the antiskyrmion core are presented in Fig. 2(b). For linear motion, we observe that the spin configuration of the core is slightly deformed but remains close to its equilibrium static configuration. Above the first threshold , the trajectory is linear at long times but exhibits a large transient phase in which the motion is curved. The rotation ceases when a new steady state regime is reached, which then allows for linear motion to proceed indefinitely (albeit with a different Hall angle with respect to the linear case ). More interestingly, the core undergoes trochoidal motion for which comprises an average displacement along a line that is accompanied by oscillations resulting in loops along the trajectory. The onset of these oscillations results in the sharp decrease in the average velocity shown in Fig. 2(a). The phase diagram of the different behaviour is shown in Fig. 2(c) for different values and ratios of and . We used algorithms based on machine learning to classify the three types of trajectories (linear, deflected, and trochoidal), which exhibit a wide range of velocities and propagation directions (see Methods). We note that the trochoidal motion occurs over a wide range of SOT parameters.

The deflected and trochoidal motion are driven by deformations to the antiskyrmion core. This deformation is characterised by the emergence of a dynamical variable that describes the helicity of the skyrmion and antiskyrmion [Fig. 3(a)].

For skyrmions, describes the continuous transition between Bloch and Néel states of opposite chirality, while for antiskyrmions it describes the rotation of the Bloch or Néel axes. In our system, the deformation is driven by the SOT, which results in a tilt in the magnetization in the film plane, characterised by an amplitude and the azimuthal angle , that depends on the relative strength between the field-like and damping-like terms. This tilt is uniform for the background spins, while it varies within the antiskyrmion core depending on the orientation of . By assuming a suitable ansatz for the deformation profile (see Methods and Supplementary Information), we can derive an equation of motion for using a Lagrangian approach,

(4) |

where is a damping structure factor, is an SOT efficiency factor, is the azimuthal angle of the spin polarisation vector ( in the simulations), and is the internal magnetic energy. Because the effective SOT force acting on in Eq. (2) can be written as

(5) |

the dynamics of determines the time dependence of the force and therefore the overall trajectory of the antiskyrmion as shown in Fig 2(b), and results from the interplay between the SOT term and the restoring force governed by . We verified this interpretation by computing the spatially-resolved forces from the atomistic spin simulations (see Supplementary Information). In Fig. 3(b) and 3(c), we present extracted from the spin dynamics simulations for different SOT strengths. We find that the potential can be described accurately by the function , where and for antiskyrmions, which is consistent with predictions from the model (see Methods). The term represents a lattice effect that accounts for the underlying hexagonal lattice structure Rózsa et al. (2017) and is found to be largely independent of the SOT. The position of the energy minimum is largely independent of the SOT for because the tilt remains almost constant for this torque ratio [Fig. 3(b)]. However, cases where lead to different values of , which results in a shift in the minimum (see Supplementary Information).

From , we can understand the salient features of Eq. (4) as follows. For low amplitudes of the SOT, the restoring force due to the lattice term dominates and the steady state value of remains close to its equilibrium value , resulting in the simple linear motion expected from Eq. (2) alone. As the strength of the SOT is increased, the deformation-induced contribution , which is also governed by the DMI, increases and leads to a change in stability, where a new steady state value is reached. This results in the deflected motion, which is characterised by large transients in leading to the stationary value at long times. As the SOT is further increased, the trochoidal regime is attained when the SOT contribution exceeds the maximum value of the restoring force , which results in a periodic solution in . In this light, the transition toward the trochoidal regime is analogous to Walker breakdown in domain wall motion Slonczewski (1973), where the magnetization angle at the domain wall centre plays the role of here. By using the fits in Fig. 3(b), we computed the dynamics of using Eq. (4) to determine the antiskyrmion trajectories in the three regimes [Fig. 3(d)]. We note that the predicted dynamics of and accurately reproduce the behaviour obtained from the atomistic spin dynamics simulations described by Eq. (3). These results also illustrate why such transitions are not seen for the skyrmion under similar conditions; Fig. 3(c) shows that similar variations in can only be obtained for a skyrmion if the DMI constant is reduced by a factor of , which indicates that the equilibrium skyrmion state is robust and remains largely unperturbed under SOT in this geometry.

Two other regimes beyond the single-particle description are also identified in Fig. 2(c). First, under large field-like and damping-like torques, the propagating antiskyrmion is no longer stable and becomes annihilated. Second, and more interestingly, a transition toward another dynamical regime is found under small and large , where deflected or trochoidal motion leads to a periodic generation of skyrmion-antiskyrmion pairs. An example of this process is given in Figure 4.

This regime is strongly nonlinear and represents a complete breakdown of the single-particle picture described by Eqs. 2 and 4 (see Supplementary Information). The pair is generated as follows: As the antiskyrmion undergoes its trochoidal trajectory, it is accompanied by a large deformation which represents an elongation of the core [i.e., at ps in Fig. 4(a)], similar to the dynamics seen for gyrating magnetic vortices close to the core reversal transition Van Waeyenberge et al. (2006); Yamada et al. (2007); Gaididei et al. (2010). This elongation, which represents a skyrmion-antiskyrmion pair with a net charge of , then separates from the core itself ( to ps). The corresponding topological charge density for these processes is shown in Fig. 4(b). Once nucleated, the pair itself separates since the SOTs lead to different motion for the skyrmion and antiskyrmion constituents. The skyrmion propagates away from the nucleation site by undergoing rectilinear motion, while the nucleated antiskyrmion executes trochoidal motion and becomes itself a new source of pair generation. This remarkable phenomenon leads to the generation of a gas of skyrmions and antiskyrmions [Fig. 4(c)]; the relative population of the two species varies in time as collisions between skyrmions and antiskyrmions lead to annihilation, while pair generation continues for antiskyrmions that survive. This process suggests that it is possible to generate an indefinite number of skyrmions and antiskyrmions from a single antiskyrmion ‘seed’. Combined with the attractive interaction between cores made possible by the frustrated exchange, this dynamics can eventually lead to a skyrmion ‘crystallite’ that condenses from the disordered gas phase (see Supplementary Information). This behaviour is very different to skyrmion generation reported previously, where single pairs are nucleated from static defects through the coupling between local magnetization gradients and spin transfer torques, in systems where skyrmions or antiskyrmions are unstable (depending on the choice of DMI) Everschor-Sitte et al. (2017); Stier et al. (2017).

Recall that we only observed deflected and trochoidal motion for antiskyrmions because the energy barriers in the helicity are orders of magnitude larger for skyrmions than for antiskyrmions for the same DMI constant (Fig. 3). As such, the asymmetry between opposite topological charges is related to the form of the underlying DMI, rather than the sign of the charge itself. To test this hypothesis we conducted simulations in which an anisotropic form of the DMI is used instead, whereby the original six-fold symmetry is retained for the exchange interactions while a two-fold symmetry is used for the DMI, the DMI strength along the two axes being different as shown in Fig. 1(e). This mimics the symmetry of the DMI induced at a (110) interface Hoffmann et al. (2017). Since the amplitudes of the magnetic interactions are unchanged but only the symmetry, the stability of the magnetic textures is only qualitatively affected. Most importantly, antiskyrmions are favoured energetically over skyrmions for this anisotropic form of the DMI.

Figure 5 summarises the current-driven dynamics of skyrmions and antiskyrmions with the anisotropic DMI in Fig. 1(e).

In Fig. 5(a), the current-dependence of the velocity is shown for the skyrmion and antiskyrmion, whose static profiles are shown in Fig. 1(f) and 1(g), respectively. In contrast to the behaviour shown in Fig. 2, the antiskyrmion undergoes only rectilinear motion while the skyrmion exhibits deflected and trochoidal motion as the strength of the SOTs is increased. We note that the associated thresholds, and , are also much higher, where velocities beyond 1 km/s can be reached in the linear regime for both skyrmions and antiskyrmions. In Fig. 5(b), examples of trajectories for the linear, deflected, and trochoidal motion for skyrmions are shown. We note that the overall skyrmion Hall angles are different to the antiskyrmion case shown in Fig. 2(b), which originates from different stationary values of for the skyrmion. This is a consequence of for the skyrmion with the anisotropic DMI, which possesses a different -dependence than the case shown in Fig. 3(b). This difference is also reflected in the phase diagram in Fig. 5(c); while the same phases are identified, the overall shape of the phase boundaries differs and certain transitions are absent, such as the transition between deflected motion and pair generation. Nevertheless, the order in which the phases appear with increasing SOT is similar.

The importance of the DMI symmetry can be highlighted further by examining the current-driven dynamics in the absence of DMI altogether. Skyrmions and antiskyrmions remain metastable states because of the frustrated exchange interactions, resulting in the equilibrium profiles shown in Fig. 6(a) and 6(b).

Since the absence of a chiral interaction results in Bloch and Néel states being degenerate. is a constant in Eq. (4), so the internal mode becomes a Goldstone mode of the system that can be excited with vanishingly small torques. The Bloch-like skyrmion profile in Fig. 6(a) is therefore only one possible realization of the metastable state. For a finite deformation , according to Eq. (4) and generates a harmonic SOT force in Eq. (2), which results in circular motion. This is confirmed in the spin dynamics simulations as shown in Fig. 6(c), where circular motion is indeed found with an opposite sense of rotation for opposite topological charges, as expected from the sign of the gyrovector in Eq. 2. Since the deformation is linear in SOT for the range of values considered, we expect a quadratic variation in the gyration frequency as a function of SOT from Eq. (4). This is confirmed in Fig. 6(d), where the simulated frequencies are well described by a quadratic function. A similar analysis predicts that the radius of gyration should be inversely proportional to the SOT, which is again confirmed by simulations as shown in Fig. 6(e). Finally, the SOT dependence of the tangential velocity is presented where Fig. 6(f), where a linear variation is found in agreement with theory. Besides the opposite sense of gyration, these results show that the skyrmion and antiskyrmion trajectories agree quantitatively within the numerical accuracy of the simulations.

These results highlight the rich dynamical behaviour that is possible under SOTs in ultrathin ferromagnetic films, especially for metastable chiral states that are not necessarily the most energetically favourable. The work also links skyrmion dynamics to other known phenomena in micromagnetism, namely Walker breakdown in domain wall motion (trochoidal motion) and vortex core reversal (pair generation). Given the primacy of the DMI symmetry in governing the particle dynamics, our work may spur new avenues of research in materials science where specific surface or interface orientations could be chosen to tailor particular dynamical properties, such as deflected or trochoidal motion, which is absent in most approaches where the focus is on quantifying and controlling rectilinear motion for skyrmion memory and logic applications. The prospect of generating different dynamics with a variety of metastable states within the same material system could also offer new possibilities for studying particle interactions and developing new application paradigms, notably skyrmion generation with a single antiskyrmion ‘seed’. Recent theoretical work shows that such seeds are likely to appear at finite temperatures Böttcher et al. (2017) and therefore offer a reliable and efficient means of producing skyrmions and antiskyrmions readily.

### Methods

#### Hamiltonian

The magnetic Hamiltonian studied is given by

(6) |

where the first term represents the Heisenberg exchange interaction, the second term the Dzyaloshinskii-Moriya interaction (DMI), the third the uniaxial anistropy along the -axis, and the last term the Zeeman energy associated with an external field . The indices in the summation for the exchange and DMI terms indicate that single site terms are neglected. The moments are assumed to reside on a hexagonal lattice and everywhere. The parameters are extracted from density functional theory (DFT) calculations of the bilayer PdFe system on Ir(111) Dupé et al. (2014, 2016), in which we consider a fcc stacking for the Pd layer. For the Heisenberg exchange, represents the exchange constant between the magnetic moments and , where up to 10 nearest-neighbours are taken into account: meV, meV, meV, meV, meV, meV, meV, meV, meV, and meV. We treat the DMI in the nearest-neighbour approximation as shown in Fig. 1(b) and 1(e), where a magnitude of meV for is obtained from DFT calculations. The anisotropy constant is meV and we used an applied magnetic field of 20 T along the direction. The magnetic moment of the Fe atoms is given by , with being the Bohr magneton. For the given parameters the system is in a ferromagnetic ground state close to the transition point to the skyrmion lattice phase, where isolated skyrmions and antiskyrmions can be stabilised. The applied magnetic field is only slightly larger than the critical field , with .

#### Atomistic spin dynamics simulations

The simulation geometry comprises a hexagonal lattice of 100100 spins with periodic boundary conditions. The ferromagnet is assumed to be one monolayer thick. The dynamics of the spin system described by Eq. 6 is solved by numerical time integration of the Landau-Lifshitz equation with Gilbert damping and spin-orbit torques given in Eq. 3. We used a Gilbert damping constant of for all the simulations presented here. The numerical time integration is performed using the Heun method. At the start of each simulation, an equilibrium skyrmion or antiskyrmion profile is first computed by relaxing the system in the absence of the SOT terms. This procedure produces the profiles shown in Figs. 1(c), 1(d), 1(f), 1(g), 6(a), and 6(b). The simulations are then executed over several ns with a fixed time step in the range of fs.

#### Extension to Thiele model

The extension to the Thiele model, expressed by Eq. 4, is based on the idea that spin-orbit torques (SOT) lead to a significant deformation of the skyrmion/antiskyrmion core. The model is based on two assumptions. First, we assume that all spins in the system are canted toward the film plane under the combined action of the field-like and damping-like SOT. The deformation is assumed to take the form , where the relaxed ground state is and

with representing the amplitude of the deformation and describing the azimuthal component of the background spins that tilt away from the -axis as a result of the SOT. Second, in addition to the core position , we elevate the helicity parameter to a dynamical variable, which is defined through the azimuthal angle , where is the topological charge. Based on this deformation ansatz, we derive the equation of motion for using a Lagrangian approach Kim (2012), which involves a continuum approximation for the magnetization, , with . By neglecting coupling terms proportional to , we derive the Euler-Lagrange equations leading to Eqs. (2) and (4), where the gyrovector term is given by,

(7) |

the damping factors are

(8) | ||||

(9) |

and the SOT efficiency factors are

(10) | ||||

(11) |

Here, the equilibrium (anti-)skyrmion core profile is assumed to possess a cylindrical symmetry, with being the radial variable in cylindrical coordinates.

Expressions for the helicity-dependent energy, , can be found in a similar way by using the continuum approximation of Eq. (6). The dominant contribution comes from the DMI. For the symmetry considered in Fig. 1(b), we use the form Bogdanov & Rößler (2001); Thiaville et al. (2012)

(12) |

where is the DMI constant. For skyrmions, we find

(13) |

where

(14) | ||||

(15) |

Here, is the dominant term, while the deformation-induced contribution provides a correction that increases quadratically with the deformation. Only the deformation-induced contribution appears for the antiskyrmion,

(16) |

where

(17) |

As noted in the main text, an additional energy term is required to describe the atomistic simulations, which arises from discretisation effects due to the underlying hexagonal lattice. This lattice term is not present in the continuum description.

#### Classification of skyrmion and antiskyrmion trajectories

The deflected and trochoidal motion for antiskyrmions [with the DMI in Fig. 1(b)] and skyrmions [with the DMI in Fig. 1(e)] can involve a wide range of speeds, propagation directions, and gyration frequencies. Classifying these behaviours efficiently from simulation data in order to construct phase diagrams shown in Figs. 2(c) and 5(c) is therefore a challenging task. We employed algorithms based on machine-learning to classify these trajectories, which were then used with adaptive meshing to identify the different phase boundaries. First, we excluded the annihilation and pair-generation states from the simulation data, which could be identified directly from the magnetization state. Second, the velocity orientations within each simulation run for the remaining data were mapped onto the unit circle, which then served as inputs for classification. The linear motion results in a small cluster of points on the circle, the deflected motion gives a partially filled circle, while the trochoidal motion results in a fully filled circle. A subset of these images (5-10 per state) were then used as learning rules to train the Classify function in the technical computing software Mathematica (version 11.2), which was then used to classify the remaining states. The target resolutions of the phase boundaries in Figs. 2(c) and 5(c) are 0.01 meV and 0.02 meV, respectively. A brute force search would therefore require 22 500 () simulation runs for each DMI symmetry, while our iterative method combined with machine learning required only 1831 [Fig. 2(c)] and 2736 [Fig. 5(c)] runs, respectively. Given that each simulation run takes 5 to 10 hours of computation time on a single central processing unit (CPU) core, our method provides a more efficient way to explore the parameter space of the dynamical system.

### References

## References

- Bogdanov & Yablonskii (1989) Bogdanov, A. & Yablonskii, D. Contribution to the theory of inhomogeneous states of magnets in the region of magnetic-field-induced phase transitions. Mixed state of antiferromagnets. Zh. Eksp. Teor. Fiz 69, 142–146 (1989).
- Bogdanov & Hubert (1999) Bogdanov, A. & Hubert, A. The stability of vortex-like structures in uniaxial ferromagnets. Journal of Magnetism and Magnetic Materials 195, 182–192 (1999).
- Hagemeister et al. (2015) Hagemeister, J., Romming, N., Von Bergmann, K., Vedmedenko, E. Y. & Wiesendanger, R. Stability of single skyrmionic bits. Nature Communications 6, 8455 (2015).
- Rohart et al. (2016) Rohart, S., Miltat, J. & Thiaville, A. Path to collapse for an isolated Néel skyrmion. Physical Review B 93, 665–6 (2016).
- Stosic et al. (2017) Stosic, D., Mulkers, J., Van Waeyenberge, B., Ludermir, T. B. & Milošević, M. V. Paths to collapse for isolated skyrmions in few-monolayer ferromagnetic films. Physical Review B 95, 214418–12 (2017).
- Koshibae & Nagaosa (2016) Koshibae, W. & Nagaosa, N. Theory of antiskyrmions in magnets. Nature Communications 7, 10542 (2016).
- Everschor-Sitte et al. (2017) Everschor-Sitte, K., Sitte, M., Valet, T., Abanov, A. & Sinova, J. Skyrmion production on demand by homogeneous DC currents. New Journal of Physics 19, 092001 (2017).
- Stier et al. (2017) Stier, M., Häusler, W., Posske, T., Gurski, G. & Thorwart, M. Skyrmion–Anti-Skyrmion Pair Creation by in-Plane Currents. Physical Review Letters 118, 267203 (2017).
- Thiele (1973) Thiele, A. A. Steady-state motion of magnetic domains. Physical Review Letters 30, 230 (1973).
- Guslienko et al. (2002) Guslienko, K. Y. et al. Eigenfrequencies of vortex state excitations in magnetic submicron-size disks. Journal of Applied Physics 91, 8037 (2002).
- Choe et al. (2004) Choe, S. B. et al. Vortex Core-Driven Magnetization Dynamics. Science 304, 420–422 (2004).
- Ivanov & Zaspel (2007) Ivanov, B. & Zaspel, C. Excitation of Spin Dynamics by Spin-Polarized Current in Vortex State Magnetic Disks 99, 247208 (2007).
- Mistral et al. (2008) Mistral, Q. et al. Current-Driven Vortex Oscillations in Metallic Nanocontacts. Physical Review Letters 100, 257201 (2008).
- Sampaio et al. (2013) Sampaio, J., Cros, V., Rohart, S., Thiaville, A. & Fert, A. Nucleation, stability and current-induced motion of isolated magnetic skyrmions in nanostructures. Nature Nanotechnology 8, 839–844 (2013).
- Nagaosa & Tokura (2013) Nagaosa, N. & Tokura, Y. Topological properties and dynamics of magnetic skyrmions. Nature Nanotechnology 8, 899–911 (2013).
- Lin et al. (2013) Lin, S.-Z., Reichhardt, C., Batista, C. D. & Saxena, A. Driven Skyrmions and Dynamical Transitions in Chiral Magnets. Physical Review Letters 110, 207202 (2013).
- Lin & Hayami (2016) Lin, S.-Z. & Hayami, S. Ginzburg-Landau theory for skyrmions in inversion-symmetric magnets with competing interactions. Physical Review B 93, 064430 (2016).
- Leonov & Mostovoy (2017) Leonov, A. O. & Mostovoy, M. Edge states and skyrmion dynamics in nanostripes of frustrated magnets. Nature Communications 8, 14394 (2017).
- Büttner et al. (2015) Büttner, F. et al. Dynamics and inertia of skyrmionic spin structures. Nature Physics 11, 225–228 (2015).
- Fert et al. (2013) Fert, A., Cros, V. & Sampaio, J. Skyrmions on the track. Nature Nanotechnology 8, 152–156 (2013).
- Fert et al. (2017) Fert, A., Reyren, N. & Cros, V. Magnetic skyrmions: advances in physics and potential applications. Nature Reviews Materials 2, 17031 (2017).
- Nayak et al. (2017) Nayak, A. K. et al. Magnetic antiskyrmions above room temperature in tetragonal Heusler materials. Nature 548, 561–566 (2017).
- Hoffmann et al. (2017) Hoffmann, M. et al. Antiskyrmions stabilized at interfaces by anisotropic Dzyaloshinskii-Moriya interaction. Nature Communications 8, 308 (2017).
- Fert & Levy (1980) Fert, A. & Levy, P. M. Role of Anisotropic Exchange Interactions in Determining the Properties of Spin-Glasses. Physical Review Letters 44, 1538–1541 (1980).
- Dupé et al. (2014) Dupé, B., Hoffmann, M., Paillard, C. & Heinze, S. Tailoring magnetic skyrmions in ultra-thin transition metal films. Nature Communications 5, 4030 (2014).
- Romming et al. (2013) Romming, N. et al. Writing and Deleting Single Magnetic Skyrmions. Science 341, 636–639 (2013).
- Dupé et al. (2016) Dupé, B., Kruse, C. N., Dornheim, T. & Heinze, S. How to reveal metastable skyrmionic spin structures by spin-polarized scanning tunneling microscopy. New Journal of Physics 18, 055015 (2016).
- Böttcher et al. (2017) Böttcher, M., Heinze, S., Sinova, J. & Dupé, B. Thermal formation of skyrmion and antiskyrmion density arXiv:1707.01708 [cond–mat.mtrl–sci] (2017).
- von Malottki et al. (2017) von Malottki, S., Dupé, B., Bessarab, P. F., Delin, A. & Heinze, S. Enhanced skyrmion stability due to exchange frustration. Scientific Reports 7, 12299 (2017).
- Leonov & Mostovoy (2015) Leonov, A. O. & Mostovoy, M. Multiply periodic states and isolated skyrmions in an anisotropic frustrated magnet. Nature Communications 6, 8275 (2015).
- Zhang et al. (2017) Zhang, X. et al. Skyrmion dynamics in a frustrated ferromagnetic film and current-induced helicity locking-unlocking transition. Nature Communications 8, 1717 (2017).
- Hu et al. (2017) Hu, Y., Chi, X., Li, X., Liu, Y. & Du, A. Creation and Annihilation of Skyrmions in the Frustrated Magnets with Competing Exchange Interactions. Scientific Reports 7, 16079 (2017).
- Rózsa et al. (2016) Rózsa, L. et al. Skyrmions with Attractive Interactions in an Ultrathin Magnetic Film. Physical Review Letters 117, 157205 (2016).
- Sondheimer (2001) Sondheimer, E. H. The mean free path of electrons in metals. Advances in Physics 50, 499–537 (2001).
- Rózsa et al. (2017) Rózsa, L. et al. Formation and stability of metastable skyrmionic spin structures with various topologies in an ultrathin film. Physical Review B 95, 094423 (2017).
- Slonczewski (1973) Slonczewski, J. C. Theory of domain-wall motion in magnetic films and platelets. Journal of Applied Physics 44, 1759–1770 (1973).
- Van Waeyenberge et al. (2006) Van Waeyenberge, B. et al. Magnetic vortex core reversal by excitation with short bursts of an alternating field 444, 461–464 (2006).
- Yamada et al. (2007) Yamada, K. et al. Electrical switching of the vortex core in a magnetic disk. Nature Materials 6, 270–273 (2007).
- Gaididei et al. (2010) Gaididei, Y., Kravchuk, V. P. & Sheka, D. D. Magnetic vortex dynamics induced by an electrical current. International Journal of Quantum Chemistry 110, 83–97 (2010).
- Kim (2012) Kim, J.-V. Spin-torque oscillators. In Solid State Physics (eds. Camley, R. E. & Stamps, R. L.), 217–294 (Academic Press, 2012).
- Bogdanov & Rößler (2001) Bogdanov, A. N. & Rößler, U. K. Chiral Symmetry Breaking in Magnetic Thin Films and Multilayers. Physical Review Letters 87, 037203 (2001).
- Thiaville et al. (2012) Thiaville, A., Rohart, S., Jué, É., Cros, V. & Fert, A. Dynamics of Dzyaloshinskii domain walls in ultrathin magnetic films. Europhysics Letters 100, 57002 (2012).
- Vansteenkiste et al. (2014) Vansteenkiste, A. et al. The design and verification of MuMax3. AIP Advances 4, 107133 (2014).
- Braun (1994) Braun, H.-B. Fluctuations and instabilities of ferromagnetic domain-wall pairs in an external magnetic field. Physical Review B 50, 16485–16500 (1994).
- Romming et al. (2015) Romming, N., Kubetzka, A., Hanneken, C., von Bergmann, K. & Wiesendanger, R. Field-Dependent Size and Shape of Single Magnetic Skyrmions. Physical Review Letters 114, 177203 (2015).

### Acknowledgements

This work was partially supported by the Horizon2020 Framework Programme of the European Commission under Grant No. 665095 (MAGicSky). JK acknowledges support from the Deutscher Akademischer Austauschdienst under Award No. 57314019. UR acknowledges support from the Deutsche Forschungsgemeinschaft (grant RI2891/1-1). UR, BD and JS acknowledge the Alexander von Humboldt Foundation, the Deutsche Forschungsemeinschaft (grant DU1489/2-1), the Graduate School Materials Mainz, the ERC Synergy Grant SC2 (No. 610115), the Transregional Collaborative Research Center (SFB/TRR) 173 SPIN+X, and the Grant Agency of the Czech Republic grant no. 14-37427G.

### Author Contributions

BD and SH initiated the project. UR and SvM developed the atomistic spin dynamics code and UR performed the atomistic spin dynamics simulations. UR and JVK interpreted the simulation results and developed the analytical model. SH, BD, JVK and UR wrote the manuscript. All the authors discussed the data.

### Competing financial interests

The authors declare no competing financial interests.

### Materials & Correspondence

Correspondence and requests for materials should be addressed to Ulrike Ritzmann

(ulrike.ritzmann@physics.uu.se).

## Supplementary Information

## Appendix S1 Force density

To obtain further insights in the role of the deformation on the dynamics of skyrmions and antiskyrmions, we calculate the force density and the total force acting on the magnetic texture during the motion following the methodology of Thiele ^{9}. The force density is given by , where the field term includes contributions from the effective field , the Gilbert damping and SOT:

(S1) |

In the case of a relaxed state in the absence of SOT terms, the force density is negligible. Internal force contributions from DMI, exchange or the magnetic field are large inside the texture, but all contributions are compensating each other. In the presence of SOT, the force density inside the magnetic texture is no longer vanishing and the motion due to the total force can be described by Eq. (2) in the linear regime. An example of the force density is shown in Fig. S1a) for the linear motion of an antiskyrmion due to SOT with the strength meV at a time of ps. In that case, the total force points towards - direction while the antiskyrmion moves in -direction.

The absolute value of the total force dependent on the time is shown in Fig. S1b). In the case of linear motion, the force density remains constant, which results in a constant velocity. In agreement with our analytical model, the force density becomes time-dependent above the threshold for deflected motion. In the regime of the deflected motion, in which the antiskyrmion deforms until it reaches a steady state, the total force decreases at the beginning and remains constant afterwards. For even higher torques in the regime of trochoidal motion, the force is changing periodically in time as illustrated in Fig. S1b) for SOT of the strength of meV. In this case, the total force and the speed during the trochoidal motion are linked to each other and change periodically with the same frequency.

## Appendix S2 Total forces and torques for pair creation

The rigid body approximation, which is used to describe the motion in the regime of linear, deflected and trochoidal motion, is no longer valid for pair creation. Nevertheless, the calculation of the acting forces can be used to explain the details of the process of pair creation. Fig. S2 shows the absolute value of the total force and the sum of the appearing torques for pair creation for SOT with meV and meV, complementing to the magnetization profiles shown in Fig. 4. A strong deformation of the antiskyrmion is observed in the simulations around 3ps leading to the creation of the pair around 4ps. This is reflected in a clear increase of the total force above 3ps, which then decreases after the pair is formed around 4ps. Moreover, we study the sum of the torques acting on the magnetic moments. In the presence of SOT after an initial relaxation of the ferromagnetic background, mainly torque terms are contributing that are acting on the magnetic moment around the texture. In the lower part of Fig. S2, the absolute value of the total torque is shown for pair creation. Here, the torque is reduced during the pair creation process.

## Appendix S3 Exchange and DMI energy density distribution during the pair creation

To understand the creation of a skyrmion-antiskyrmion from a single antiskyrmion, we show on figure S3 the exchange and the DMI energy density distribution for three different time: 3.0 ps, 3.8 ps and 4.2 ps corresponding to Fig. 4. The panel (a) shows the exchange energy at 3.0 ps. The skyrmion profile is already non-spherical and the exchange energy shows two distinct maximum. This Exchange energy profile does not correspond to the one of an antiskyrmion which only shows a single maximum in the exchange energy density. The presence of 2 maxima can be viewed as the beginning of the creation process. The panel (b) corresponds to the DMI energy density at 3.0 ps. The DMI energy density corresponds to the profile of an antiskyrmions. The antiskyrmion is composed of left- and right-rotating magnetic regions depending on the axis corresponding to a positive and a negative DMI energy density, respectively. However, the region of left-rotating magnetic texture is showing a lobe and extends away from the antiskyrmion core. This shows that the SOT is enhancing the stability of the left rotating region when is aligns to the polarization of the current.

At 3.8 ps (panel c), the exchange energy shows two distinct maxima which are enclosed within 2 distinct isolines corresponding to . At this time, the DMI energy density distribution exhibits a region characteristic for the antiskyrmion (centered at coordinate (0;-6)) and a region corresponding to a skyrmion-antiskyrmion pair centered at (-2.5;4). The center of the skyrmion corresponds to the negative DMI energy distribution present at (1,5) visible on panel d. The antiskyrmion starts to form below at (-6;2.5) and a right-rotating magnetic texture start to form between this antiskyrmions and the skyrmions.

At 4.2 ps (panel e), the exchange energy density shows 2 regions of positive contribution corresponding to the presence of 2 antiskyrmions and one region of negative contribution corresponding to the skyrmion. The isolines are now forming 3 distinct regions in space. The DMI energy density distribution (panel f) confirms the presence of 2 antiskyrmions (exhibiting a butterfly contrast) and one homogeneous negative energy density corresponding to the skyrmion. It is interesting to notice that the left-rotating magnetic regions of the antiskyrmions are aligned but the skyrmion faces one of the right-rotating magnetic region of the antiskyrmion.

## Appendix S4 Total force density distribution during pair creation

The previous section describes the evolution of the energy density distribution during the creation of a skyrmion-antiskyrmion pair. In this section, the pair creation is described in term of forces as defined in section S1.

At 3.0 ps (panel a), the forces are maximum within the antiskyrmion region corresponding to the isolines . Although the total energy density distribution shows 2 distinct maxima (panel a), the forces show a strong maximum for one half of the deformed antiskyrmion (panel b).

At 3.8 ps, the total energy distribution (panel c) shows two distinct maxima in 2 separate regions of space. The direction of the forces are pointing towards the regions showing a maximum of total energy distribution. This shows that two distinct particle-like magnetic textures are present. The amplitude of the forces (panel d) also shows 2 distinct regions. One is corresponding to the antiskyrmion and shows a strong maximum of the amplitude of the forces density distribution. In the region where the skyrmion-antiskyrmion pair is created, the amplitude of the force starts to rise.

At 4.2 ps, the total energy density shows 3 separated maxima which are enclosed within separate region of space corresponding to the isolines . These 3 distinct maxima show that 3 particle-like non-collinear magnetic textures are now present. The direction of the forces is interesting because it can characterize the displacement. In this region of the phase diagram, the antiskyrmion have a trochroidal trajectory which means that they should have a non-zero angular momentum. Indeed, the forces are pointing toward the cores of the antiskyrmions which means that . On the contrary, the skyrmions which have a linear motion, show a force density distribution which is collinear and pointing in the same direction. A closer look on the amplitude of the force density distribution (panel e) shows that there are 3 different regions indicating that the different skyrmions start to get into motion.

## Appendix S5 Analysis of for antiskyrmions

In this section, we present some details of the numerical analysis of the energy extracted from spin dynamics simulations. This discussion complements Fig. 3 of the main text. In Figure S5, we present supplementary data for for different values of the SOT terms.

In Fig. S5(a) where the case is considered, the data have been fitted with the function

(S2) |

where the first term is the deformation-induced term that depends on the DMI and the second term is a lattice term, as discussed in the main text. As the SOT is increased, we observe that the energy landscape becomes predominantly governed by the cosine term, as expected from the theory. We note that the energies are well described by this function for the values of the SOT considered. Moreover, the position of the global extrema do not depend on the SOT, which reflects the fact that the overall tilt varies little along the contour . A shift in the positions of the extrema, along with the overall shape of the energy landscape, can be seen for cases where as shown in Fig. S5(b). This is consistent with the extended Thiele model.

The SOT dependence of the coefficients , , and for is shown in Figs. S5(c)-(e). is found to be well described by a quadratic function, which is consistent with the model prediction

(S3) |

if we assume that the deformation is linear in SOT. Here is the Dzyaloshinskii-Moriya coefficient and is a numerical prefactor that depends on the details of the shape of the antiskyrmion core (see Methods section of the main text). We note that the linear dependence of on the SOT can be expected from simple field arguments and is consistent with the behaviour discussed in Fig. 6 of the main text. The lattice term is found to be largely independent of the SOT, with deviations only seen at small (0.04 meV) and larger values (0.2, 0.3 meV) of the SOT. The discrepancy for the former may arise from uncertainties in the fit, since the linear motion does not lead to a large range of values being sampled. For the larger values, we note that SOT-induced core deformations become more important, where the core itself loses much of its cylindrical symmetry. Nevertheless, the general trends are consistent with the model. Similarly, the fitted coefficient is also found to be largely independent of the SOT, which is consistent with the fact that varies little with the SOT for .

## Appendix S6 Pair generation and skyrmion crystallisation

In this section, we provide further details of the skyrmion-antiskyrmion pair generation, as shown in Fig. 4 in the main text. In Fig. S6, we illustrate a particular example in which the generation leads to a skyrmion ‘crystallite’ that condenses from a dilute gas of skyrmions.

Here, the SOT is dominated by the damping-like term ( meV, meV). The initial state of the simulation comprises an antiskyrmion at the centre of the simulation grid. Under the applied SOT, the antiskyrmion undergoes a trochoidal motion, leading first to a large core deformation and then a nucleation of the skyrmion-antiskyrmion pair, as shown at ns in Figs. S6(a) and S6(b). This pair then separates under the action of the SOT and is followed by the periodic generation of other pairs. Interestingly, the nucleated antiskyrmions do not appear to escape the region of generation and are annihilated more readily than their skyrmion counterparts, leading to an excess in the skyrmion population as shown at and ns. This is also quantified in Fig. S6(c), where the total topological charge is shown as a function of time. An event occurs after ns that annihilates the seed antiskyrmion, resulting in a dilute gas of pure skyrmions. Because of the attractive interaction made possible by the frustrated exchange, this gas progressively condenses into a ‘crystallite’ as time evolves, where successive collisions lead to the ordered configuration shown at ns in Fig. S6.

## Appendix S7 SOT-induced tilts and model for core deformation

In this section, we describe the magnetization tilts induced by the spin-orbit torques and how these are used to motivate the ansatz for the core deformation profile. The influence of the field-like and damping-like torques on the uniform ferromagnetic state is presented in Fig. S7.

The equilibrium configuration was computed for different values of the SOT (with ) in the absence of a skyrmion or antiskyrmion with the magnetic Hamiltonian given in the Methods section of the main text. Overall, the combined action of the field-like and damping-like torques can be assimilated to the presence of a fictitious SOT magnetic field applied in the film plane, which tilts the magnetization away from the -axis, characterised by the polar angle , along an azimuthal orientation given by the angle taken from the -axis in the usual way. We note that the polar tilt is largely independent of the ratio, which is expected for an applied magnetic field in the film plane [Fig. S7(a)]. On the other hand, the azimuthal tilt depends primarily on this ratio, rather than the magnitudes of the SOT terms, which can also be expected from this field argument [Fig. S7(b)]. The polar tilt is found to vary linearly as a function of SOT for the range of values considered [Fig. S7(c)], while the azimuthal tilt is largely independent of the SOT strength above a threshold for the case [Fig. S7(d)].

Our ansatz for the core deformation is based on this idea that the SOT acts like an effective magnetic field on the equilibrium state. We first discuss the macrospin approximation. Let be a unit vector representing the magnetization orientation, which we can write in terms of the spherical polar angles in the usual way, . Let us consider the effect of an applied magnetic field in the film plane, . Under the action of the Landau-Lifshitz damping torque, the time rate of change in can be written as

(S4) |

By assuming and , this leads to the two linearly independent equations,

(S5) | ||||

(S6) |

Therefore, we can assume that applied in-plane fields will lead to tilts of the form

(S7) | ||||

(S8) |

where is a parameter that describes the tilt amplitude. We then extend this approach to describe the tilting of any arbitrary magnetization field by assuming that and . Let us consider deformations about an equilibrium configuration, , due to the same uniform in-plane applied field. To linear order in the deformations, the magnetization field can be written as

By using the expression for the tilts in Eqs. (S7) and (S8), we obtain

(S9) | ||||

(S10) | ||||

(S11) |

which is the model described in the main text, where we have substituted the SOT polarization angle in the place of here. Note that this deformation ansatz does not conserve the length of the magnetization vector,

(S12) |

This means that needs to be rescaled when making comparisons with simulation data.

To explore the validity of this ansatz, we performed micromagnetics simulations of a static skyrmion in the presence of an applied field along the direction. Micromagnetics simulations were used (rather than the atomistic code) as they provide a better correspondence with the continuum approximation in which this deformation model is developed. We used the MuMax3 code ^{43} and considered a nm film that was discretised with finite difference cells. Periodic boundary conditions were used. We assumed an exchange constant of pJ/m, a uniaxial anisotropy of MJ/m, a saturation magnetization of MA/m, and an interfacial DMI of mJ/m. The relaxed skyrmion profile can be well described with the double soliton model ^{44, 45},

(S13) |

where is the radial coordinate, is a measure of the skyrmion radius, and is the wall width parameter. This profile was used to define the functional form of the polar angle in the model. An example of fits to the simulation data with the deformation ansatz is presented in Fig. S8.

We observe excellent quantitative agreement between the deformation described by the model and the micromagnetics simulations for all three components of the magnetization. From the fits, we can correlate the deformation amplitude, , with the background magnetization tilt, . This is shown in Fig. S9, where these parameters are shown as a function of applied in-plane magnetic field.

From this data, we deduce the empirical relationship, , where the tilt angle is expressed in degrees. These results also suggest that we can expect the deformation to be a linear function of the SOT strength, since the SOT plays a similar role to a magnetic field for a static profile.