Current-driven Dynamics of Magnetic Skyrmion in Chiral Ferromagnetic Film with Spatially Modulated Dzyaloshinskii-Moriya Interaction

TL;DR

This paper theoretically investigates how spatially modulated Dzyaloshinskii-Moriya interaction influences the current-driven dynamics of magnetic skyrmions in chiral ferromagnetic films, revealing controllable skyrmion behavior for spintronics applications.

Contribution

It introduces a theoretical framework for understanding skyrmion dynamics under spatially varying Dzyaloshinskii-Moriya interaction, enabling manipulation via material engineering.

Findings

Skyrmion size and velocity can be increased or decreased by different modulation patterns.

The Thiele equation effectively explains the observed phenomena.

Spatial modulation of Dzyaloshinskii-Moriya interaction offers a way to control skyrmion motion.

Abstract

The dynamics of magnetic skyrmion driven by spin-polarized current is theoretically studied in the chiral ferromagnetic film with spatially modulated Dzyaloshinskii-Moriya interaction. Three cases including linear, sinusoidal, and periodic rectangular modulations have been considered, where the increase, decrease, and the periodic modification of the size and velocity of the skyrmion have been observed in the microscopic simulations. These phenomena are well explained by the Thiele equation, where an effective force on the skyrmion is induced by the inhomogeneous Dzyaloshinskii-Moriya interaction. The results here suggest that the dynamics of skyrmion can be manipulated by artificially tuning the Dzyaloshinskii-Moriya interaction in chiral ferromagnetic film with material engineering methods, which will be useful to design skyrmion-based spintronics devices.

Figures

Equations

$$\mathcal{E}[\mathbf{m}]=A(\nabla \mathbf{m}{)}^2+{\mathcal{E}}_{DM}[\mathbf{m}]-K_u({\mathbf{e}}_z\cdot \mathbf{m}{)}^2.$$

$$\frac{d\mathbf{m}}{dt}=-{\gamma }_0\mathbf{m}\times {\mathbf{H}}_{eff}+\alpha \mathbf{m}\times \frac{d\mathbf{m}}{dt}-\beta \mathbf{m}\times (\mathbf{m}\times {\mathbf{m}}_p).$$

$$\hat{\mathcal{G}}\mathbf{v}-\alpha \hat{\mathcal{D}}\mathbf{v}-\mathbf{F}-{\mathbf{F}}^s=0.$$

$${\hat{\mathcal{G}}}_{{\xi }^{\prime }\xi }$$

$${\hat{\mathcal{D}}}_{{\xi }^{\prime }\xi }$$

$${\mathbf{F}}_{\xi }$$

$${\mathbf{F}}_{\xi }^s$$

$$v_x$$

$$v_y$$

$$\mathcal{D}=\pi {\int }_0^{\infty }[r({\Theta }^{\prime }{)}^2+\frac{1}{r}{\sin }^2\Theta ]dr.$$

$$\mathcal{F}=-\pi {\int }_0^{\infty }[r{\Theta }^{\prime }+\sin \Theta \cos \Theta ]dr,$$

$$v_x=\frac{\mathcal{G}\beta +\frac{2{\gamma }_0}{M_s}\alpha \nabla D\mathcal{D}}{{\mathcal{G}}^2+{\alpha }^2{\mathcal{D}}^2}\mathcal{F},v_y=\frac{\frac{2{\gamma }_0}{M_s}\mathcal{G}\nabla D-\alpha \beta \mathcal{D}}{{\mathcal{G}}^2+{\alpha }^2{\mathcal{D}}^2}\mathcal{F}.$$

$$\mathbf{m}\times \frac{d\mathbf{m}}{dt}$$

$$-{\nabla }_{{\mathbf{\xi }}^{\prime }}\mathbf{m}\cdot (\mathbf{m}\times \mathbf{v}\cdot {\nabla }_{\mathbf{\xi }}\mathbf{m})={\gamma }_0{\nabla }_{{\mathbf{\xi }}^{\prime }}\mathbf{m}\cdot {\mathbf{H}}_{eff}$$

