Dynamics and stability of skyrmions in a bent nano-beam

Magnetic skyrmions are attracting much attention due to their nontrivial topology and high mobility to electric current. Nevertheless, suppression of the skyrmion Hall effect and maintaining high velocity of skyrmions with low energy cost are two major challenges concerning skyrmion-based spintronic devices. Here we show theoretically that in a nano-beam suffering appropriate bending moment, both Bloch-type and Néel-type skyrmions move with a vanishing Hall angle under a current density smaller than that required when the bending is absent. Moreover, bending alone can be used to move skyrmions, whose velocity is solved analytically from the Thiele equation. Generally speaking, inhomogeneous elastic fields affect the stability and dynamics of skyrmions, where the local stability is dominantly determined by the local bulk stress. These findings throw new light on how to drive skyrmions in a straight line with lower energy cost, which is vital for utilizing skyrmions as information carriers.


Introduction
A magnetic skyrmion is a nanoscale vortex-like spin structure usually stabilized by the Dzyaloshinskii-Moriya interaction (DMI) [1,2]. Magnetic skyrmions and their crystalline form called skyrmion crystals (SkX), have been observed in bulk materials [3,4], thin films [5,6], nano-sized discs [7,8] and nanowires [9,10] of P2 1 3 − type chiral magnets such as MnSi, Fe 1−x Co x Si and Cu 2 OSeO 3 . The nontrivial topology of skyrmions gives rise to exotic properties such as topological Hall effect [11], skyrmion magnetic resonance [12] and thermally induced ratchet motion [13]. Besides the intriguing physical properties, skyrmions have been attracting great interest for their potential applications in memory technology and storage devices [14,15] due to their small size [14], topological stability [16] as well as ultrahigh mobility to electric currents [17][18][19]. Specifically, skyrmions can be manipulated by electric current [20][21][22][23] with a density 10 5 times smaller than that required to move other magnetic textures such as domain walls [19,24]. However, current-driven dynamics of skyrmions suffer from the skyrmion Hall effect (SkHE), that is, the skyrmion trajectories deviate from the driving current direction due to the Magnus force [25,26]. To suppress the SkHE, various approaches have been proposed, including the skyrmionium [27], magnetic bilayer-skyrmion [28], antiferromagnetic skyrmion [29] and modifying the magnetic structure [30]. Besides the SkHE, the current density required to maintain a relatively high speed of skyrmions [15] is still quite high. To facilitate the application of skyrmion-based devices, a major challenge is to suppress or control the effect of SkHE, while simultaneously reducing the current density to maintain high speed motion of skyrmions.
Strain engineering of skyrmions, that is, manipulation of skyrmions by mechanical loads, is a promising direction. To begin with, it has been proved that homogeneous elastic fields change the stability and field patterns of skyrmions significantly [31][32][33][34]. Moreover, effectiveness of gradient fields on driving skyrmion motion has been confirmed, including magnetic field gradient [35][36][37], temperature gradient [13,[38][39][40] and electric field gradient [41]. In this sense, it is of great interest to examine the possibility of skyrmion motion driven by elastic field gradient.
In this work, we study the stability and dynamics of Bloch-type and Néel-type skyrmions in bent nano-beams, i.e., nano-beams subjected to a stress gradient, based on a free-energy functional [42,43] which has proved its effectiveness by quantitatively explaining various magnetoelastic phenomena of skyrmions. We find that skyrmion motion can be directly driven by bending the nano-beam, leading to a current-free driving method. The bending-driven skyrmion moves toward an oblique direction with respect to the direction of stress gradient, and the analytical expression of its velocity is derived based on the Thiele equation. When both current and bending are applied, the skyrmion velocity can be regarded as the resultant velocity of the current-driven motion and the bending-driven motion. We find that by applying an appropriate bending moment, one can effectively suppress the SkHE while simultaneously accelerating the moving skyrmion. In this way, the required current density can be significantly reduced to attain the original skyrmion velocity when the skyrmion is driven by electric current alone. The dynamic study is based on the extended Landau-Lifshitz-Gilbert (LLG) equation including the Slonczewski-like spin-transfer torque (STT) [44], which has been widely used in micromagnetic simulations and spintronic research [20,25,28,[45][46][47].

Results
Based on the Ginzburg-Landau functional theory, the equilibrium state of a helimagnet nano-beam suffering pure bending with a moment M on the left and right ends (see, e.g., figure 1(a1)) is determined by the rescaled elastic Gibbs free-energy density functional: where m is the rescaled magnetization, x, y, z represent components of the rescaled Cartesian coordinates, T is the rescaled temperature, b is the rescaled magnetic field, k eff term represents the uniaxial anisotropy, g 0 includes the terms independent of m,g DM denotes the DMI and g ME denotes the magnetoelastic terms. For Bloch-type skyrmions stabilized in B20 compounds, we have g DMI = 2m · (∇ × m) and where M and h denote respectively the bending moment and the beam size along y axis, K * , L * 1 , L * 2 , L * Oi (i = 1, 2, 3) are the magnetoelastic parameters, ξ is the coefficient of the fourth-order Landau expansion term, G = D 2 DM 4A where D DM is the DMI constant and A is the exchange stiffness, and m i,j = ∂m i ∂r j with (r 1 , r 2 , r 3 ) = x, y, z .
For Néel-type skyrmions stabilized in thin films with C 6v symmetry, we have g DMI = 2[m 1 m 3,1 + m 2 m 3,2 − m 3 m 1,1 + m 2,2 ] and where Γ and Λ are the magnetoelastic parameters. Equations (2) and (3) are obtained by substituting the solution of stresses for a bent beam into the elastic Gibbs free-energy density [43,48]. For given M, b and T, we perform free-energy minimization to determine m with the nonlinear conjugate gradient method. Details of the derivation of elastic Gibbs free-energy density and the numerical method used are shown in the 'methods' section. For the skyrmion driven by the spin current perpendicular to the plane (CPP) in ferromagnetic nano-beam, the magnetic spin dynamics is governed by the extended LLG equation including the Slonczewski-like STT [44], given as: where h eff = − ∂ g ∂m is the rescaled effective field derived from equation (1), t is the rescaled time, α is the dimensionless damping constant, i s is the rescaled spin current density and p is the polarization direction of the spin current. We numerically solve the LLG equation with Runge-Kutta 4 algorithm detailed in the 'methods' section.

Stability of skyrmions in bent nano-beam
For the two types of skyrmion, we apply pure bending with a moment M on nano-beams of three different sizes as shown in figure 1. Figures 1(a1)-(a5) and (b1)-(b5) present respectively the variation of stability of Bloch-type and Néel-type skyrmions in nano-beams containing one row of skyrmion at T = 0.5, suffering increasing bending moment. When the bending moment exceeds a threshold value M c , instability of the skyrmions occurs, initiated from the lower half of the beam. We determine the occurrence of instability of skyrmions by calculating the skyrmion number Q = 1 4π m · ( ∂m ∂x × ∂m ∂y )dx dy [49] and M c corresponds to the strength of the bending moment where the Q reduces sharply to approximately zero. M c equals 3 × 10 −7 N m (figure 1(a5)) and 3 × 10 −5 N m (figure 1(b5)) respectively for the nano-beam containing one row of Bloch-type skyrmions and one row of Néel-type skyrmions. As the bending moment increases, no significant deformation is observed, yet the distribution of magnetization is asymmetric at the upper half and the lower half of the beam.
Effects of bending on the nano-beams containing multi-row skyrmions are revealed in figures 1(ci)-(fi) (i = 1, 2, 3, 4, 5). Similar phenomenon is revealed for Bloch-type and Néel-type skyrmions: as the bending moment increases, annihilation of skyrmions at the lower half of the beam initiates from the bottom, where the skyrmions evolve to the ferromagnetic state. As the bending moment further increases, instability of the skyrmions at the upper half of the beam starts from the top, where the skyrmions evolve to the in-plane single-Q (IPSQ) phase.

Bending-driven dynamics of skyrmion in nano-beam
Similar to section 2.1, we apply pure bending with different moment M on the ferromagnetic nano-beam without current applied. At the initial moment t = 0 ps, a single skyrmion are created at the center of the nano-beam. We study the temporal evolution of the magnetization structure. The moving direction of skyrmion for positive and negative bending are shown in figure 2(a). In the nano-beam suffering positive (resp., negative) bending moment, a single Bloch-type skyrmion can be driven by the elastic field toward the upper left (resp., lower right), while Néel-type skyrmion is driven in the opposite direction. Here we define the skyrmion Hall angle θ through tan θ = V y /V x , where V x and V y denote the initial velocity of skyrmion along x and y directions. Figures 2(b) and (c) show the initial velocity and Hall angle of the bending-driven skyrmion motion, indicating that for both types of skyrmion, the moving direction is mainly independent of the absolute value of the bending moment, while the initial velocity of skyrmion increases with the bending strength.
For both types of skyrmion, the motion of a rigid magnetic skyrmion driven by bending in a ferromagnetic nano-beam can be expressed as follows based on the Thiele equation: for Bloch-type skyrmion and C = 3ξγR 2G 2 h 3 (RΛ − πΓ) for Néel-type skyrmion, where Q is the skyrmion number, M, R and γ denote respectively the bending moment, skyrmion radius and gyromagnetic ratio. D is the integral 1 4π ∂m ∂x · ∂m ∂x dx dy = 1 4π ∂m ∂y · ∂m ∂y dx dy on the area occupied by a single skyrmion. Detailed deduction of equation (5) is given in 'methods' section. Thus, we can deduce the bending-induced velocity of a single skyrmion: Here we have Q equals −1, D is positive, C is negative for Bloch-type skyrmion and positive for Néel-type skyrmion. Therefore, the sign of bending moment M determines the motion direction. These formulations explain well the modulus and direction of velocity for different bending moment presented in figures 2(b) and (c). Due to existence of the boundary and possible deformation of skyrmions, a slight deviation between the results from the simulations and the results from the Thiele equation exists.
Bending-driven skyrmion motion may be difficult to be applied in the design of skyrmion racetrack memory due to the significant Hall angle. However, previous researches reveal that skyrmion ratchet effects, which can be generated by alternating currents (AC) in a racetrack with a broken symmetry [50][51][52] or by asymmetric 2D periodic pinning [53], may lead to net unidirectional transport based on the SkHE. Similarly, alternating bending could also be used to drive the ratcheting motion as shown in figure 2(d), which provides new ideas for bending-driven skyrmion ratchet device.

Current-driven dynamics of skyrmion in bent nano-beam
Here we study the skyrmion dynamics when both bending and current are presented. Spin current perpendicular to the plane induced by the charge current in the heavy-metal layer with various current density is applied. For the Bloch-type skyrmion, charge current is applied along the short side of the nano-beam in heavy-metal layer as shown in figure 3(a), while for the Néel-type skyrmion, charge current is applied along the long side of the nano-beam as shown in figure 3(e).
For the stress-free nano-beam (M = 0 in figures 3(a) and (e)), the skyrmion exhibits an oblique trajectory as a result of the SkHE. When bending is applied to the nano-beam, it is found that the skyrmion Hall angle can be reduced. Figure 3(a) presents the trajectories from t = 0 ps to t = 2219.11 ps of a Bloch-type skyrmion driven by the current with a density of 9.48 × 10 11 A m −2 when bending with various bending moment is applied, and figure 3(e) presents the trajectories from t = 0 ps to t = 518.96 ps of a Néel-type skyrmion driven by the current with a density of 1.30 × 10 11 A m −2 . The skyrmion drawn is located at its initial position (t = 0 ps). From the trajectories, we find that the skyrmion Hall angle decreases as the bending moment increases and can even reach a negative value (e.g., M = 5 × 10 −6 N m in figure 3(a) and M = 3 × 10 −4 N m in figure 3(e)). Figures 3(b) and (f) present the skyrmion Hall angle as functions of bending moment for various current density. We notice that the skyrmion in a bent nano-beam can move straightly along the long side of the nano-beam instead of bending away in a stress-free nano-beam. Here we define the critical bending moment M c (J) as the value of the bending moment at which the skyrmion Hall angle equals zero, which is a function of the current density J. Figures 3(c) and (g) indicate that the strength of critical bending moment increases linearly with the current density. Figures 3(d) and (h) reveal the skyrmion velocity as functions of current density for various bending moment. It is generally known that one can augment the skyrmion velocity by applying spin-polarized current with a higher current density. Here we find that bending the nano-beam is another efficient approach to accelerate the moving skyrmion. As shown in figures 3(d) and (h), skyrmion velocity in the bent nano-beam under critical bending moment, which eliminates exactly the SkHE (the skyrmion Hall angle equals zero), is evidently higher than that in the stress-free nano-beam at one certain value of current density. To be more precise, for a Bloch-type skyrmion, current with density J = 1 × 10 12 A m −2 should be applied to reach a velocity of 20 m s −1 in a stress-free nano-beam. However, in the bent nano-beam under critical bending moment, current with J = 7 × 10 11 A m −2 can drive the skyrmion to the same velocity. This result reveals that by bending the nano-beam, the required current density to achieve the same speed of skyrmion can be significantly decreased. Hence, one can realize the suppression of SkHE and the acceleration of skyrmion at the same time by applying an appropriate bending moment.

Localized instability of skyrmions induced by inhomogeneous elastic fields
For the bending condition studied, we observe localized instability of skyrmions due to presence of inhomogeneous elastic fields. From equations (2) and (3), the influence of elastic fields to the magnetization is induced by three types of magnetoelastic terms: the term with a coefficient K * , the terms with coefficients L * i in equation (2) or Λ in equation (3), the terms with coefficients L * Oi in equation (2) or Γ in equation (3). Concerning the localized instability, the K * term has a dominating effect. To be more specific, the K * term and the L * i or Λ terms both affect the stability of skyrmions by renormalizing the coefficient of the second-order Landau expansion term, while we generally have K * L * i andK * Λ. The L * Oi or Γ terms, on the other hand, lead to DMI anisotropy, which mainly causes deformation and size change of the skyrmions. In this case, it will be convenient to consider only the K * term in the magnetoelastic interaction when studying the localized stability of skyrmions, for which we will have the condition T + 3 K * σ m (r) = 1 to determine if a skyrmion at location r is stable or not. Here σ m (r) = 1 3 σ ii (r) denotes the bulk stress at location r. At a given rescaled temperature T, a critical value of bulk stress σ mc can be derived as σ mc = 1− T 3 K * . The condition σ m (r) = σ mc roughly determines the phase boundary between the skyrmion phase and the ferromagnetic phase in a material suffering inhomogeneous elastic fields, where the skyrmion phase is located in the area σ m (r) > σ mc and the ferromagnetic phase is located in the area σ m (r) < σ mc . In figure 1, we have plotted the curve σ m (r) = σ mc with a black line, and it is shown that this line describes very well the phase boundary between the skyrmion phase and the ferromagnetic phase. This explanation is universal for both Bloch-type and Néel-type skyrmion suffering inhomogeneous elastic fields.

Variation of skyrmion Hall angle and skyrmion velocity induced by inhomogeneous elastic fields
To explain the bending-induced inhibition of SkHE and the acceleration of the moving skyrmion, we can regard the skyrmion velocity as the resultant vector of current-induced velocity and bending-induced velocity. Thiele equation for the current-bending-induced motion is expressed as follows: with i x and i y represent the two components of rescaled charge current, I r i r j = 1 Then equation (7) can be simplified as for both Bloch-type and Néel-type skyrmion in the cases studied. Therefore, bending represented by the term CM in equation (8) can be regarded as current applied in the direction perpendicular to the current i s , which explains the suppression of SkHE. To eliminate the SkHE which means V y = 0, the bending with bending moment M c (i s ) = i s IQ αCD should be applied, corresponding to the linearity in figures 3(c) and (g) since i s is proportional to the current density. In addition, at the same current density, V x equals −i s I αD in a bent nano-beam with bending moment exactly counteracting the SkHE, while V x equals −i s I αD+ Q 2 αD in a stress-free nano-beam. Since Q 2 αD is positive, the moving skyrmion accelerates due to the bending. The strain induced by the critical bending considered in this work is acceptable in most of the cases studied. For example, for the critical bending corresponding to the movement in a straight line in figures 3(a) and (e), the maximum local strain of the material equals 4.7% and 2.8% respectively, within the scope of the misfit strains that can be applied by epitaxial growth, and also below the ultimate strain values of most high-strength reinforcing steels. Moreover, if we apply a current with lower density, the critical bending reduces as shown in figures 3(c) and (g) so that the maximum strain also decreases.

Summary
In conclusion, we have revealed the bending-induced static and dynamic behaviors of Bloch-type and Néel-type skyrmions in nano-beams. The local bulk stress determines the local stability of a skyrmion against the ferromagnetic state when the material is suffering an arbitrary inhomogeneous elastic field. Furthermore, we find that skyrmion can be directly moved by bending the nano-beam without applying the spin current. For the current-driven motion of skyrmion, bending can be used to suppress the SkHE and simultaneously accelerate the moving skyrmion. In this way, we are able to control the moving direction of the skyrmion with less energy cost. Generally speaking, inhomogeneous elastic fields that break the rotational symmetry of a skyrmion drive its motion. Further investigation is needed to elucidate the effect of inhomogeneous elastic fields around inhomogeneities of materials or pinning sites on skyrmion dynamics.

Free-energy density functional
The equilibrium state of a nano-beam in an inhomogeneous stress field is determined by the rescaled Gibbs free-energy density functional, which can be derived from the rescaled Helmholtz free-energy density functional proposed by Hu [43]: g me (m, σ). (9) In equation (9), g 0 contains the terms of exchange energy density, the DMI, the Zeeman energy density, the uniaxial anisotropy energy density, the second-and fourth-order Landau expansion terms: where m is the rescaled magnetization, g DM represent the DMI, b is the rescaled magnetic field, T is the rescaled temperature and r = [r 1 , r 2 , r 3 ] T = [x, y, z] T represents the vector of rescaled Cartesian coordinates (in the 'methods' section, we sometimes replace x, y, z with r i (i = 1, 2, 3) for convenience). The rescaled form g 0 is derived from where and For Bloch-type skyrmion, g el denotes the elastic energy density of materials with cubic symmetry: where σ = [ σ 11 , σ 22 , σ 33 , σ 12 , σ 13 , σ 23 ] T = ξ G 2 σ is the rescaled stress and S is the inverse matrix of the rescaled stiffness matrix C of a cubic system: g me denotes the magnetoelastic energy density. For Bloch-type skyrmion, we have For pure bending, only σ 11 exists. Then equation (9) can be simplified as For Néel-type skyrmion, we have where Here, m 1 , m 2 and m 3 denote the components of the rescaled magnetization m. K, L 1 , L 2 , L 3 , L Oi (i = 1, 2, . . . , 6), λ n (n = 1, 2, ..., 6) and γ n (n = 1, 2, ..., 10) are the thermodynamics parameters [43]. The For Pt/Co/Ta multilayer film, the parameters used are listed as follows: We did not find the experimental value of K * for Pt/Co/Ta multilayer film in any previous works. Therefore, we take the same value as MnSi.

Nonlinear conjugate gradient algorithm
We perform nonlinear conjugate gradient algorithm [56,57] to minimize the free-energy of the system. To perform the calculation, we discretize the free-energy density functional shown in equation (9) with 2D grids in Cartesian coordinates. For a given external stress field, the energy density functional can be written as g = g(m 1 (x, y), m 2 (x, y), m 3 (x, y)).
We use the central difference method to calculate the value of partial derivatives: for k = 1, 2, 3. At each time, we minimize the local free-energy of one cell (i, j). The gradient −∇ m g indicates the direction of the steepest descent of g. Consider an initial local magnetization m 0 . Let Δm 0 = −∇ m g(m 0 ). By performing a line search in this direction, we solve α 0 = min which is a polynomial extremum problem. Let K eff . Now consider the nth step of the minimization with a local magnetization m n . Let Δm n = −∇ m g(m n ). With Polak-Ribiere model [56], we compute the parameter ξ PR n as We update the conjugate direction s n = Δm n + ξ PR n s n−1 , and we perform again the line research α n = min α ( g(m n + αs n )). Then m n+1 = m n + α n s n becomes the new local magnetization of this cell. We repeat the steps above until the free-energy density of the local cell converges. This process of energy minimization is applied to each cell of the grid repeatedly until the free-energy density of the whole system converges.

Landau-Lifshitz-Gilbert equation
The dimensionless LLG equation including the spin torque from the spin Hall effect can be written as: where h eff = − δ g δm is the rescaled effective field, γ is the gyromagnetic ratio, α is the damping coefficient originating from the spin relaxation, i s is the STT coefficient, i s is the field-like out-of-plane STT coefficient and p represents the electron polarization direction. We have i s = μ 0 e J|p| 2dM s with μ 0 the vacuum magnetic permittivity, d = 1 nm the film thickness of ferromagnet layer, J the current density and M s the saturation magnetization equaling 1.63 × 10 5 A m −1 for MnSi [43] and 6.97 × 10 5 A m −1 for Pt/Co/Ta multilayer film [55].
The field-like out-of-plane STT coefficient i s is set to zero. Substituting t = γt in equation (25), we have

Runge-Kutta 4 algorithm
In numerical analysis, the Runge-Kutta method, developed around 1900 by the German mathematicians Carl Runge and Wilhelm Kutta, are one of the most famous algorithms to solve the problems in form of In this explicit method, y = f (x, y) can be written as Y n+1 = Y n + hf (X n , Y n ) where h is the time step for the iteration. With the mean value theorem, there exist 0 < s < 1 so that Y n+1 = Y n + hf (X n + sh, Y (X n + sh)). K = f (X n + sh, Y (X n + sh)) is called the mean slope, and Runge-Kutta is an algorithm aiming at calculating K. The Runge-Kutta with a truncation error equaling O(h 5 ) is formulated as follows: In this way, we can evaluate Y n+1 from Y n .
Substituting equations (30) and (31) back in equation (29), we have Replacing m × ∂m ∂t in equation (26) For the temporal discretization, we set k as the index of the discretized time and h as the time step. Then equation (33) can be written as Then we are able to perform the Runge-Kutta 4 algorithm.

Thiele equation for bending-induced skyrmion motion
We commence with the LLG equation without current applied: Applying the operation of f · ∂m ∂r i (integration on the area occupied by a single skyrmion) on both side of equation (37), we have Assuming the magnetization structure of the skyrmion in rigid motion, we have ∂m ∂t = − j V j ∂m ∂r j with j = 1, 2, 3 indicating the three directions. We further assume that m 3 varies linearly in radial direction. We can therefore deduce: and with D = 1 4π ∂m ∂x · ∂m ∂x dx dy = 1 4π ∂m ∂y · ∂m ∂y dx dy. For the left side of equation (38), we have since δ g 0 δm · ∂m ∂r i = 0 due to the translational invariance. For Bloch-type skyrmion, the integral of equation (41) with C = 12ξ G 2 h 3 . For Néel-type skyrmion, the integral of equation (41) Substituting equations (39), (40) and (42) or equation (43) into equation (38), we get with C = 3ξγR 4G 2 h 3 L * 1 R − 2 L * 2 R + π L * O2 for Bloch-type skyrmion and C = 3ξγR 2G 2 h 3 (RΛ − πΓ) for Néel-type skyrmion.