Dynamic behavior of skyrmion collision: spiral and breath

A magnetic skyrmion is a particle-like topological soliton, which is an ideal candidate for developing high-density storage and logic devices due to its nonvolatility and tunability. In view of the particle motion characteristics of skyrmion, different skyrmions in a material inevitably interact in the form of short-range repulsion and long-range attraction. In this work, the dynamic characteristics of skyrmion collision in a ferromagnetic Co thin film are investigated by using micromagnetic simulations. It is found that the dynamic behavior of skyrmion after collision is highly dependent on the size of the strip, the initial velocity of skyrmion and magnetic damping constant. For the collision of two skyrmions, when the strip width exceeds the critical value, the skyrmions form a pair and rotate counterclockwise in the form of spiral and breath. It is interesting that the rotation and breath of skyrmions keep the same periodicity under the negligible damping, and the frequency increases with the increase of the initial velocity of skyrmion. Further, the collision of a system of three skyrmions reveals that they interact in pairs to form closed periodic trajectories. The results of the present work not only give an insight into the multi-skyrmion dynamics, but also provide guidance for the development of spintronic devices based on multi-skyrmion motion.

Without external fields, the dynamic behaviors of skyrmions are closely related to the number of skyrmions. For individual skyrmion, it is only need to consider its initial state and interaction with the boundary [38]. For multi-skyrmions, however, it is necessary to consider how skyrmions interact with each other. Recent theoretical and experimental studies have shown that there are short-range repulsions and long-range attractions between skyrmions [39][40][41][42][43][44]. Capic et al [45] demonstrated that static skyrmion pairs with the same chirality repelled each other and unfolded spirally. Reichhardt and Reichhardt [46] found that particles with the same Magnus force formed stable pairs, triplets and higher-order clusters, or exhibited chaotic motion.
In a system of multi-skyrmions, because of the particle motion properties, collision is common, especially the moving skyrmion colliding with the static ones. Iwasaki et al [47] investigated the elastic collisions between the emitted skyrmion and the static one without an applied external field. Interestingly, after separation, the skyrmions accelerate along the boundary, but not along the center, which may be the small size limits the full interaction between skyrmions. However, little is known about how strip size affects the collisions between two skyrmions, let alone collisions with different initial velocities, different magnetic damping constants, and three skyrmions.
Here, we investigate the collision dynamics between the moving skyrmion with different initial velocities and the static skyrmion in ferromagnetic Co thin film by using micromagnetic simulations. It is found that the strip width, the velocity of skyrmion and magnetic damping play an important role in the dynamic behavior of skyrmion collision. For the strip with small width, the skyrmions bounce off after collision due to the effect of boundary. For the strip with large width, in the absence of boundary effect, skyrmions form a pair and rotate counterclockwise in the form of spiral and breath, in which the rotation speed and breath frequency increase with the increase of the initial velocity. For the film with large magnetic damping, the rotating skyrmion pair gradually separates into isolated skyrmion because of the energy dissipation. In addition, the collision of three skyrmions tends to interact in pairs and forms the peanut-like trajectories. The present work gives an insight into the collision dynamics between multi-skyrmions under different conditions, which is beneficial to develop new magnetic memory and logic devices based on multi-skyrmions.