$$+\alpha (\mathbf{v}\cdot {\nabla }_{\mathbf{\xi }}\mathbf{m})\cdot {\nabla }_{{\mathbf{\xi }}^{\prime }}\mathbf{m}+\beta {\nabla }_{{\mathbf{\xi }}^{\prime }}\mathbf{m}\cdot (\mathbf{m}\times {\mathbf{m}}_p).$$

$$\mathcal{Q}\mathbf{v}-\alpha \mathcal{D}\mathbf{v}-\mathbf{F}-{\mathbf{F}}^s=0.$$

$$\mathcal{Q}$$

$$\mathcal{D}$$

$$\mathbf{F}$$

$${\mathbf{F}}^s$$

$$\frac{\partial m_x}{\partial x}$$

$$\frac{\partial m_x}{\partial y}$$

$$\frac{\partial m_y}{\partial x}$$

$$\frac{\partial m_y}{\partial y}$$

$$\frac{\partial m_z}{\partial x}$$

$$\frac{\partial m_z}{\partial y}$$

$$\frac{\partial r}{\partial x}=\cos \phi ,\frac{\partial r}{\partial y}=\sin \phi ,\frac{\partial \phi }{\partial x}=-\frac{1}{r}\sin \phi ,\frac{\partial \phi }{\partial y}=\frac{1}{r}\cos \phi .$$

$$\mathcal{E}[\mathbf{m}]$$

$$D({\Theta }^{\prime }+\frac{1}{r}\sin \Theta \cos \Theta )\sin (\Phi -\phi )=0,$$

$${\Theta }^{\prime \prime }+\frac{1}{r}{\Theta }^{\prime }-(\frac{1}{2r^2}+\frac{K_u}{2A})\sin 2\Theta +\frac{D}{Ar}{\sin }^2\Theta \cos (\Phi -\phi )=0.$$

$$\mathcal{E}[\mathbf{m}]$$

$$D({\Theta }^{\prime }+\frac{1}{r}\sin \Theta \cos \Theta )\cos (\Phi -\phi )=0,$$

$${\Theta }^{\prime \prime }+\frac{1}{r}{\Theta }^{\prime }-(\frac{1}{2r^2}+\frac{K_u}{2A})\sin 2\Theta +\frac{D}{Ar}{\sin }^2\Theta \sin (\Phi -\phi )=0.$$

$${\mathcal{Q}}_{xy}=-{\mathcal{Q}}_{yx}={\int }_0^{\infty }dr{\int }_0^{2\pi }d\phi {\Theta }^{\prime }\sin \Theta =-4\pi ,$$

$$\frac{\partial \mathbf{m}}{\partial x}\cdot \frac{\partial \mathbf{m}}{\partial x}=({\Theta }^{\prime }{)}^2{\cos }^2\phi +\frac{1}{r^2}{\sin }^2\Theta {\sin }^2\phi ,$$

$$\frac{\partial \mathbf{m}}{\partial y}\cdot \frac{\partial \mathbf{m}}{\partial y}=({\Theta }^{\prime }{)}^2{\sin }^2\phi +\frac{1}{r^2}{\sin }^2\Theta {\cos }^2\phi ,$$

$$\frac{\partial \mathbf{m}}{\partial x}\cdot \frac{\partial \mathbf{m}}{\partial y}=(({\Theta }^{\prime }{)}^2-\frac{1}{r^2}{\sin }^2\Theta )\sin \phi \cos \phi .$$

$${\mathcal{D}}_{xx}={\mathcal{D}}_{yy}=\pi {\int }_0^{\infty }[r({\Theta }^{\prime }{)}^2+\frac{1}{r}{\sin }^2\Theta ]dr.$$

Full text

I Introduction

In recent years, the discovery of magnetic skyrmion in chiral ferromagnetic materials has triggered extensive research activities,JMMM1994 ; PRL2001 ; Nat2006 ; Sci2009 ; Nat2010 ; NM2010 ; NP2011 ; NP2012 ; PRL2009 ; PRL2011-1 ; PRL2011-2 ; Sci2013 ; NN2013-1 ; NN2013-2 ; PRL2013 ; NC2014 ; Sci2015 ; NP2017-1 ; NP2017-2 ; NN2017 mainly due to its potential applications in spintronics devices.NN2013-3 ; NRM1 ; NRM2 ; Proc ; NatEle2018 This type of topological defect is stabilized by the finite Dzyaloshinskii-Moriya interaction (DMI),Dzya ; Moriya which origins from the broken inversion symmetry in the materials. The dynamics of the skyrmion can be flexibly manipulated by spin-polarized current,NN2013-1 ; NN2013-2 and the trajectory of its center can be simply described by the Thiele equation.NN2013-3 It shows that a transverse component of the velocity will be generated by an effective force on the skyrmion, i.e. “the skyrmion Hall effect”,NN2013-2 which has been confirmed in recent experimental observations.NP2017-1 ; NP2017-2

Meanwhile, the physical factors influencing the DMI in chiral ferromagnetic materials have been carefully investigated. PRB2015 ; SR2015 ; PRB2016 ; SR2016 ; PRL2017-1 ; PRL2017-2 ; NL2018 ; PRB2018 ; Nanoscale2018 ; SR2018 For instance, it has been shown that the bulk DMI can be modified by the chemical composition,PRB2015 carrier density and the strain,SR2015 , band filling,PRL2017-2 etc. The interfacial DMI in the multilayer thin films, on the other hand, can be tuned by the thickness of certian layer,PRB2016 ; Nanoscale2018 ; SR2018 oxygen coverage,SR2016 ion irradiation,PRL2017-1 etc. Furthermore, the effect of electric field has also been studied,NL2018 ; SR2018 where 130% variation of DMI in Ta/FeCoB/TaOx trilayer has been demonstrated experimentally.NL2018 These developments thereby lay the foundation to control the dynamics of magnetic skyrmion by artificially modulate the DMI in chiral ferromagnetic materials.

By far, the current-driven dynamics of magnetic skyrmion has mainly been studied in chiral ferromagnetic materials with homogeneous DMI. In this paper, we consider a chiral ferromagnetic film with spatially modulated DMI, and show that the trajectory of the skyrmion will be drastically modified due to the inhomogeneous distribution of DMI, which can be well understood by analyzing the Thiele equation. The results here will be beneficial to design skyrmion-based devices by controlling the DMI with material engineering techniques.

The rest of paper is organized as follows. In Section II, we describe the theoretical model of the current-driven skyrmion with inhomogeneous DMI, and give the corresponding Thiele equation. In section III, we perform micromagnetic simulations with linearly or periodically modulated DMI, and compare the results with Thiele equation. The conclusive remarks of our work will be given in Section IV.

II Theoretical Model

As shown in Fig. 1, the model system is a chiral ferromagnetic film with length $L$, width $W$, and thickness $d$. The main interactions among the magnetic moments here include the Heisenberg exchange interaction, DMI, and the perpendicular magnetic anisotropy, which thus give the energy density of the film asNN2013-1 ; NN2013-2 ; NN2013-3

[TABLE]

Here, $\mathbf{m}(\mathbf{r})=\mathbf{M}(\mathbf{r})/M_{s}$ is the normalized magnetization distribution, and $M_{s}$ is the saturation magnetization; $A$ is the exchange stiffness; $K_{u}$ is the perpendicular magnetic anisotropy coefficient, and $\mathbf{e}_{z}$ is the anisotropic axis along $z$-direction. $\mathcal{E}_{DM}[\mathbf{m}]$ represents the energy density due to DMI interaction, where $\mathcal{E}_{DM}[\mathbf{m}]=D(m_{z}\nabla\cdot\mathbf{m}-\mathbf{m}\cdot\nabla m_{z})$ for interfacial DMI and $\mathcal{E}_{DM}[\mathbf{m}]=D\mathbf{m}\cdot\nabla\times\mathbf{m}$ for bulk DMI.NN2013-1 ; NN2013-2 ; NN2013-3 The coefficient $D$ determines the strength of DMI and is assumed to be spatially inhomogeneous.

A magnetic skyrmion will be prepared near the upper-left corner of the film, and a spin-polarized current perpendicular to the film will be applied. The amplitude, degree and direction of spin polarization of the applied current are denoted by $I$, $p$, $\mathbf{m}_{p}$, respectively. Then the magnetization dynamics of the film will be described by the Landau-Lifshitz-Gilbert-Slonczewski (LLGS) equation