Methodology
In order to construct the dynamic model of skyrmions collision, micromagnetic simulations based on the MuMax3 program [48] were performed to investigate the dynamic properties. The magnetization dynamics of skyrmions follow the rule of the Landau-Lifshitz-Gilbert (LLG) equation with the spin transfer torque τ stt associated with the electric current along the surface as follows [47,49,50] where m represents the spatial distribution of the local normalized magnetization vector with a saturation magnetization M s . α is the Gilbert damping parameter, γ 0 is the gyromagnetic ratio, e is electronic charge, P is spin polarization rate of the current, j is spin-polarized current density, β is no-adiabatic coefficient, and µ B is Bohr magneton. The effective field H eff is defined by where µ 0 is the vacuum permeability. E tot is the total energy of the system, which can be written as a function of magnetization The magnetic exchange density term f exch is determined by the gradient of the local magnetization and the exchange coupling constant A ex , which describes the isotropic Heisenberg exchange interaction: The magnetostatic energy density f mag is given by A finite difference discretization allows the magnetostatic field to be evaluated as a (discrete) convolution of the magnetization with a demagnetizing kernelK, here, theK is Magnetostatic kernel [48].
The perpendicular magnetic anisotropy (PMA) energy density derived from the top and bottom surfaces of the ultrathin magnetic thin film is given by where K 1 and K 2 are the PMA coefficients. The interfacial DMI energy density is where D is the DMI constant, which describes the strength of interfacial DMI. The material parameters used in the present work are obtained from literatures [21,[51][52][53][54], which are listed as follows, M s = 0.58 × Here, for simplicity of the simulation, the spin polarization rate is assumed to be P = 1. The in-plane spin current density is uniformly set to j = −1 × 10 12 Am −2 for investigating the effects of size, damping and dynamic behavior of three skyrmions. In an ideal condition, in order to explore the effect of the emissive velocity of skyrmion, we vary the current density from j = −0.5 × 10 12 Am −2 to j = −8 × 10 12 Am −2 to obtain different initial velocities of the emitted skyrmion varying from 50.8 m s −2 to 774.3 m s −2 as shown in figure 5(d).
On account of the properties of particle motion, the skyrmion can be emitted from a region with a transverse spin current to another without external fields. Therefore, we carried out the collision between the emitted skyrmion and the stationary isolated skyrmion as shown in the figure 1, where the Co/Pt bilayer structure film was uniformly placed on the surface of the Si substrate. To be specific, two isolated skyrmions were created and stabilized on the Co film by a spin current perpendicular to the film. The skyrmion located in the green region is marked with S 1 , driven by the in-plane transverse current, and moved to the S 1 ′ in the gray region. The current pulse was withdrawn at the moment that the emitted skyrmion left the drive region and entered the collision region. The skyrmion without the applied field is marked with S 2 , interacting with S 1 ′ and moving to S 2 ′ . The dashed arrows in the figure 1 schematically indicate the possible trajectories of the skyrmions. The yellow region is strip edges or constructed defect area, used to constrain the motion of the skyrmions. The inset illustrates the magnetization distribution of the skyrmion, where the colors of arrows changed from blue to red indicate the distribution of m z from −1 to +1. The model size is 400nm × 60nm × 1nm, and the finite difference cell size is 1nm × 1nm × 1nm.
In the present model, the applied spin current just provides the initial velocity for S 1 . There are no external fields acted on the collision process of the skyrmions that we focus. Therefore, the torque term of τ stt in equation (1) in the part of Methodology in this paper can be ignored in the dynamic solution process after the spin current is withdrawn. By take the derivation of the LLG equation without τ stt , the standard Thiele equation can be obtained to describe the skyrmion motion process [55] where G is the gyrovector, D is the dissipative force tensor and v is the velocity of skyrmion. G = Ge z in which e z is the unit vector along z direction and G = −4π for the skyrmion. The nonzero components D ij of the tensor D are D xx = D yy = αD 0 . F ss is the repulsive force between skyrmions and can be described from the following equation [39,56]: , where K 1 is the modified Bessel function and H k is the perpendicular anisotropy field. This follows from the derivation of the interaction potential for the overlapping spin fields of the two bodies based on the Thiele model. The interaction potential U at different skyrmions distances r d and the rule of change are be discussed in supplementary figure 2 [57]. F edge is the repulsive force given by the strip edge. By the analysis of the Thiele equation, the forces acted on the skyrmions during the movement can be described, denoted with arrows in figure 1. Here, F g = G × v is the Magnus force, perpendicular to the velocity direction, and F d = D × v is the damping force, antiparallel to the velocity direction.
In order to accurately determine the position of the skyrmion in the ferromagnetic system, we introduce the skyrmion topological number Q and topological charge density q, which are expressed as [55] And the position of skyrmion R = (x c , y c ) can be calculated from the equations (10) and (11) by The trajectories of the skyrmions during motion can be described by the R = (x c , y c ) based on the congruent relationship between the central spin and the position of a skyrmion [58,59]. However, some errors occur when calculating trajectories for multi-skyrmions, we select the central part with the opposite direction (S z < 0) as the position of the skyrmions. In order to track the trajectory, multi-skyrmions centers are also pre-set to facilitate the calculation.

Size effects on collisions
Based on the above model, the collision of skyrmions can be simulated. Iwasaki et al [47] used an alternating spin current to drive the skyrmion with the initial velocity. However, the driving force of the alternating spin current exerted on the skyrmion is not enough to actuate it with the desired high velocity. Beyond that, the finite strip size plays a role in the collision and interaction process of the skyrmions. Here, we adopt the natural edge of the strip in the simulation model for simplicity. Therefore, due to the restriction of strip size, the size effects on the collision process of skyrmions with the higher initial velocity are analyzed in the following subsection.
In this section, the collision between two skyrmions with the Gilbert damping constant of α = 0.01 was carried out in a strip with the width of 60 nm, and the simulation results are shown in figure 2. The emissive region is on the left of the yellow dotted line in figure 2(c), with the location coordinate x ⩽ 100 nm. The other region, x > 100nm, is the collision region without external fields. The spin current density, J = 1e 12 A m − 2 , is applied on the emissive region. After S 1 was accelerated along the x-direction by the spin current, it moved to the collision region and collided with S 2 . During this process, the trajectories of S 1 and , the moving S 1 gradually got close to S 2 before T 1 . When interacting, S 2 moved toward the upper edge of strip. Meanwhile, S 1 changed the direction, moving away from S 2 along the lower edge of strip, indicating that S 1 and S 2 scattered after the collision.
In addition, the magnetization distributions of strip at several typical moments were snapshotted in detail in figure 2(c). These moments marked with T from T 0 to T 4 correspond to those in figures 2(a) and (b), respectively. T 0 = 0 ns represents the moment of initial state when the two skyrmions kept stable in the strip. After the application of the spin current, S 1 started to move and then interact with S 2 at the moment of T 1 = 0.25 ns. When T 2 = 0.6 ns, S 1 entered the collision region. Simultaneously, the spin current applied to drive S 1 was removed, making the entire strip in a state without any external fields. At T 3 = 1.1 ns, S 1 and S 2 were just out of interaction. Then, they would not interfere with each other and moved along the strip edge as shown in the figure 2(c) of T 4 = 2.5 ns. Continuous dynamic process of the skyrmions collision can be observed in supplementary movie 1. The distance between the two skyrmions is defined as r d , which can be used to determine whether the two skyrmions are separated since r d = 50 nm is the limitation of the overlapping distance between the spin fields of two skyrmion. The details can be obtained from supplementary figure 1.
In the narrow strip, the trajectories of the skyrmions always deviate from the center and move along the boundary. Obviously, the width of strip plays an important role in the interaction of the two skyrmions. Thus, the width of the strip was set to 80 nm and 100 nm, respectively, with the trajectories of the skyrmions recorded in the figures 3(a) and (b). Due to the weakening of the repulsion force applied to the skyrmion from the strip edge, the collision process of the skyrmions is different from that in a narrow strip: S 1 and S 2 rotated counterclockwise, forming an unfolded spiral trajectory. The continuous dynamic process to form the spiral trajectory can be observed in supplementary movie 2. The counterclockwise rotation is caused by the combined effect of the Magnus force F g and the repulsive force F ss between the skyrmions. In an ideal undamped system, the skyrmions keep the counterclockwise rotation all the time. However, with the damping loss, the velocity of skyrmions is attenuated according to the Thiele equation [55], resulting in a decrease in F g and F ss , which increases the distance r d between skyrmions. Since there are no external forces applied during the interactions of the collision, the velocity of the skyrmion can be derived from the Thiele equation of equation (9): . Here, F is the repulsive force between the skyrmions. The first term in the equation describes the skyrmion pair moving circularly, and the second dissipative term converts the motion of skyrmion from the circle to a spiral expansion.
As shown in figure 3(a), when r d reaches to about half of the strip width, the skyrmions will hit the strip edge and speed up [55]. In figure 3(b), when the width of strip increases to 100 nm, the motion trajectory of the skyrmions shows continuous rotation. Due to the damping loss, the distance r d becomes larger and larger, bringing about the ultimate separation between skyrmions after a long period of time. The existence of boundary conditions not only constrains the counterclockwise rotation of skyrmions, but also accelerates the separation of skyrmions. Figure 3(c) displays the variation of separation time with the strip width, indicating different processes of collision and movement under different widths of strips. As shown in figure 3(c), direct scattering occurs when a strip width increasing from 60 to 80 nm denoted by cyan region, leading to the separation of skyrmions after a short time of interaction. When a strip width reaching to 80 -150 nm denoted by yellow region, the trajectory of skyrmions is spiral and counterclockwise. Once hitting the strip edge, the skyrmion scatters and is separated from the other as shown in figure 3(a). During this period, it will take longer time to separate with the increasing width of strip. When it is greater than 150 nm denoted by white region, the skyrmions maintain continuous counterclockwise and spiral trajectory as shown in figure 3(b), regarded as the condition without any boundary constraint.

Damping effects on collisions
According to the equation (9), dissipative force tensor D derived from the damping has a great influence on the magnetization and the velocity of the skyrmions, bringing about the variation of trajectory in the collision. Thus, the magnetic damping effects on the skyrmions collision need to be further considered. In order to explore the effects of damping, the boundary effects should be eliminated. So, the width of strip is set as 400 nm. Different from infinite rotation and direct scattering, the motion of the skyrmions under the high damping is more concise. Figure 4(a) is snapshots of magnetization distribution of the skyrmion collision at the moments of T 0 -T 4 with α = 0.1. Here, T 0 -T 4 corresponds to the specific time of 0 ns, 0.4 ns, 0.6 ns, 15.5 ns, 500 ns, respectively. Combination with the trajectory depicted in figure 4(b), the movement process of skyrmions collision can be analyzed in detail. From T 0 to T 1 , S 1 , driven by the spin current, kept moving and interacted with S 2 . After the collision, the skyrmions counterclockwise rotated with the increasing radius of rotation. This is mainly because of the interaction, that is, the long-range attraction gradually dominates compared with the short-range repulsion, affecting the direction of skyrmions. Due to the excessive damping dissipation, the skyrmion velocity dropped rapidly after the withdrawal of spin current. Once reaching to the separation condition that the distance is so far that the interaction can be ignored, the skyrmions would stop rotating and then move separately. The continuous and dynamic motion of skyrmions can be observed in supplementary movie 3. Figure 4(c) is the variation of separation time with the damping constant α, manifesting that the increase of α accelerates the magnetization precession and sharply decreases the separation time of skyrmions.

Collisions under ideal conditions
To further explore the interaction between the two skyrmions, we carried out the collisions under ideal conditions, namely neglecting the effects of damping and boundary conditions. Thus, α = 0.001 and w = 400 nm are set in the following simulation. Due to the neglection of boundary and damping effects, the terms of F d and F edge in equation (9) can be ignored. As a result, the movement of skyrmions during the collision is determined by the combined effect of F g and F ss , showing the counterclockwise rotation in a spiral. As shown in figures 5(a)-(c), the cyan region represents the counterclockwise rotation of skyrmions in a spiral completely. The special moments are marked, of which 0.53 ns is the moment when S 1 left the emissive region. 0.63 ns is the moment when the two skyrmions were closest. 0.72 ns is the moment when the spin current was removed. Figure 5(a) presents the variation of F g = ( F x , F y ) acted on the skyrmions with time, in which the phase difference of the Magnus force F g of S 1 and S 2 is 180 • , indicating that S 1 and S 2 completely perform counterclockwise rotation. Because of F g = G × v, the variation of velocity directions of S 1 and S 2 is the same with that of the Magnus force as presented in figure 5(b).
When the damping effects can be ignored or at very low level, the sharp magnetization precession of skyrmions significantly changes the radius of skyrmions. During this process, the variation of radius with time is called as breathing mode [44][45][46][47]. As shown in figure 5(b). The rangeability and frequency of radius variation stay the same compared with S 1 and S 2 . Besides, the rotational period reflected by the linear velocities along the x-direction is consistent with that of the breathing mode as a result of the strong spin field between the skyrmions during the rotation. Furthermore, emissive velocity of the skyrmion also affects the movement and breathing mode. As shown in figure 5(d), when the emissive velocity is below 200 m s −1 , the rotational velocity and breathing frequency of the skyrmions increase dramatically. While above that, the rate of increase gradually slows down.
To further explain the interaction behavior of skyrmion collisions, we use the formulation of the domain wall energy proposed by Hu et al [60] where σ w is the local domain wall energy, U w is the total domain wall energy, A ex is the exchange coupling constant, D is the DMI constant, K is the PMA coefficient, R is the skrymion radius, w d is the wall width, and h is the thickness of film. According to equation (13), the domain wall energy is composed of the magnetic interface anisotropy energy, exchange energy and interfacial DMI energy. In this work, w d = 4 nm. Hence, the magnetization distribution and radius of the skyrmions are closely related to the domain wall energy. Figure 5(c) shows the variation of the domain wall energy with time. According to equation (14), due to the constancy of wall width w d and the thickness of h, the total domain wall energy U w is related to the individual skyrmion radius R [60]. Compared with the change of radius in figure 5(b), the domain wall energy keeps the same periodicity. Obviously, the velocity of S 1 originates from the applied spin current. The interaction with S 1 generates the variation of domain wall energy, contributing to the movement of S 2 . When counterclockwise rotating in a spiral, the two skyrmions are in a state of dynamic equilibrium, the conversion of kinetic and domain wall energy contributes to the consistent frequency of breathing mode and rotational velocity.

Three skyrmions collisions
Actually, compared with a system of two skyrmions, multi-skyrmions are more likely to be generated in the experiment because of the complicated external condition and material structure. Thus, it is essential to study the collision of multi-skyrmions for practical applications. We implemented the simulation with neglecting the effects of boundary and setting the damping constant α to 0.001, regulating a skyrmion with an initial velocity to hit the other two paratactic and stationary skyrmions, and investigated the interaction among multi-skyrmions. Figure 6(a) shows the snapshots of magnetization distribution of ferromagnetic films. T 0 = 0 ns represents the moment of initial magnetization state. At T 1 = 0.5 ns, S 1 moved to the collision zone and interacted with S 2 and S 3 . It can be seen that the radii of S 1 and S 2 became smaller, while S 3 remained unchanged. At the moment of T 2 = 1 ns, S 1 stopped interacting with S 2 , and started to interact with S 3 . At T 3 = 1.5 ns, S 3 kept interacting with S 1 , then finished the interaction at T 4 = 2 ns. From T 1 to T 4 , S 1,2,3 realizes the position replacement of S 3,1,2 , and the cycle process is repeated thereafter. Corresponding trajectories of the three skyrmions are shown in figure 6(b), where we generally observe closed periodic trajectories, like peanuts. Due to inertia and damping, the periodically changing positions of the skyrmions do not match perfectly [46], but the periodic interaction behavior can be obtained. After 30 ns, the peanut-like trajectories rotate at a slight angle step by step, forming a flower-like pattern depicted in figure 6(c). The long-time and dynamic motion of three skyrmions collision can be acquired in supplementary figure 3 and movie 4. As described, the investigation in the collision of a system of three skyrmions contributes to in-depth understanding the interaction and getting the movement rules of multi-skyrmions, prospectively realizing the quantitative regulation of multi-skyrmions in ferromagnetic materials.

Conclusion
In summary, the dynamic behavior of skyrmion collision in a ferromagnetic Co thin film has been systematically studied by utilizing micromagnetic simulations. The simulation results show that the size of the strip film, the magnetic damping constant and the initial velocity of the skyrmion have a great influence on the dynamic behaviors of skyrmion collision. Different movement modes of skyrmions during collisions are predicted as the strip width increases. First, the skyrmions bounce off once colliding due to the effect of boundary in the strip with small width. Second, due to the absence of boundary effect in the strip with large width, skyrmions form a pair and rotate counterclockwise in the form of spiral and breath. It is interesting that the rotation and breath of skyrmions keep the same periodicity under the negligible damping, in which the frequency increases with the increase of the initial velocity of skyrmion. When the magnetic damping becomes large, the rotating skyrmion pair gradually separates into isolated skyrmion because of the energy dissipation. While, rather different from two skyrmions, a novel collision mode, interacting in pairs, is predicted in the system of three skyrmions. These new findings and mechanisms on dynamics of skyrmion collision will provide theoretical guidance for experiments, and hopefully promote the development of magnetic storage and spintronic devices.

Data availability statement
All data that support the findings of this study are included within the article (and any supplementary files).