[TABLE]

Here, $\gamma_{0}$ is the gyromagnetic ratio; $\mathbf{H}_{eff}=-\frac{1}{M_{s}}\nabla_{\mathbf{m}}\mathcal{E}[\mathbf{m}]$ is the effective magnetic field; $\alpha$ is the Gilbert damping coefficient. The last term in Eq. (2) is the Slonczewski spin transfer torque with the amplitude $\beta=\frac{\gamma_{0}\hbar Ip}{2dM_{s}e}$, where $\hbar$ and $e$ are the reduced Planck constant and elementary charge, respectively. The solution of Eq. (2) will give the current-driven dynamics of the skyrmion. Before going to the detailed micromagnetic simulations, we will first analyze the motion of skyrmion based on the Thiele equation. Book

If the spatial variation of the DMI is small enough in the scale of the skyrmion size, we can approximately describe the structure of the skyrmion with fixed DMI strength $D$ at its center. Its trajectory then can be described by the Thiele equation (Appendix A)

[TABLE]

Here, $\bm{v}$ is the drift velocity of the skyrmion center; $\hat{\mathcal{G}}$ and $\hat{\mathcal{D}}$ denote the gyromagnetic coupling tensor and dissipative force tensor; $\mathbf{F}$ and $\mathbf{F}^{s}$ are the effective forces acting on the skyrmion due to the inhomogeneous energy density and the spin transfer torque, respectively. They are explicitly related to the magnetization distribution $\mathbf{m}(\mathbf{r})$ as

[TABLE]

The variables $\xi,\xi^{\prime}$ here denote the spatial coordinates $x,y$. Then the velocity $\bm{v}$ of the skyrmion will be obtained from the Thiele equation as

[TABLE]

which then gives the Hall angle $\vartheta_{H}=\arctan(v_{y}/v_{x})$.

The interfacial and bulk DMI will result in Néel-type and Bloch-type skyrmion, respectively. In the polar coordinates taking the skyrmion center as the origin, the magnetization configuration of the skyrmion can be expressed asNN2013-3 $\mathbf{m}(r,\varphi)=(\sin\Theta(r)\cos\Phi(\varphi),\sin\Theta(r)\sin\Phi(\varphi),\cos\Theta(r))$, where $\Phi=\varphi+\gamma$. Depending on $D>0$ or $D<0$, the helicity $\gamma$ will be [math] or $\pi$ for the Néel-type skyrmion, and will be $\pi/2$ or $-\pi/2$ for the Bloch-type skyrmion (Appendix B). Further calculations (Appendix C) give the gyromagnetic coupling tensor \hat{\mathcal{G}}=\left(\begin{array}[]{cc}0&-\mathcal{G}\\ \mathcal{G}&0\end{array}\right) with $\mathcal{G}=4\pi$, and the dissipative force tensor \hat{\mathcal{D}}=\left(\begin{array}[]{cc}\mathcal{D}&0\\ 0&\mathcal{D}\end{array}\right), where

[TABLE]

The force $\mathbf{F}$ due to the inhomogeneous DMI strength $D$ is expressed as $\mathbf{F}=-\frac{2\gamma_{0}}{M_{s}}g_{f}{\mathcal{F}}\nabla D$. Here, we have defined the integral (Appendix C)

[TABLE]

and the coefficient $g_{f}=\cos\gamma$ for Néel-type skyrmion and $g_{f}=\sin\gamma$ for Bloch-type skyrmion.

When $\mathbf{m}_{p}=(\sin\theta\cos\phi,\sin\theta\sin\phi,\cos\theta)$, the force acting on the Néel-type skyrmion by the spin transfer torque is $\mathbf{F}^{s}=\cos\gamma(\sin\theta\sin\phi,-\sin\theta\cos\phi)\beta\mathcal{F}$, while the force on the Bloch-type skyrmion will be $\mathbf{F}^{s}=-\sin\gamma(\sin\theta\cos\phi,\sin\theta\sin\phi)\beta\mathcal{F}$.

The parameters in the Thiele equation (3) can be further calculated if the radial profile of the skyrmion $\Theta(r)$ is known, which then can predict the velocity of the skyrmion based on Eq. (8) and (9). Here, we take the ansatz $\Theta(r)=2\arctan[\frac{\sinh(R/w)}{\sinh(r/w)}]$ proposed recently,CP2018 with the skyrmion radius $R=\pi D\sqrt{\frac{A}{16AK_{u}^{2}-\pi^{2}D^{2}K_{u}}}$ and the skyrmion width $w=\frac{\pi D}{4K_{u}}$.CP2018 In Fig. 2, $R$, $w$, $\mathcal{D}$, $\mathcal{F}$ are shown as the functions of DMI strength $D$ by fixing the material parameters $A=15$ pJ/m and $K_{u}=0.8$ MJ/m3. All of them will monotonically increase when $D$ varies from $3$ mJ/m2 to $4$ mJ/m2.

III Micromagnetic Simulations

In this section, we study the current-driven dynamics of the magnetic skyrmion in the film by simulating the LLGS equation (2), and the simulated results are compared with the theoretical predictions from Thiele equation. The simulation parameters are set as $L=500$ nm, $W=100$ nm, $d=0.4$ nm, $M_{s}=5.8\times 10^{5}$ A/m, $A=15$ pJ/m, $K_{u}=0.8$ MJ/m${}^{3},\alpha=0.05$. The current amplitude is taken as $I=3\times 10^{11}$ A/m, with the degree of spin polarization $p=0.4$. We restrict the simulations here to the case of Néel-type skyrmion, and set the direction of polarization $\mathbf{m}_{p}=(-1,0,0)$. The film is discretized into grids with size $1.0\times 1.0\times 0.4$ nm3, and the time step during the simulations is taken as $1$ fs for convergence. Three cases will be considered for the spatial modulation of DMI strength $D$ below.

III.1 Linear Modulation

We first consider that the DMI strength is linearly modulated along the $x$-direction. In Fig. 3(a), the trajectories of the skyrmion are presented when $D$ increases from 3 to 4 mJ/m2 and decreases from 4 to 3 mJ/m2, respectively. For comparison, the trajectory for homogeneous $D=3.5$ mJ/m2 is also shown. The results reveal that the Hall angle $\vartheta_{H}$ can be controlled by introducing the gradient of $D$. Besides, the size of the skyrmion will also vary when $D$ is changed during the motion. In Fig. 2(b) and (c), the radius $R$ and width $w$ of the skyrmion are obtained by fitting the simulation results to the ansatz $\Theta(r)=2\arctan[\frac{\sinh(R/w)}{\sinh(r/w)}]$. We see that the fitted $R$ is generally larger than the theoretical value, which can be about $2$ nm for large $D$. On the other hand, the fitted $w$ is close to the theoretical values with a small difference about $0.3$ nm. Therefore, we can take the proposed model in Ref. CP2018, to approximately estimate the profile of the skyrmion.

The velocity of the skyrmion is further determined from the trajectory of the skyrmion, as shown in Fig. 3(d) and (e). As $D$ is linearly increased from $3$ mJ/m2 to $4$ mJ/m2, both components of the velocity $v_{x}$ and $v_{y}$ will become larger. Based on the results in Section II, the force due to the linearly modulated $D$ is along the $x$-direction and is given as $\mathbf{F}_{x}=-\frac{2\gamma_{0}}{M_{s}}\mathcal{F}\nabla D$, while the force due to the spin transfer torque is along the $y$-direction with the value $\mathbf{F}_{y}^{s}=\beta\mathcal{F}$. The Thiele equation will then give the velocity

[TABLE]

Fig. 2 has shown that $\mathcal{F}$ will quickly increase from $40$ to $90$ nm during the motion, thus both of $v_{x}$ and $v_{y}$ will increase along with the DMI strength $D$. Meanwhile, $\mathcal{D}$ will also slowly increase from $15$ to $18$ (see Fig. 2(c)), which thus explains the slight increase of the magnitude of Hall angle $\vartheta_{H}=\arctan(\frac{\frac{2\gamma_{0}}{M_{s}}\mathcal{G}\nabla D-\alpha\beta\mathcal{D}}{\mathcal{G}\beta+\frac{2\gamma_{0}}{M_{s}}\alpha\nabla D\mathcal{D}})$. Similarly, when $D$ decreases from $4$ mJ/m2 to $3$ mJ/m2, the velocity $v_{x}$ and $v_{y}$ will both decrease, while the magnitude of the Hall angle $\vartheta_{H}$ will continuously decrease. Besides, the magnitude of Hall angle will be small for positive $\nabla D$ and be large for negative $\nabla D$ at fixed $D$, as shown in Fig. 3(a) and (f). Fig. 3(d)(e)(f) suggest that the Thiele equation can correctly capture the trend of the velocity and Hall angle with linearly modulated $D$, although some quantitative differences exist.

III.2 Sinusoidal Modulation

We further study the current-driven dynamics of skyrmion with a sinusoidal modulation of the DMI strength, where $D=3.5+0.5\sin(2\pi x/100)$ mJ/m2 with $x\in[0,500]$ nm. Fig. 4(a) shows that the trajectory of the skyrmion will be periodically distorted by the modulated DMI. Specially, the skyrmion will be pushed away but can always move back to the trajectory if $D$ is set as $3.5$ mJ/m2. The behavior can be understood from the periodic modulation of the Hall angle, the magnitude of which will become large for negative $\nabla D$ and vice versa. The turning points $x_{0}$ of the trajectory will then be determined by the condition $\nabla D=0$, which gives $x_{0}=25,75,\cdots,475$ nm. The skyrmion size will also be periodically modulated as expected. In Fig. 4(b) and (c), the radius $R$ and width $w$ of the skyrmion obtained by fitting the simulation results are compared to the theoretical values. The size of the skyrmion will expand and shrink periodically together with the DMI strength. In spite of quantitative differences, the theoretical calculations do can predict the variation of the skyrmion size satisfactorily.

The velocity $\mathbf{v}$ and Hall angle $\vartheta_{H}$ have also been obtained from the micromagnetic simulations and calculated from the Thiele equation separately. As shown in Fig. 4(d)(e)(f), $\mathbf{v}$ and $\vartheta_{H}$ will oscillate around the reference values obtained for fixed $D=3.5$ mJ/m2. Besides, $v_{x}$ will achieve it maximal value when $D$ is largest and vice versa, as determined by the factor $\mathcal{F}$ in Thiele equation. However, $v_{y}$ in this case can be either positive or negative, since $\nabla D$ here is large enough to be dominant. In the same reason, $\vartheta_{H}$ can also be positive or negative, which results in the “push-back” trajectory given in Fig. 4(a). Moreover, the Hall angle with maximally negative value will be achieved when $\nabla D$ is negatively largest at $x=50,150,250,350,450$ nm, while the Hall angle with maximally positive value will be achieved at $x=100,200,300,400$ nm.

III.3 Periodic Rectangular Modulation

Finally, we consider the effect of the periodic rectangular modulation of DMI on the dynamics of skyrmion. In Fig. 5(a), the trajectory of skyrmion has similar “push-back” feature, where $D$ takes the values $4$ mJ/m2 or $3$ mJ/m2 alternately for each $50$ nm along $x$-direction. Unlike the sinusoidal modulation case, the skyrmion size will be suddenly changed because of the abrupt increase or decrease of DMI; otherwise, the radius $R$ and width $w$ will be nearly the same in the regions with constant $D$ (see Fig. 5(b)(c)). Once again, the agreement between the theoretical and simulated results will be better when $D$ is small.

We have assumed that the modulation of DMI should be slow in the scale of the skyrmion size in order to get the force $\mathbf{F}$ in Section II. This assumption will break down for the case here, since $D$ is a step function of the location and $\nabla D$ will give the $\delta$-pulsive force. Considering the finite size of the skyrmion, we exploit the Gaussian distribution $\frac{1}{\sqrt{\pi}\sigma}e^{-\frac{x^{2}}{\sigma^{2}}}$ to replace the $\delta(x)$ function, where $\sigma$ defines the broadening degree of the force $\mathbf{F}$. As shown in Fig. 5(d)(e)(f), the velocity $\mathbf{v}$ and Hall angle $\vartheta_{H}$ of the skyrmion are calculated from the Thiele equation by setting $\sigma=15$ nm, which can reasonably explain the simulation results. Not surprisingly, $\mathbf{v}$ and $\vartheta_{H}$ will repeatedly modulated, and the direction of $v_{y}$ will be changed to give the “push-back” motion in Fig. 5(a).

IV Conclusion

In conclusion, we have investigated the current-driven dynamics of magnetic skyrmion in chiral ferromagnetic film with spatially modulated DMI based on the Thiele equation and the micromagnetic simulations. We show that the trajectory and size of skyrmion can be manipulated by tuning the Dzyaloshinskii-Moriya interaction. Specially, the motion of the skyrmion can be sped up or slowed down by the gradient of DMI, and periodic DMI will cause the “push-back” trajectory of the skyrmion. These phenomena could be realized by the material engineering technologies nowadays, e.g. applying a spatially modulated gate voltage. The results here will not only advance our understanding of the rich dynamics of magnetic skyrmion further, but also can be useful to develop novel skyrmion-based spintronics devices.

Acknowledgements.

This work has been supported by NSFC Project No. 61674083 and No. 11604162 and by the Fundamental Research Funds for the Central Universities, Nankai University (Grant No. 7540).

Appendix A Derivation of Thiele Equation

Here, we give the derivation of the Thiele equation from the LLGS equation (2). Taking the cross product of $\mathbf{m}$ and Eq. (2), we will get

[TABLE]

which can further reduce to

[TABLE]

Here, we have applied the rule of vector triple product $\mathbf{a}\times(\mathbf{b}\times\mathbf{c})=\mathbf{b}(\mathbf{a}\cdot\mathbf{c})-\mathbf{c}(\mathbf{a}\cdot\mathbf{b})$ and the relations $\mathbf{m}\cdot\mathbf{m}=1,\mathbf{m}\cdot\frac{d\mathbf{m}}{dt}=0$. Assuming that the magnetization distribution of the skyrmion is not deformed during its motion, which is given as $\mathbf{m}(\bm{\xi})$ with the local coordinates $\bm{\xi}$ relative to its center $\mathbf{R}(t)$, then the magnetization distribution of the whole film in the lab coordinates will be $\mathbf{m}(\mathbf{r}-\mathbf{R}(t))=\mathbf{m}(\bm{\xi})$, and we will have $\frac{d\mathbf{m}}{dt}=-\bm{v}\cdot\nabla_{\bm{\xi}}\mathbf{m}$ with the skyrmion velocity $\bm{v}=\frac{d\mathbf{R}}{dt}$. Substituting $\frac{d\mathbf{m}}{dt}$ into Eq. (13) and taking the dot product with $\nabla_{\bm{\xi}^{\prime}}\mathbf{m}$, we further get

[TABLE]

After integrating over the whole film, Eq. (14) will reduce to the Thiele equation as

[TABLE]

Here, $\mathcal{Q}$ is the gyromagnetic coupling tensor and $\mathcal{D}$ is the dissipative force tensor; $\mathbf{F}$ and $\mathbf{F}^{s}$ are the forces on the skyrmion due to the inhomogeneous energy density $\mathcal{E}[\mathbf{m}]$ and spin transfer torque, respectively. They are defined as the functions of the magnetization configuration $\mathbf{m}(\mathbf{r})$ as

[TABLE]

Appendix B Magnetic Skyrmion Configuration

In the polar coordinates, the magnetization configuration the skyrmion is expressed by $\mathbf{m}(\mathbf{r})=(\sin\Theta(r)\cos\Phi(\varphi),\sin\Theta(r)\sin\Phi(\varphi),\cos\Theta(r))$. The components of $\nabla_{\bm{\xi}}\mathbf{m}$ are calculated as

[TABLE]

Here, we have utilized the relations

[TABLE]

Then the energy density $\mathcal{E}[\mathbf{m}]$ with the interfacial DMI reads

[TABLE]

The stable magnetization configuration is determined by minimizing the total energy $E[(31Dec2018)\Theta,\Phi]=\int d\mathbf{r}\mathcal{E}[\mathbf{m}]$ with the variational principle $\delta E=0$, which results in the following equations

[TABLE]

Therefore, the helicity $\gamma=\Phi-\varphi$ should be [math] or $\pi$, namely, Néel-type skyrmion. Besides, considering that $\Theta|_{r=0}=\pi$ and $\Theta|_{r\rightarrow\infty}=0$, the energy due to DMI will be minimized if $\gamma=0$ for $D>0$ or $\gamma=\pi$ for $D<0$.

Similarly, the energy density $\mathcal{E}[\mathbf{m}]$ with bulk DMI is

[TABLE]

The variational equation $\delta\int d\mathbf{r}\mathcal{E}[\mathbf{m}]=0$ then results in the following equations

[TABLE]

In this case, the skyrmion will be Bloch-type with the helicity $\gamma=\pi/2$ for $D>0$ or $\gamma=-\pi/2$ for $D<0$.

Appendix C Expressions of $\mathcal{Q},\mathcal{D},\mathbf{F},\mathbf{F}^{s}$

We now determine the expressions of $\mathcal{Q},\mathcal{D},\mathbf{F},\mathbf{F}^{s}$ for a magnetic skyrmion in terms of the polar coordinates. First, It is easy to see that the diagonal matrix elements $\mathcal{Q}_{xx}=\mathcal{Q}_{yy}=0$. While for its non-diagonal matrix element, straightforward calculations give $\frac{\partial\mathbf{m}}{\partial x}\times\frac{\partial\mathbf{m}}{\partial y}=\frac{1}{r}\Theta^{\prime}\sin\Theta\mathbf{m}$, therefore we can get

[TABLE]

which is the well-known skyrmion number with a $-4\pi$ multiplier.NN2013-3

In order to calculate $\mathcal{D}$, we first have

[TABLE]

Therefore, the non-diagonal matrix elements $\mathcal{D}_{xy}=\mathcal{D}_{yx}=0$ due to the integral over $\varphi$ from [math] to $2\pi$. The diagonal matrix elements of $\mathcal{D}$ will be

[TABLE]

The force $\mathbf{F}$ can be calculated from the energy density $\mathcal{E}[\mathbf{m}]$ given in Appendix B. Since we have assumed that the skyrmion is not deformed, the only non-zero contribution to the force $\mathbf{F}$ will come from the inhomogeneous $D$, which is a function of the skyrmion center $\mathbf{R}$. Since $\nabla_{\bm{\xi}}=-\nabla_{\mathbf{R}}$, we will have

[TABLE]

Here, the factor $g_{f}=\cos\gamma$ for Néel-type skyrmion and $g_{f}=\sin\gamma$ for Bloch-type skyrmion.

Finally, we calculate the force $\mathbf{F}^{s}$ due to the spin transfer torque. We first have

[TABLE]

[TABLE]

For Néel-type skyrmion where $\Phi=\varphi$ or $\Phi=\varphi+\pi$, only the $y$-component of $\int d\mathbf{r}(\partial_{x}\mathbf{m}\times\mathbf{m})$ and the $x$-component of $\int d\mathbf{r}(\partial_{y}\mathbf{m}\times\mathbf{m})$ are not zero after integral over $\varphi$. While for Bloch-type skyrmion where $\Phi=\varphi\pm\frac{\pi}{2}$, only the $x$-component of $\int d\mathbf{r}(\partial_{x}\mathbf{m}\times\mathbf{m})$ and the $y$-component of $\int d\mathbf{r}(\partial_{y}\mathbf{m}\times\mathbf{m})$ are not zero after the integral. Therefore, if $\mathbf{m}_{p}=(\sin\theta\cos\phi,\sin\theta\sin\phi,\cos\theta)$, the force acting on the Néel-type skyrmion will be

[TABLE]

while the force acting on the Bloch-type skyrmion will be

[TABLE]

Here, we have introduced the notation

[TABLE]

which has already appeared in the expression for $\mathbf{F}$.