Phys. Rev. B 105, 214415 (2022).

# Nonhermitian adiabatic perturbation theory of topological quantization of
the average velocity of a magnetic skyrmion under thermal fluctuations

Shan-Chang Tang  Department of Physics & State Key Laboratory of Surface
Physics, Fudan University,
Shanghai 200433, China  Yu Shi<EMAIL_ADDRESS>Department of Physics &
State Key Laboratory of Surface Physics, Fudan University,
Shanghai 200433, China

###### Abstract

We study the two-dimensional motion of a magnetic skyrmion driven by a
ratchetlike polarized electric current that is periodic in both space and
time. Some general cases are considered, in each of which, in the low
temperature and adiabatic limit, regardless of the details of the driving
current, the time and statistical average velocity along any direction is
topologically quantized as a Chern number, multiplied by a basic unit. We make
two approaches, one based on identifying the drift direction, and the other
based on the nonhermitian adiabatic perturbation theory developed for the
Fokker-Planck operator. Both approach applies in the case of periodicity along
the direction of the driving current and homogeneity in the transverse
direction, for which the analytical result is confirmed by our numerical
simulation on the constituent spins, and a convenient experiment is proposed.

##  I introduction

Magnetic skyrmion is a kind of noncollinear spin texture with topological
stability, and has attracted a lot of interest ever since it had been
theoretically proposed [1] and experimentally observed [2, 3, 4, 5, 6]. On a
large scale, a magnetic skyrmion behaves as a pointlike object moving on a
two-dimensional space. The presence of both driving force and thermal
fluctuations motivates us to consider using magnetic skyrmion to implement the
thermal ratchet model, which is an important topic with a wide range of
interests [7]. We proposed to use a magnetic skymion adiabatically driven by a
ratchet-like spin-polarized electric currents to implement the adiabatic
thermal ratchet[8], in addition to the stochastic force representing thermal
fluctuations [9], which was the first realization of a thermal ratchet in
terms of skyrmions in a uniform temperature. Other realizations of ratchet
motions of skyrmions include unidirectional rotation driven by thermal
fluctuations in presence of a temperature gradient [10], Magnus-induced
ratchet effects for skyrmions interacting with asymmetric substrates [11, 12],
ratchet motion induced by a biharmonic in-plane magnetic field [13].

For the adiabatially driven skyrmion thermal ratchet [9], the thermal
fluctuations are represented as a stochastic force, and the dynamics is
described by Langevin equation, which was treated by using Fokker-Planck
equation. It was shown that if the driving electric current is periodic along
a specific direction, which is different from that of the current itself, and
is determined by a quantity which is a function of several parameters of the
system, in other words, the periods along and orthogonal to the direction of
the driving current are locked in a specific way, then the time and
statistical average velocity of the skyrmion is proportional to a closed
integral of a curvature of an eigenfunction of an hermitian operator, which is
a similarity transformation of the Fokker-Planck operator. Hence the average
velocity is topologically quantized as a Chern number multiplied by a basic
unit. The result implemented a generalization of a one-dimensional adiabatic
thermal rathet model [8], and provides a novel method of manipulating magnetic
skyrmions. Interesting as it is, this result was under the special condition
concerning the direction of the periodicity, which needs very careful
arrangement in the experiment.

In this paper, we make two new approaches and consider several extensions of
this problem. First, we make an approach based on identifying the drift
direction, along which the component of the velocity depends on the driving
electric current while independent along the direction perpendicular to it.
This approach can apply to the case studied in the previous work [9], which
assumes the periodicity in the drift direction, and is here generalized to a
more general form of locking between periods in the direction of the driving
current and the orthogonal direction, so that there is periodicity along the
drift direction.

The second approach is based on a perturbation theory for the nonhermitian
operator which is a similarity transformation of the Fokker-Planck operator.
This approach applies to the case that the periods along the direction of the
driving electric current and along the orthogonal direction are independent.

These two approaches both apply to a special case, which is also most
practical, that the electric current is periodic along or perpendicular to the
direction of the current itself while homogeneous in the direction orthogonal
to it. This is because the homogeneity can be regarded as the periodicity with
period $0$. Except for this special case, the first approach cannot apply to
the case that the periods in the longitudinal and transverse directions are
independent, as the driving current is then not periodic along the drift
direction.

In each of these two approaches, we find that the time and statistical average
velocity is topologically quantized. For the special case in which the driving
electric current is periodic along its own direction while homogeneous on the
orthogonal direction, we also perform a numerical simulation in terms of the
constituent spins, using the stochastic Landau-Lifschitz-Gilbert equation,
confirming the topological quantization. We also propose a convenient
experimental setup for this special case.

On the theoretical aspect, we find a Chern number in a nonhermitian system.
From the two previous papers [8, 9], to the present paper, this line of
research has been inspired by the analogy with the adiabatic transport of
quantum particle in a periodic potential [14, 15], known as Thouless pump.
Indeed, the theoretical framework based on the perturbation theory has been
inspired by that for the Thouless pump [14, 15]. But the systems considered in
our line of research are classical stochastic systems with thermal
fluctuations, rather than quantum fluctuations, as in Thouless pump. It has
been noted that the skyrmions can manifest quantum behavior at low
temperatures [16]. It is interesting, as the future work, to consider the
coexistence of quantum and thermal fluctuations and combine elements of
Thouless pump and thermal ratchet.

The rest of the paper is organized as the following. In Section II, we
introduce the stochastic motion of the magnetic skyrmion, its description in
terms of Langevin and Fokker-Planck equations, and the similarity
transformation of the Fokker-Planck operator. In Section III, we present the
approach based on identifying the drift direction. In Section IV, we develop a
nonhermitian perturbation theory based on the eigenfunctions. Some details are
given in Appendices. Especially, we discuss the case with independent periods
along and orthogonal to the direction of the driving current, as well as the
special case that the driving electric current is periodic along its own
direction while homogeneous on the orthogonal direction, or vice versa. For
the first special case, we also make the numerical simulation by using the
stochastic Landau-Lifschitz-Gilbert equation, and propose an experiment. A
summary is made in Section V.

##  II Stochastic motion of a magnetic skyrmion

###  II.1 Fokker-Planck Equation and Probability Current

Consider a magnetic skyrmion driven by a spin-polarized electric current in a
two-dimensional space. At a finite temperature, it is subject to both the
driving force and a stochastic force representing the thermal fluctuations.
Its stochastic motion at a finite temperature can be described in terms of a
Langevin equation with a stochastic term [9, 17]

$\alpha_{d}\left[\bm{\dot{q}}-\frac{\beta}{\alpha}\bm{v}_{s}\right]+\alpha_{m}\bm{\hat{z}}\times[\bm{\dot{q}}-\bm{v}_{s}]=\bm{\nu}(t),$
(1)

where the stochastic variable $\bm{q}=(q_{x},q_{y})$ represents the position
of the skyrmion as a whole, $\alpha_{d}\equiv\alpha\iint
dxdy\left(\frac{\partial\bm{n}}{\partial x}\right)^{2}$, where $\alpha$ is the
Gilbert damping coefficient, $\bm{n}$ represents the direction of each
constituent spin, $\alpha_{m}\equiv\iint
dxdy\bm{n}\cdot\left(\frac{\partial\bm{n}}{\partial
x}\times\frac{\partial\bm{n}}{\partial y}\right)$, $\beta$ is the non-
adiabatic coefficient, usually $\beta\ll\alpha$.
$\bm{v}_{s}=-\frac{a^{3}}{2e}\bm{j}$, where $\bm{j}$ is the spin-polarized
electric current density multiplied by its spin polarization and divided by
the magnetic saturation. $\bm{\nu}=(\nu_{x},\nu_{y})$ is the stochastic force
due to the finite temperature, satisfying

$\langle\nu_{i}(t)\rangle=0,\,\langle\nu_{i}(t)\nu_{j}(t^{\prime})\rangle=2\frac{\alpha_{d}k_{B}Ta^{2}}{\hbar}\delta_{ij}\delta(t-t^{\prime}),$
(2)

where $\langle\cdots\rangle$ denotes the statistical ensemble average, $k_{B}$
is the Boltzmann constant; $\hbar$ is the Planck constant, $T$ is the
temperature, $a$ is the lattice constant of the lattice of spins. For the time
being, let us assume that that the electric current is periodic and asymmetric
in $x$ and $y$ directions, with periods $L_{x}$ and $L_{y}$ respectively. It
is also periodic in time with period $\mathcal{T}$.

The instantaneous velocity of the skyrmion is $\dot{\mathbf{q}}$. In view of
its stochastic nature, we will study its statistical ensemble average

$\langle\dot{\mathbf{q}}\rangle,$

which is also called particle current [7].

The statistical nature of the skyrmion can be described in terms of the
probability density $\rho(\bm{r},t)$, that is, $\rho(\bm{r},t)dxdy$ is the
probability that the skyrmion is located in the region $x\sim x+dx,y\sim
y+dy$. $\rho(\bm{r},t)$ can be obtained as the statistical ensemble average of
the constraint that the actual position of the skyrmion $\mathbf{q}(t)$ as a
function of $t$, determined from the Langevin equation, is $\mathbf{r}$, that
is [7],

$\rho(\mathbf{r},t)=\langle\delta(\mathbf{r}-\mathbf{q}(t)\rangle.$ (3)

Then from the Langevin equation, it can be obtained the continuity equation

$\frac{\partial\rho}{\partial t}+\nabla\cdot{\mathcal{J}}=0,$ (4)

which is nothing but the continuity equation, with the probability current
density

${\mathcal{J}}(\mathbf{r},t)=\langle\dot{\mathbf{q}}(t)\delta(\mathbf{r}-\mathbf{q}(t)\rangle.$
(5)

Consequently the total probability current is

$\bm{J}=\iint d^{2}\bm{r}{\mathcal{J}}=\langle\dot{\mathbf{q}}\rangle.$ (6)

That is, the probability current is just the statistical average of the
instantaneous velocity, i.e. the particle current [7].

From the Langevin equation (1), one can derive the Fokker-Planck equation [9,
18, 7]

$-\frac{\partial\rho(\bm{r},t)}{\partial t}={\cal
D}\mathcal{O}\rho(\bm{r},t),$ (7)

where

${\cal
D}\equiv\frac{\alpha_{d}k_{B}Ta^{2}}{\hbar(\alpha_{m}^{2}+\alpha_{d}^{2})}.$
(8)

is the diffusion coefficient,

$\mathcal{O}=-\nabla^{2}+\frac{\partial}{\partial
x}(C_{1}v_{sx}+C_{2}v_{sy})+\frac{\partial}{\partial
y}(-C_{2}v_{sx}+C_{1}v_{sy})$ (9)

is the Fokker-Planck operator, with

$C_{1}\equiv\hbar\frac{\frac{\beta}{\alpha}\alpha_{d}^{2}+\alpha_{m}^{2}}{\alpha_{d}k_{B}Ta^{2}},\,\,C_{2}\equiv\hbar\frac{(\frac{\beta}{\alpha}-1)\alpha_{m}}{k_{B}Ta^{2}}.$
(10)

For simplicity, we define a 2-vector $\bm{G}(x,y,t)$, whose components are

$\begin{array}[]{rl}G_{x}&\equiv C_{1}v_{sx}+C_{2}v_{sy},\\\ G_{y}&\equiv-
C_{2}v_{sx}+C_{1}v_{sy}.\end{array}$ (11)

It is clear that $\bm{G}(x,y,t)$ is periodic in time while it is periodic and
asymmetric in the two space dimensions. Then

$\mathcal{O}=-\nabla^{2}+\nabla\cdot\mathbf{G}=-\nabla\cdot(\nabla-\mathbf{G}).$
(12)

The Fokker-Planck equation can be rewritten as the continuity equation (4),
with a different form of the probability current density where

${\mathcal{J}}={\cal D}(\bm{G}-\nabla)\rho.$ (13)

Therefore the probability current can be obtained as

$\bm{J}=\iint d^{2}\bm{r}{\mathcal{J}}=\iint d^{2}\bm{r}{\cal
D}(\bm{G}-\nabla)\rho,$ (14)

which we will use in the following.

Consider the eigenfunction $\Psi_{n}$ of $\mathcal{O}$, with eigenvalue
$E_{n}$,

$\mathcal{O}\Psi_{n}=E_{n}\Psi_{n}.$ (15)

The real part of each $E_{n}$ is nonnegative, the smallest one being $E_{0}=0$
[18]. The corresponding “ground state” eigenfunction is $\Psi_{0}=\rho_{0}$,
which satisfies

$\mathcal{O}\rho_{0}=0.$ (16)

Hence

$\nabla\cdot(\nabla-\mathbf{G})\rho_{0}=0,$ (17)

which implies

$(\nabla-\mathbf{G})\rho_{0}=\nabla\times{\mathbf{A}},$ (18)

where ${\mathbf{A}}$ is some function, and can be chosen to be
${\mathbf{A}}=A{\mathbf{e}}_{z}$, therefore

$\nabla\rho_{0}=\bm{G}\rho_{0}+\hat{\bm{e}}_{x}\partial_{y}A-\hat{\bm{e}}_{y}\partial_{x}A.$
(19)

Note that the instantaneous eigenfunctions themselves are not solutions to the
time-dependent Fokker-Planck equation, as the Fokker-Planck operator ${\cal
O}$ itself is time-dependent. So $\partial_{t}\Psi_{0}(t)\neq 0$ though
$E_{0}(t)=0$.

###  II.2 Similarity Transformation

We make a similarity transformation, under which each eigenfunction $\Psi_{n}$
is transformed as

$\Psi_{n}\rightarrow\psi_{n}\equiv\rho_{0}^{-\frac{1}{2}}\Psi_{n},$ (20)

with the eigenvalue $E_{n}$ unchanged, satisfying

$\tilde{\mathcal{O}}\psi_{n}=E_{n}\psi_{n},$ (21)

where

$\tilde{\mathcal{O}}\equiv\rho_{0}^{-1/2}\mathcal{O}\rho_{0}^{1/2}$ (22)

is the transformed operator. For the “ground state” $\Psi_{0}=\rho_{0}$, the
transformed eigenfunction is

$\psi_{0}\equiv\rho_{0}^{-\frac{1}{2}}\Psi_{0}=\rho_{0}^{\frac{1}{2}}.$ (23)

Therefore, the similarity transformation can rewritten as

$\begin{split}\psi_{n}&\equiv\rho_{0}^{-\frac{1}{2}}\Psi_{n},\\\
\tilde{\mathcal{O}}&\equiv\psi_{0}^{-1}\mathcal{O}\psi_{0}\\\
&=-\nabla^{2}-\rho_{0}^{-1}(\hat{\bm{e}}_{x}\partial_{y}A-\hat{\bm{e}}_{y}\partial_{x}A)\cdot\nabla+U,\end{split}$
(24)

with

$U\equiv-\Psi_{0}^{-1}\nabla^{2}\Psi_{0}+\nabla\cdot\bm{G}+\Psi_{0}^{-1}\bm{G}\cdot(\nabla\Psi_{0}).$
(25)

Note that when $A$ is independent of $x$ and $y$, $\tilde{\mathcal{O}}$
becomes hermitian, as in our previous work, otherwise, $\tilde{\mathcal{O}}$
is nonhermitian.

Let us define, in general,

$\psi\equiv\rho_{0}^{-\frac{1}{2}}\rho,$ (26)

of which (20) is the case for eigenfunctions. The Fokker-Planck equation (7)
can be rewritten in terms of $\Psi$ and $\tilde{\mathcal{O}}$,

$-\frac{\partial\psi}{\partial t}=\left({\cal
D}\tilde{\mathcal{O}}+\frac{\partial\ln\psi_{0}}{\partial t}\right)\psi.$ (27)

By substituting (26) into (14), the probability current $\bm{J}$ can be
obtained as

$\bm{J}=-2{\cal D}\iint
d^{2}\bm{r}\rho_{0}^{\frac{1}{2}}\left[\nabla+\frac{1}{2}\rho_{0}^{-1}(\hat{\bm{e}}_{x}\partial_{y}A-\hat{\bm{e}}_{y}\partial_{x}A)\right]\psi,$
(28)

##  III  Approach based on identifying the drift direction

Let’s use the orthogonal coordinate system with the direction of the electric
current as the $x$ direction, and the direction orthogonal to it as the $y$
direction. The initial position of the skyrmion is the origin. From the
Langevin equation (1), one can obtain

$\displaystyle\dot{q}_{x}$
$\displaystyle=\frac{\frac{\beta}{\alpha}\alpha_{d}^{2}+\alpha_{m}^{2}}{\alpha_{d}^{2}+\alpha_{m}^{2}}v_{sx}(q_{x},q_{y},t)+\frac{\alpha_{d}}{\alpha_{d}^{2}+\alpha_{m}^{2}}\nu_{x}+\frac{\alpha_{m}}{\alpha_{d}^{2}+\alpha_{m}^{2}}\nu_{y},$
(29) $\displaystyle\dot{q}_{y}$
$\displaystyle=\frac{(-\frac{\beta}{\alpha}+1)\alpha_{d}\alpha_{m}}{\alpha_{d}^{2}+\alpha_{m}^{2}}v_{sx}(q_{x},q_{y},t)+\frac{-\alpha_{m}}{\alpha_{d}^{2}+\alpha_{m}^{2}}\nu_{x}+\frac{\alpha_{d}}{\alpha_{d}^{2}+\alpha_{m}^{2}}\nu_{y}.$
(30)

It is easy to find the direction $v$ in which the velocity component $q_{v}$
is independent of the driving electric current, and the orthogonal direction
$u$,

$\begin{pmatrix}u\\\ v\end{pmatrix}=W\begin{pmatrix}x\\\ y\end{pmatrix}.$ (31)
$\begin{pmatrix}q_{u}\\\ q_{v}\end{pmatrix}=W\begin{pmatrix}q_{x}\\\
q_{y}\end{pmatrix},$ (32)

where

$W=\frac{1}{\sqrt{1+\kappa^{2}}}\begin{pmatrix}1&-\kappa\\\
\kappa&1\end{pmatrix},$ (33)

and

$\kappa=\frac{(\frac{\beta}{\alpha}-1)\alpha_{d}\alpha_{m}}{\frac{\beta}{\alpha}\alpha_{d}^{2}+\alpha_{m}^{2}}.$
(34)

In $uv$ coordinate system, the Langevin equations read

$\displaystyle\dot{q}_{u}$ $\displaystyle={\cal
D}\zeta\sqrt{1+\kappa^{2}}v_{sx}\left(\frac{1}{\sqrt{1+\kappa^{2}}}q_{u}+\frac{\kappa}{\sqrt{1+\kappa^{2}}}q_{v},\frac{-\kappa}{\sqrt{1+\kappa^{2}}}q_{u}+\frac{1}{\sqrt{1+\kappa^{2}}}q_{v},t\right)+\frac{\alpha_{d}}{\alpha_{d}^{2}+\alpha_{m}^{2}}\nu_{u}+\frac{\alpha_{m}}{\alpha_{d}^{2}+\alpha_{m}^{2}}\nu_{v},$
(35) $\displaystyle\dot{q}_{v}$
$\displaystyle=0+\frac{-\alpha_{m}}{\alpha_{d}^{2}+\alpha_{m}^{2}}\nu_{u}+\frac{\alpha_{d}}{\alpha_{d}^{2}+\alpha_{m}^{2}}\nu_{v},$
(36)

where
$\zeta\equiv\hbar\frac{\frac{\beta}{\alpha}\alpha_{d}^{2}+\alpha_{m}^{2}}{\alpha_{d}k_{B}Ta^{2}}$,

$\begin{pmatrix}\nu_{u}\\\ \nu_{v}\end{pmatrix}=W\begin{pmatrix}\nu_{x}\\\
\nu_{y}\end{pmatrix}.$ (37)

It is clear that

$\langle\dot{q}_{v}\rangle=0,$ (38)

as there is no driving term in (36), So we can approximately omit the thermal
drift in $v$ direction, setting

$q_{v}(t)\approx q_{v}(0)=0,$ (39)

and call $u$ direction as the drift direction. Then Eq. (35) becomes

$\dot{q}_{u}\approx{\cal
D}\zeta\sqrt{1+\kappa^{2}}v_{sx}\left(\frac{1}{\sqrt{1+\kappa^{2}}}q_{u},\frac{-\kappa}{\sqrt{1+\kappa^{2}}}q_{u},t\right)+\frac{\alpha_{d}}{\alpha_{d}^{2}+\alpha_{m}^{2}}\nu_{u}+\frac{\alpha_{m}}{\alpha_{d}^{2}+\alpha_{m}^{2}}\nu_{v}.$
(40)

From Eq. (40) and Eq. (36), one can obtain the Fokker-Planck equation [18]

$-\frac{\partial\rho(u,v,t)}{\partial t}={\cal
D}\left\\{-\nabla^{2}\rho(u,v,t)+\zeta\sqrt{1+\kappa^{2}}\left[\frac{\partial}{\partial
u}v_{sx}\left(\frac{1}{\sqrt{1+\kappa^{2}}}u,\frac{-\kappa}{\sqrt{1+\kappa^{2}}}u,t\right)\rho(u,v,t)\right]\right\\}.$
(41)

Then we separate the variables as

$\rho(u,v,t)=\rho_{1}(u,t)\rho_{2}(v).$ (42)

Consequently, we can obtain two equations for $\rho_{1}$ and $\rho_{2}$
respectively

$\displaystyle-\frac{\partial\rho_{1}(u,t)}{\partial t}=-{\cal
D}\frac{\partial^{2}\rho_{1}(u,t)}{\partial u^{2}}+{\cal
D}\zeta\sqrt{1+\kappa^{2}}\frac{\partial}{\partial
u}\left[v_{sx}\left(\frac{1}{\sqrt{1+\kappa^{2}}}u,\frac{-\kappa}{\sqrt{1+\kappa^{2}}}u,t\right)\rho_{1}(u,t)\right]+\lambda\rho_{1}(u,t),$
(43) $\displaystyle\frac{d^{2}\rho_{2}(v)}{dv^{2}}+\frac{\lambda}{{\cal
D}}\rho_{2}(v)=0,$ (44)

where $\lambda$ is an arbitrary constant. The second equation is not important
since we have already obtained the average velocity along the $v$ direction
(38). The first equation can be made simpler by defining
$\rho_{1}^{\prime}\equiv\rho_{1}e^{\lambda t}$,

$-\frac{\partial\rho_{1}^{\prime}(u,t)}{\partial t}={\cal
D}\left[-\frac{\partial^{2}}{\partial u^{2}}+\frac{\partial}{\partial
u}\zeta\sqrt{1+\kappa^{2}}v_{sx}\left(\frac{1}{\sqrt{1+\kappa^{2}}}u,\frac{-\kappa}{\sqrt{1+\kappa^{2}}}u,t\right)\right]\rho_{1}^{\prime}(u,t),$
(45)

which is the same as the Fokker-Planck equation for the one-dimensional
adiabatic particle transport in a periodic ratchet potential [8].

We now suppose the space period of $v_{sx}$ is periodic along $u$ direction,
with period $\mathcal{L}$ and time period $\mathcal{T}$. Then the average
velocity of the skyrmion along the $u$ direction is

$\overline{\langle\dot{q}_{u}\rangle}=\mathcal{C}\frac{\mathcal{L}}{\mathcal{T}}.$
(46)

where $\mathcal{C}$ is the Chern number. From the average velocity along the
$u$ and $v$ direction (46) and (38), we can obtain that along the $x$ and $y$
direction

$\displaystyle\overline{\langle\dot{q}_{x}\rangle}=\frac{1}{\sqrt{1+\kappa^{2}}}\mathcal{C}\frac{\mathcal{L}}{\mathcal{T}},$
(47)
$\displaystyle\overline{\langle\dot{q}_{y}\rangle}=-\frac{\kappa}{\sqrt{1+\kappa^{2}}}\mathcal{C}\frac{\mathcal{L}}{\mathcal{T}}.$
(48)

This recovers the result in our previous work [9], where the hermitian
condition leads to $v_{sx}=v_{sx}(x-\kappa y,t)$, which means that in $uv$
coordinates, $v_{sx}$ only depends on $u\equiv(x-\kappa
y)/\sqrt{1+\kappa^{2}}$, while independent of $v\equiv(\kappa
x+y)/\sqrt{1+\kappa^{2}}$, and it was assumed that the period in $u$ is ${\cal
L}$.

(47) and (48) also apply to a generalized case that
$v_{sx}=v_{sx}(\kappa_{1}x+\kappa_{2}y,t)$ depends on $x$ and $y$ as a
function of $\kappa_{1}x+\kappa_{2}y$. The hermitian case above is its special
case with $\kappa_{1}=1$ and $\kappa_{2}=-\kappa$. As a consequence of (31),
$\kappa_{1}x+\kappa_{2}y=\frac{\kappa_{1}-\kappa_{2}\kappa}{\sqrt{1+\kappa^{2}}}u+\frac{\kappa_{2}+\kappa_{1}\kappa}{\sqrt{1+\kappa^{2}}}v$.
Hence the generalized case can be written as
$v_{sx}=v_{sx}\left(\frac{\kappa_{1}-\kappa_{2}\kappa}{\sqrt{1+\kappa^{2}}}u+\frac{\kappa_{2}+\kappa_{1}\kappa}{\sqrt{1+\kappa^{2}}}v,t\right)$.
Now suppose $v_{sx}$ is also periodic in $x$ and $y$, with periods $L_{x}$ and
$L_{y}$ respectively. In order that $v_{sx}=v_{sx}(\kappa_{1}x+\kappa_{2}y,t)$
is periodic in $\kappa_{1}x+\kappa_{2}y$, it is required that the periods in
$x$ and $y$ are locked as

$\kappa_{1}L_{x}=\kappa_{2}L_{y},$ (49)

which is just the period of $v_{sx}$ in $\kappa_{1}x+\kappa_{2}y$. We are now
considering the generalized case that $v_{sx}(\kappa_{1}x+\kappa_{2}y,t)$ is
periodic in $\kappa_{1}x+\kappa_{2}y$ and is approximately independent on $v$.
Remember
$v_{sx}(\kappa_{1}x+\kappa_{2}y,t)=v_{sx}\left(\frac{\kappa_{1}-\kappa_{2}\kappa}{\sqrt{1+\kappa^{2}}}u+\frac{\kappa_{2}+\kappa_{1}\kappa}{\sqrt{1+\kappa^{2}}}v,t\right)$.
Hence it is periodic in
$\frac{\kappa_{1}-\kappa_{2}\kappa}{\sqrt{1+\kappa^{2}}}u$ with period
$\kappa_{1}L_{x}$. In other words, it is periodic in $u$ with period

$\mathcal{L}=\frac{\kappa_{1}\sqrt{1+\kappa^{2}}}{\kappa_{1}-\kappa_{2}\kappa}L_{x}.$
(50)

Then (47) and (48) can be rewritten as

$\displaystyle\overline{\langle\dot{q}_{x}\rangle}=\frac{\kappa_{1}}{\kappa_{1}-\kappa_{2}\kappa}\mathcal{C}\frac{L_{x}}{\mathcal{T}},$
(51)
$\displaystyle\overline{\langle\dot{q}_{y}\rangle}=-\frac{\kappa_{2}\kappa}{\kappa_{1}-\kappa_{2}\kappa}\mathcal{C}\frac{L_{y}}{\mathcal{T}}.$
(52)

For the case with period-locking but with (49) unsatisfied, $v_{sx}$ is not
periodic in $u$, consequently the present approach does not apply. In general,
the present approach applies to all cases in which the driving current is
periodic along $u$ direction, including the special case that the driving
current is periodic along one of the longitudinal and transverse directions
while homogeneous along the other. This approach does not apply to the case
that the periods are independent and both nonzero along these two directions,
as the driving current now is not periodic along $u$ direction.

##  IV  Approach based on Nonhermitian Adiabatic Perturbation Theory

###  IV.1 Nonhermitian Adiabatic Perturbation Theory

Now we consider another generalization, namely, the case that
$\tilde{\mathcal{O}}$ is nonhermitian.

For this purpose, we develop a nonhermitian adiabatic perturbation theory for
Eq. (27). First we define the instantaneous eigenfunctions of
$\tilde{\mathcal{O}}$. Since $\tilde{\mathcal{O}}$ is not Hermitian, its
eigenfunctions do not necessarily constitute an orthonormal set, that is,
$\iint d^{2}\bm{r}\psi_{m}^{*}\psi_{n}$ is not necessarily equal to
$\delta_{mn}$. Instead, we define the dual of the original eigenfunctions

$\phi_{m}^{*}\equiv\sum_{l}\left(T^{-1}\right)_{ml}\psi_{l}^{*},$ (53)

where $T^{-1}$ is the inverse of $T$, which is defined as
$T_{mn}\equiv\langle\psi_{m}|\psi_{n}\rangle=\iint
d^{2}\bm{r}\psi_{m}^{*}\psi_{n}$. Clearly,

$\langle\phi_{m}|\psi_{n}\rangle\equiv\iint
d^{2}\bm{r}\phi_{m}^{*}\psi_{n}=\delta_{mn}.$

It can be easily confirmed that $\phi_{m}$ is the eigenfunctions of the
operator

$\tilde{\mathcal{O}}^{\dagger}=-\nabla^{2}-\nabla\cdot\rho_{0}^{-1}(\hat{\bm{e}}_{x}\partial_{y}A-\hat{\bm{e}}_{y}\partial_{x}A)+U,$
(54)

with the eigenvalue $E_{m}^{*}$ (see Appendix A). Another important relation
is

$\tilde{\mathcal{O}}^{\dagger}\psi_{0}=0,$ (55)

which indicates that $\psi_{0}$ is also an eigenfunction of
$\tilde{\mathcal{O}}^{\dagger}$ with eigenvalue $0$, so we can define
$\phi_{0}\equiv\psi_{0}$.

The transformed probability density $\psi$ can then be expanded by the
instantaneous eigenfunctions

$\psi=\sum_{n}c_{n}\psi_{n}e^{-{\cal
D}\int_{0}^{t}E_{0}(t^{\prime})dt^{\prime}}.$ (56)

Substitute this into the transformed Fokker-Planck equation (27), calculate
the inner products with $\phi$’s, then we obtain the coefficients through
adiabatic perturbation theory. The final result is

$\psi=\psi_{0}+\sum_{n\neq
0}\frac{2\langle\phi_{n}|\dot{\psi}_{0}\rangle}{{\cal
D}(E_{0}-E_{n})}\psi_{n}.$ (57)

Now we discuss the adiabatic condition. We consider the case that the
potential term dominates the Fokker-Planck operator, that is, the amplitude of
$\bm{G}$,

$G_{0}\gg\frac{1}{a},$ (58)

$U$ in (24) can also be written as

$U=\frac{1}{4}(\nabla\ln\rho_{0})^{2}+\frac{1}{2}(\nabla\cdot\bm{G}).$ (59)

Since generically $\rho_{0}$ is a periodic function, the first term of $U$
possesses a double-well structure. Consequently, the lowest two eigenstates of
the system is degenerate approximately with a small eigenvalue difference due
to the second term of $U$. Thus the band gap of the system can be estimated to
be

$\Delta E\sim\frac{G_{0}}{L},$ (60)

where $L=\max\\{L_{x},L_{y}\\}$. The adiabatic condition is

${\cal T}\gg\frac{1}{{\cal D}\Delta E},$ (61)

where ${\cal T}$ is the time period of the electric current. Hence by
substituting (60) into (61), we obtains

${\cal T}\gg\frac{L}{{\cal D}G_{0}}$ (62)

Substituting (57) into the probability current (28), we obtain

$\bm{J}=-4\sum_{n\neq
0}\frac{\langle\phi_{0}|\left[\nabla+\frac{1}{2}\rho_{0}^{-1}(\hat{\bm{e}}_{x}\partial_{y}A-\hat{\bm{e}}_{y}\partial_{x}A)\right]\psi_{n}\rangle\langle\phi_{n}|\dot{\psi}_{0}\rangle}{E_{0}-E_{n}}.$
(63)

###  IV.2 Topologically Quantized Velocity

$\bm{G}$ is a periodic function, as a linear combination of the two components
of the driving current. For the time being, suppose that periodicities in $x$
and $y$ directions are independent. So $\psi_{n}$ and $\phi_{n}$ must be Bloch
waves satisfying

$\psi_{n\bm{k}}(\bm{r})=e^{i\bm{k}\cdot\bm{r}}w_{n\bm{k}}(\bm{r}),\phi_{n\bm{k}}(\bm{r})=e^{i\bm{k}\cdot\bm{r}}v_{n\bm{k}}(\bm{r}),$
(64)

where $w_{n\bm{k}}(\bm{r})$ and $v_{n\bm{k}}(\bm{r})$ are both periodic
functions. The probability current can be regarded as

$\bm{J}=\bm{J}_{\bm{k}=0},$ (65)

where

$\bm{J}_{\bm{k}}\equiv-2\sum_{n\neq
0}\left(\frac{\langle\phi_{0\bm{k}}|\left[\nabla+\frac{1}{2}\rho_{0}^{-1}(\hat{\bm{e}}_{x}\partial_{y}A-\hat{\bm{e}}_{y}\partial_{x}A)\right]|\psi_{n\bm{k}}\rangle\langle\phi_{n\bm{k}}|\dot{\psi}_{0\bm{k}}\rangle}{E_{0\bm{k}}-E_{n\bm{k}}}+\text{c.c.}\right)\equiv-2(\bm{J}_{\bm{k}}^{h}+\bm{J}_{\bm{k}}^{h*}).$
(66)

We can rewrite $\bm{J}_{\bm{k}}^{h}$ as

$\bm{J}_{\bm{k}}^{h}=\frac{i}{2}\langle\partial_{\bm{k}}v_{0\bm{k}}|\partial_{t}w_{0\bm{k}}\rangle-\frac{i}{2}\langle\partial_{\bm{k}}v_{0\bm{k}}|w_{0\bm{k}}\rangle\langle
v_{0\bm{k}}|\partial_{t}w_{0\bm{k}}\rangle,$ (67)

the derivation of which is given in Appendix C, where it can be seen that $A$
disappears because it is contained in the derivative of an operator with
respect to $\mathbf{k}$.

If the temperature of the system is very low, the potential term dominates the
transformed Fokker-Planck operator $\tilde{\mathcal{O}}$. As a result, the
eigenvalues and the eigenfunctions are insensitive to $\bm{k}$, which means
$\psi_{0\bm{k}}\approx\psi_{0}$ and $\phi_{0\bm{k}}\approx\phi_{0}$. Thus we
obtain $\psi_{0\bm{k}}\approx\phi_{0\bm{k}}$ and $w_{0\bm{k}}\approx
v_{0\bm{k}}$. Consequently $\bm{J}_{\bm{k}}^{h}$ can be approximated by

$\bm{J}_{\bm{k}}^{h}\approx\frac{i}{2}\langle\partial_{\bm{k}}w_{0\bm{k}}|\partial_{t}w_{0\bm{k}}\rangle-\frac{i}{2}\langle\partial_{\bm{k}}w_{0\bm{k}}|w_{0\bm{k}}\rangle\langle
w_{0\bm{k}}|\partial_{t}w_{0\bm{k}}\rangle.$ (68)

From this, we can calculate the total probability current

$\bm{J}_{\bm{k}}=-i(\langle\partial_{\bm{k}}w_{0\bm{k}}|\partial_{t}w_{0\bm{k}}\rangle-\langle\partial_{t}w_{0\bm{k}}|\partial_{\bm{k}}w_{0\bm{k}}\rangle).$
(69)

$\bm{J}_{\bm{k}}$ is insensitive to $\bm{k}$, as demonstrated in Appendix B.
It can also be qualitatively understood in the following way. The dependence
of $J_{\bm{k}}$ on $\bm{k}$ mainly originates from the spatial derive in the
Fokker-Planck operator, which is proportional to temperature, hence is
dominated by other terms at low temperatures.

Then the probability current can be written as

$\displaystyle J_{x}(t)\approx\frac{L_{x}}{2\pi i}\int
dk_{x}(\langle\partial_{k_{x}}w_{0\bm{k}}|\partial_{t}w_{0\bm{k}}\rangle-\langle\partial_{t}w_{0\bm{k}}|\partial_{k_{x}}w_{0\bm{k}}\rangle),$
(70) $\displaystyle J_{y}(t)\approx\frac{L_{y}}{2\pi i}\int
dk_{y}(\langle\partial_{k_{y}}w_{0\bm{k}}|\partial_{t}w_{0\bm{k}}\rangle-\langle\partial_{t}w_{0\bm{k}}|\partial_{k_{y}}w_{0\bm{k}}\rangle).$
(71)

According to (6), the probability current is just the probabilistic average of
the instantaneous velocity of the magnetic skyrmion [7]. Since the driving
electric current is periodic in time, the time average of probabilistic
average of the velocity of the skyrmion is

$\displaystyle\overline{\langle\dot{q}_{x}\rangle}=\frac{L_{x}}{\mathcal{T}}\frac{1}{2\pi
i}\iint
dtdk_{x}(\langle\partial_{k_{x}}w_{0\bm{k}}|\partial_{t}w_{0\bm{k}}\rangle-\langle\partial_{t}w_{0\bm{k}}|\partial_{k_{x}}w_{0\bm{k}}\rangle)=\frac{L_{x}}{\mathcal{T}}\mathcal{C},$
(72)
$\displaystyle\overline{\langle\dot{q}_{y}\rangle}=\frac{L_{y}}{\mathcal{T}}\frac{1}{2\pi
i}\iint
dtdk_{y}(\langle\partial_{k_{y}}w_{0\bm{k}}|\partial_{t}w_{0\bm{k}}\rangle-\langle\partial_{t}w_{0\bm{k}}|\partial_{k_{y}}w_{0\bm{k}}\rangle)=\frac{L_{y}}{\mathcal{T}}\mathcal{C}^{\prime}.$
(73)

where $\mathcal{C}$ and $\mathcal{C}^{\prime}$ are Chern numbers. The above
expressions clearly demonstrate that the average velocity of a magnetic
skyrmion is just a basic unit multiplied by an integer number. This is what we
mean by topological quantization. But notice that our system is a classical
stochastic system.

Notice that he key point is that $\bm{J}_{\bm{k}}$ is insensitive to $\bm{k}$
at low temperature. It doesn’t really matter whether the velocity is averaged
over $k_{x}$, $k_{y}$ or the whole Brillouin zone. The result remains
unchanged.

If we average the velocity over the whole Brillouin zone, the time and
probabilistic average velocity of the $x$-component velocity is

$\begin{split}\overline{\langle\dot{q}_{x}\rangle}=&\frac{L_{y}}{2\pi}\int
dk_{y}\frac{L_{x}}{\mathcal{T}}\frac{1}{2\pi i}\iint
dtdk_{x}(\langle\partial_{k_{x}}w_{0\bm{k}}|\partial_{t}w_{0\bm{k}}\rangle-\langle\partial_{t}w_{0\bm{k}}|\partial_{k_{x}}w_{0\bm{k}}\rangle)\\\
=&\frac{L_{x}}{\mathcal{T}}\frac{L_{y}}{2\pi}\int
dk_{y}\mathcal{C}(k_{y}).\end{split}$ (74)

The insensitivity of $\bm{J}_{\bm{k}}$ to $\bm{k}$ implies the insensitivity
of $\mathcal{C}(k_{y})$ to $k_{y}$, which is enhanced by the feature that the
eigenvalue spectrum is fully gapped at low temperature and that the Chern
number is a topological invariant, which does not change unless the gap is
closed. Thus $\mathcal{C}(k_{y})=\mathcal{C}$ is constant and does not depend
on $k_{y}$. As a result, the average velocity becomes

$\begin{split}\overline{\langle\dot{q}_{x}\rangle}=&\frac{L_{x}}{\mathcal{T}}\left(\frac{L_{y}}{2\pi}\int
dk_{y}\right)\mathcal{C}\\\
=&\frac{L_{x}}{\mathcal{T}}\mathcal{C}.\end{split}$ (75)

One can also start with (63), with the summation over $n$ replaced as a
summation over $n$ and $\bm{k}$, as mentioned by the referee.

Now the probability current can be written as

$\bm{J}=-\frac{4}{N_{x}N_{y}}\sum_{n\neq
0,\bm{k}}\frac{\langle\phi_{0\bm{k}}|\left[\nabla+\frac{1}{2}\rho_{0}^{-1}(\hat{\bm{e}}_{x}\partial_{y}A-\hat{\bm{e}}_{y}\partial_{x}A)\right]\psi_{n\bm{k}}\rangle\langle\phi_{n\bm{k}}|\dot{\psi}_{0\bm{k}}\rangle}{E_{0\bm{k}}-E_{n\bm{k}}},$
(76)

where $N_{x}=2\pi/L_{x}$, $N_{y}=2\pi/L_{y}$, $N_{x}N_{y}$ is the number of
different values of the two-dimensional discrete crystalline momentum.

Then following the method similar to above, one can obtain

$\begin{split}\bm{J}=&-\frac{1}{N_{x}N_{y}}\sum_{\bm{k}}(-i)(\langle\partial_{\bm{k}}w_{0\bm{k}}|\partial_{t}w_{0\bm{k}}\rangle-\langle\partial_{t}w_{0\bm{k}}|\partial_{\bm{k}}w_{0\bm{k}}\rangle)\\\
=&i\frac{L_{x}}{2\pi}\frac{L_{y}}{2\pi}\iint
dk_{x}dk_{y}(\langle\partial_{\bm{k}}w_{0\bm{k}}|\partial_{t}w_{0\bm{k}}\rangle-\langle\partial_{t}w_{0\bm{k}}|\partial_{\bm{k}}w_{0\bm{k}}\rangle).\end{split}$
(77)

Consequently, the time-averaged particle current is

$\begin{split}\overline{\langle\dot{\bm{q}}\rangle}=&\frac{1}{\mathcal{T}}\int
dt\bm{J}\\\ =&\frac{L_{x}}{2\pi}\frac{L_{y}}{2\pi}\frac{1}{\mathcal{T}}\iiint
dk_{x}dk_{y}dt(\langle\partial_{\bm{k}}w_{0\bm{k}}|\partial_{t}w_{0\bm{k}}\rangle-\langle\partial_{t}w_{0\bm{k}}|\partial_{\bm{k}}w_{0\bm{k}}\rangle)\end{split}$
(78)

Then we again arrive at the conclusion that the time average of the particle
current is topologically quantized.

###  IV.3 Discussion

Without loss of generality, suppose that the electric current is along $x$
direction. We now derive a constraint on the relation between the two
components of the average velocity. According to the Langevin equation (1), we
find the following relation

$\displaystyle\overline{\langle\dot{q}_{x}\rangle}=\frac{\frac{\beta}{\alpha}\alpha_{d}^{2}+\alpha_{m}^{2}}{\alpha_{d}^{2}+\alpha_{m}^{2}}\overline{\langle
v_{sx}\rangle},$
$\displaystyle\overline{\langle\dot{q}_{y}\rangle}=\frac{\left(-\frac{\beta}{\alpha}+1\right)\alpha_{d}\alpha_{m}}{\alpha_{d}^{2}+\alpha_{m}^{2}}\overline{\langle
v_{sx}\rangle}.$

Comparing the above two equation, we conclude that
$\overline{\langle\dot{q}_{y}\rangle}$ is proportional to
$\overline{\langle\dot{q}_{x}\rangle}$, as

$\overline{\langle\dot{q}_{y}\rangle}=-\kappa\overline{\langle\dot{q}_{x}\rangle},$
(79)

where

$\kappa\equiv\frac{\left(-\frac{\beta}{\alpha}+1\right)\alpha_{d}\alpha_{m}}{\frac{\beta}{\alpha}\alpha_{d}^{2}+\alpha_{m}^{2}}.$
(80)

This constraint is satisfied by all cases considered in this paper.

In the following, we consider three subcases. In the first subcase, the
electric current is periodic in $x$ direction while constant in $y$ direction,
which is easy to realize in the experiment, as discussed in Section IV.5. As a
result, the $x$ component of the average velocity is quantized, as given in
(72), while the argument for the velocity quantization in the preceding
section does not apply to $y$ component. However, it is obtained from (79)
that
$\overline{\langle\dot{q}_{y}\rangle}=-\kappa\frac{L_{x}}{\mathcal{T}}\mathcal{C},$
which is quantized with a more complicated unit. Hence the result for the
first case is

$\displaystyle\begin{cases}&\overline{\langle\dot{q}_{x}\rangle}=\frac{L_{x}}{\mathcal{T}}\mathcal{C},\\\
&\overline{\langle\dot{q}_{y}\rangle}=-\kappa\frac{L_{x}}{\mathcal{T}}\mathcal{C}.\end{cases}$
(81)

This result can also be obtained in the approach based on the drift direction.
Now $v_{sx}$ is independent of $y$, hence
$v_{sx}={\color[rgb]{1,0,0}v_{sx}}\left(\frac{1}{\sqrt{1+\kappa^{2}}}u,t\right)$.
If the period along the $x$ direction is $L_{x}$, that along the $u$ direction
is $\mathcal{L}=L_{x}\sqrt{1+\kappa^{2}}$. Substituting this relation into Eq.
(47) and Eq. (48), we can reproduce (83).

In the second case, the electric current is periodic in the $y$ direction
while constant in the $x$ direction. Consequently the average velocity along
the $y$ direction satisfies Eq. (73), while it is the average velocity along
the $x$ direction that is obtained from Eq. (79), as
$\overline{\langle\dot{q}_{x}\rangle}=-\frac{1}{\kappa}\frac{L_{y}}{\mathcal{T}}\mathcal{C}^{\prime},$
which is quantized with a more complicated unit. Hence the result for the
second case is

$\displaystyle\begin{cases}&\overline{\langle\dot{q}_{x}\rangle}=-\frac{1}{\kappa}\frac{L_{y}}{\mathcal{T}}\mathcal{C}^{\prime},\\\
&\overline{\langle\dot{q}_{y}\rangle}=\frac{L_{y}}{\mathcal{T}}\mathcal{C}^{\prime}.\end{cases}$
(82)

This can also be reproduced in the approach based on the drift direction, in a
way similar to the first case.

In the third case, the electric current is periodic in both $x$ and $y$
direction, and the periods are unrelated. This situation is difficult to
realize in the experiment. Since the relation between the average velocities
along the two directions satisfy (79). There are two possibilities,

$\displaystyle\begin{cases}&\overline{\langle\dot{q}_{x}\rangle}=\frac{L_{x}}{\mathcal{T}}\mathcal{C},\\\
&\overline{\langle\dot{q}_{y}\rangle}=-\kappa\frac{L_{x}}{\mathcal{T}}\mathcal{C},\end{cases}$
(83)

or

$\displaystyle\begin{cases}&\overline{\langle\dot{q}_{x}\rangle}=-\frac{1}{\kappa}\frac{L_{y}}{\mathcal{T}}\mathcal{C}^{\prime},\\\
&\overline{\langle\dot{q}_{y}\rangle}=\frac{L_{y}}{\mathcal{T}}\mathcal{C}^{\prime}.\end{cases}$
(84)

They cannot be reproduced in the the approach based on the drift direction.

What those Chern numbers are exactly, and which of the two possibilites
actually appears in the third case, are determined by the driving electric
current.

###  IV.4 Numerical Simulation

In order to confirm our theoretical result, we perform a numerical simulation
of the stochastic Landau-Lifschitz-Gilbert equation [19, 20, 21, 10, 17, 22,
23, 9]

$\frac{\partial\bm{n}}{\partial
t}+(\bm{v}_{s}\cdot\nabla)\bm{n}=-\frac{1}{\hbar}\bm{n}\times(\bm{H}_{eff}+\bm{R})+\alpha\bm{n}\times\frac{\partial\bm{n}}{\partial
t}+\beta\bm{n}\times(\bm{v}_{s}\cdot\nabla)\bm{n},$ (85)

which describes the dynamics of the constituent spins of the magnetic
skyrmion. $\bm{H}_{eff}\equiv-\frac{\partial\bm{H}_{S}}{\partial\bm{n}}$ is
the effective magnetic field, where the skyrmion Hamiltonian is

$\bm{H}_{S}=-J\sum_{\langle ij\rangle}\bm{n}_{i}\cdot\bm{n}_{j}-D\sum_{\langle
ij\rangle}\bm{\hat{e}}_{ij}\cdot{\bm{n}_{i}\times\bm{n}_{j}}-\bm{B}\cdot\sum_{i}\bm{n}_{i}-K\sum_{i}n_{iz}^{2}.$
(86)

In this equation, $J$ is the exchange interaction constant, $D$ is the
Dzyaloshinskii-Moriya interaction constant [24, 25], $\bm{B}$ is the magnetic
field, $K$ is the anisotropic constant, $\bm{R}$ is the random magnetic field,
which characterizes the effect of the finite temperature $T$, with $\langle
R_{i}(\bm{r},t)\rangle=0$, $\langle
R_{i}(\bm{r},t)R_{j}(\bm{r}^{\prime},t^{\prime})\rangle=2\alpha\hbar
k_{B}Ta^{2}\delta_{ij}\delta(\bm{r}-\bm{r}^{\prime})\delta(t-t^{\prime}),$
where $i,j=x,y,z$.

The simulation is performed on a $100\times 100$ lattice, which means
$L_{x}=L_{y}=100$. The Gilbert damping constant is $\alpha=0.1$. The non-
adiabatic spin transfer torque constant is $\beta=0$. The Dzyaloshinski-Moriya
interaction constant is $D=0.12J$. The magnetic field is $B=0.015J$. The
anisotropic energy constant is $K=0.01J$. The electric current density is
assumed to be [26, 9]

$\bm{j}=\frac{2e}{a^{2}\tau}\left[-j_{c}\left(\cos\frac{2\pi}{L_{x}}x+\frac{1}{2}\cos\frac{4\pi}{L_{x}}x\right)-A\cos\frac{2\pi}{\mathcal{T}}t\right]\hat{\bm{e}}_{x},$
(87)

where $A=0.2$ and $j_{c}=0.08\sim 0.2$, $\tau\equiv\frac{\hbar}{J}$ is the
time unit. $\bm{j}$ is periodic in $x$ direction while homogeneous in $y$
direction. We use the Runge-Kutta method of fourth order, while the time step
is chosen to be $0.1\tau$. The choice of the time period $\mathcal{T}$ must
satisfy the adiabatic condition (62), under which our adiabatic perturbation
theory applies. We have done the simulation for several values of temperature
relative to $J/k_{B}$, given as $k_{B}T/J=0.001,0.01,0.1$.

According to the definition of $\bm{G}$, the amplitude $G_{0}$ of $\bm{G}$,can
be approximated as

$G_{0}\sim\frac{\frac{\beta}{\alpha}\alpha_{d}^{2}+\alpha_{m}^{2}}{\alpha_{d}k_{B}Ta^{2}}\hbar\frac{a}{\tau}(j_{c}+A).$
(88)

In our simulation, the corresponding parameters are $\alpha_{m}=-12.2296\sim
10,\alpha_{d}=1.41767\sim 1,j_{c}+A\sim 0.1$, so $G_{0}$ is approximated by

$G_{0}\sim\frac{1}{a}\left(\frac{k_{B}T}{J}\right)^{-1}\times 10.$ (89)

In the deterministic limit, (58) must be satisfied, which means
$\frac{k_{B}T}{J}\ll 10$, namely in the low temperature regime. This is
actually the case discussed above in the theoretical sections. Then by
substituting the expressions in (8) and (89) for the certain terms in (62) and
making some approximations, we can get the explicit adiabatic condition

$\mathcal{T}\gg 10^{3}\tau.$ (90)

Therefore, $\mathcal{T}=5000\tau$ is chosen for the simulation.

From the above parameters, the theoretical values of the two components of the
average velocity can be obtained from (83) as

$\displaystyle\overline{\langle\dot{q}_{x}\rangle}=\mathcal{C}\frac{L_{x}}{\mathcal{T}}=\mathcal{C}\times
0.02\frac{a}{\tau},$ (91)
$\displaystyle\overline{\langle\dot{q}_{y}\rangle}=\frac{1.41767}{-12.2296}\mathcal{C}\times
0.02\frac{a}{\tau}=\mathcal{C}\times\left(-0.00231759\frac{a}{\tau}\right).$
(92)

In the simulation, we obtain the average velocity, which is averaged over ten
periods, versus the parameter $j_{c}$ for different temperatures, represented
as multiplies of exchange interaction constant $J$. The results are shown in
FIG. 1.

It is clear that the average velocity of the skyrmion at a low temperature is
indeed quantized as given theoretically in (91) and (92).

Figure 1: The skyrmion’s average velocities along $x$ and $y$ directions as
functions of the amplitude $j_{c}$ of the polarized electric current. The unit
of the velocity is $\frac{a}{\tau}$. Different symbols and colours represent
different values of $k_{B}T$ in unit of $J$. Black curves and squares
represent results for $k_{B}T=0.001J$; red curves and circles represent
results for $k_{B}T=0.01J$; green curves and triangles represent results for
$k_{B}T=0.1J$. The grey line represents the analytically predicted value of
the velocity.

###  IV.5 Experimental Proposal

In the above simulation, the electric current density possesses the form (87),
which is not easy to realize in the experiment since it is difficult to make
the electric current vary with position as trigonometric functions. However,
by using the method we have used in our previous work [9], we can replace the
trigonometric function with the following function

$\displaystyle f(x)=\begin{cases}1.5,&0\leqslant x<20a,\\\ -1&20a\leqslant
x<80a,\\\ 1.5,&80a\leqslant x<100a.\end{cases}$ (93)

Furthermore, $f(x+L_{x})=f(x)$. As a result, the electric current density can
be written as

$\bm{j}=\frac{2e}{a^{2}\tau}\left[-j_{c}f(x)-A\cos\frac{2\pi}{\mathcal{T}}t\right]\hat{\bm{e}}_{x}.$
(94)

In order to realize the above ratchetlike electric current, we devise the
experiment as shown in Fig. 2. The thick lines are all the electrodes with
different electric voltages. The distance between the blue and the electrodes
is $l_{1}$ while it between the red one and the green one is $l_{2}$. On the
other hand, the distance between the neighboring green and blue electrodes
must be as small as possible so that the electric current between them does
not affect the motion of the magnetic skyrmion much. In our simulation,
$l_{1}=40a$ and $l_{2}=60a$. The actual values are not essential.

Figure 2: The experimental realization of the ratchetlike electric current.
The thick lines with different colors represent the electrodes with different
electric voltages.

The red electrodes are all grounded, which means

$V_{+}=0.$ (95)

The voltage of each blue electrode is

$V_{-}=-\left(1.5j_{c}+A\cos\frac{2\pi t}{\cal
T}\right)\frac{2e}{a^{2}\tau}\frac{l_{1}}{\sigma},$ (96)

where $\sigma$ is the electrical conductivity of the material. The voltage of
each green electrode is

$V_{-}^{\prime}=-\left(1.0j_{c}-A\cos\frac{2\pi t}{\cal
T}\right)\frac{2e}{a^{2}\tau}\frac{l_{2}}{\sigma}.$ (97)

Then the electric current density in different region of the sample is as
described by Eq. (94).

In the actual experiment, we can first generate a single magnetic skyrmion on
the sample where the electrodes are mounted in advance. Then we apply the
above electric voltages to the electrodes and the magnetic skyrmion start
moving. One measures the change of the position of the skyrmion as a function
of time, from which the instantaneous velocity of the skyrmion can be
calculated. Finally, the average velocity of the skyrmion can be obtained by
averaging over several periods.

##  V Summary

We have studied in details the two-dimensional stochastic motion of a magnetic
skyrmion driven by a generic spin-polarized electric current which is periodic
in time while periodic and asymmetric in the direction of the electric current
or in the transverse direction, or in both directions. In any case, the
average velocities along the two directions are shown to be proportional, with
the proportional factor given by the drift direction.

We have considered some general cases significantly beyond the special case
considered in our previous work, in which the periods in the longitudinal and
transverse directions are locked in a special way such that the superposed
periodicity is along the drift direction, which is determined by the
parameters of the system.

We have made an approach based on identifying the drift direction, which
applies to a more general case of period-locking, of which the case treated in
our previous work is a special one. If the adiabatic condition is satisfied,
the time and probabilistic average of the velocity component along the drift
direction is the basic unit, which is the ratio between the space period along
this direction and the time period, multiplied by a Chern number. The average
velocity along the longitudinal and transverse directions can be obtained as
components. Consequently, the average velocity along any direction, as a
projection of that along the drift direction, is quantized.

We have also made a second approach and developed a formalism based on the
eigenfunctions of the nonhermitian similarity transformation of the Fokker-
Planck operator, and it is assumed that the periods along the longitudinal and
transverse direction are independent.

In case the driving current is periodic along one of these two direction while
homogeneous along the other, the average velocity along this direction is the
basic unit multiplied by a Chern number. Multiplying it by the proportional
factor mentioned above gives the average velocity along the orthogonal
direction. This result can be obtained using either of the two approaches. For
the first approach to be applicable, the periods along the longitudinal and
transverse directions should be in a way that lead to periodicity along the
drift direction. This requirement may not be satisfied if the periods along
those two directions are independent and both nonzero.

For the case that the driving current is periodic along its own direction
while homogeneous in the transverse direction, we have also performed a
numerical simulation which confirms our theoretical prediction, and have
proposed the experimental setup to realize this case, which is more convenient
than that in our previous work [9], in which the electric current must be in
the form of $f(x-\kappa y)$, where $x$ and $y$ are the spatial coordinates,
$\kappa$ is the proportional factor.

The topological quantization provides a method to robustly manipulate the
magnetic skyrmions at a low temperature, which may be useful in memory storage
and communication.

###### Acknowledgements.

We thank Qian Niu for useful discussions. This work was supported by National
Science Foundation of China (Grant No. 12075059).

##  Appendix A Eigenfunctions of nonhermitian Operators

For a nonhermitian operator $O$, define a set of orthogonal basis functions
$f_{n}$, $n=1,2,...$, with

$\langle f_{m}|f_{n}\rangle\equiv\int d\tau f_{m}^{*}f_{n}=\delta_{mn}.$ (98)

Then a matrix $O$ can be defined with the matrix elements

$O_{mn}\equiv\langle f_{m}|O|f_{n}\rangle\equiv\int d\tau f_{m}^{*}Of_{n}.$
(99)

Suppose $\det{O}\neq 0$, then the matrix can be diagonalized through the
similarity transformation

$P^{-1}OP=diag(E_{1},E_{2},...)\equiv E,$ (100)

where $\\{E_{n}\\}$ are eigenvalues. Thus

$OP=PE.$ (101)

Therefore, the eigenvectors of $O$ are

$a_{1}=\begin{pmatrix}P_{11}\\\ P_{21}\\\
\vdots\end{pmatrix},a_{2}=\begin{pmatrix}P_{12}\\\ P_{22}\\\
\vdots\end{pmatrix},\cdots,a_{n}=\begin{pmatrix}P_{1n}\\\ P_{2n}\\\
\vdots\end{pmatrix},\cdots,$ (102)

with $E_{1}$, $E_{2}$, $\cdots$, $E_{n}$, $\cdots$. That is,

$Oa_{n}=E_{n}a_{n}.$ (103)

On the other hand,

$P^{-1}O=EP^{-1},$ (104)

or

$O^{\dagger}(P^{-1})^{\dagger}=(P^{-1})^{\dagger}E^{\dagger},$ (105)

which implies that the eigenvectors of $O^{\dagger}$ are

$b_{1}=\begin{pmatrix}(P^{-1})^{\dagger}_{11}\\\ (P^{-1})^{\dagger}_{21}\\\
\vdots\end{pmatrix},b_{2}=\begin{pmatrix}(P^{-1})^{\dagger}_{12}\\\
(P^{-1})^{\dagger}_{22}\\\
\vdots\end{pmatrix},\cdots,b_{n}=\begin{pmatrix}(P^{-1})^{\dagger}_{1n}\\\
(P^{-1})^{\dagger}_{2n}\\\ \vdots\end{pmatrix},\cdots,$ (106)

with eigenvalues $E_{1}^{*}$, $E_{2}^{*}$, $\cdots$, $E_{n}^{*}$, $\cdots$.
That is,

$O^{\dagger}b_{n}=E_{n}^{\dagger}b_{n}.$ (107)

It is straightforward to confirm

$b_{m}^{\dagger}a_{n}=((P^{-1})^{\text{T}}_{1m},(P^{-1})^{\text{T}}_{2m},\cdots)\begin{pmatrix}P_{1n}\\\
P_{2n}\\\
\vdots\end{pmatrix}=(P^{-1})_{m1}P_{1n}+(P^{-1})_{m2}P_{2n}+\cdots=(P^{-1}P)_{mn}=\delta_{mn}.$
(108)

The eigenfunctions of the operator $O$ and $O^{\dagger}$ can be obtained as

$\displaystyle\psi_{n}=$ $\displaystyle\sum_{i}P_{in}f_{i},$ (109)
$\displaystyle\phi_{n}=$ $\displaystyle\sum_{i}(P^{-1})^{\dagger}_{in}f_{i}.$
(110)

By using Dirac notation, the operator $\hat{A}$ can be written as

$O=\sum_{mn}O_{mn}|f_{m}\rangle\langle f_{n}|.$

Therefore

$\displaystyle O|\psi_{n}\rangle=$
$\displaystyle\sum_{ml}O_{ml}|f_{m}\rangle\langle
f_{l}|\sum_{i}P_{in}|f_{i}\rangle$ $\displaystyle=$
$\displaystyle\sum_{iml}O_{ml}P_{in}\delta_{il}|f_{m}\rangle=\sum_{ml}O_{ml}P_{ln}|f_{m}\rangle=\sum_{ml}P_{ml}E_{ln}|f_{m}\rangle=E_{n}\sum_{m}P_{mn}|f_{m}\rangle=E_{n}|\psi_{n}\rangle,$
$\displaystyle O^{\dagger}|\phi_{n}\rangle=$
$\displaystyle\sum_{ml}O^{\dagger}_{ml}|f_{m}\rangle\langle
f_{l}|\sum_{i}(P^{-1})^{\dagger}_{in}|f_{i}\rangle$ $\displaystyle=$
$\displaystyle\sum_{iml}O^{\dagger}_{ml}(P^{-1})^{\dagger}_{in}\delta_{il}|f_{m}\rangle=\sum_{ml}O^{\dagger}_{ml}(P^{-1})^{\dagger}_{ln}|f_{m}\rangle=\sum_{ml}(P^{-1})^{\dagger}_{ml}E^{\dagger}_{ln}|f_{m}\rangle=E^{*}_{n}\sum_{m}(P^{-1})_{mn}|f_{m}\rangle$
$\displaystyle=$ $\displaystyle E^{*}_{n}|\phi_{n}\rangle,$

which confirms that $|\psi_{n}\rangle$ and $|\phi_{n}\rangle$ are indeed
eigenfunctions of $O$ and $O^{\dagger}$, respectively. Now the inner products
of these eigenfunctions can be calculated as

$\displaystyle\langle\psi_{m}|\psi_{n}\rangle=$
$\displaystyle\sum_{ij}P^{*}_{im}P_{jn}\langle
f_{i}|f_{j}\rangle=\sum_{ij}(P^{\dagger})_{mi}P_{jn}\delta_{ij}=(P^{\dagger}P)_{mn}=T_{mn},$
(111) $\displaystyle\langle\phi_{m}|\psi_{n}\rangle=$
$\displaystyle\sum_{ij}\left[(P^{-1})^{\dagger}_{im}\right]^{*}P_{jn}\langle
f_{i}|f_{j}\rangle=\sum_{ij}(P^{-1})_{mi}P_{jn}\delta_{ij}=(P^{-1}P)_{mn}=\delta_{mn}.$
(112)

These inner products help figure out whether the orthogonal partner defined in
(53) is the eigenfunctions of the hermitian conjugate operator obtained here.

$\displaystyle\sum_{l}(T^{-1})_{ml}\psi_{l}^{*}=$
$\displaystyle\sum_{ln}(P^{-1})_{mn}(P^{-1})^{\dagger}_{nl}\sum_{i}P^{*}_{il}f_{i}^{*}$
$\displaystyle=$
$\displaystyle\sum_{iln}(P^{-1})_{mn}(P^{-1})^{\dagger}_{nl}P^{\dagger}_{li}f_{i}^{*}=\sum_{in}(P^{-1})_{mn}\delta_{ni}f_{i}^{*}=\sum_{n}(P^{-1})_{mn}f_{n}^{*}$
$\displaystyle=$
$\displaystyle\left(\sum_{n}(P^{-1})^{\dagger}_{nm}f_{n}\right)^{*}=\phi_{m}^{*}.$

##  Appendix B Insensitivity of $\bm{J}_{\bm{k}}$ to $\bm{k}$

The following Hermitian and antihermitian operators can be obtained from the
transformed Fokker-Planck operator (24)

$\displaystyle L_{H}=$
$\displaystyle\frac{\tilde{\mathcal{O}}+\tilde{\mathcal{O}}^{\dagger}}{2}=-\nabla^{2}+\frac{1}{2}(\nabla\cdot\bm{M})+U,$
(113) $\displaystyle L_{A}=$
$\displaystyle\frac{\tilde{\mathcal{O}}-\tilde{\mathcal{O}}^{\dagger}}{2}=-\bm{M}\cdot\nabla-\frac{1}{2}(\nabla\cdot\bm{M}),$
(114)

where
$\bm{M}\equiv\rho_{0}^{-1}(\hat{\bm{e}}_{x}\partial_{y}A-\hat{\bm{e}}_{y}\partial_{x}A)$.
Then we can construct an operator from the above two operators [18]

$\mathcal{H}=L_{H}-i\eta L_{A}=-\nabla^{2}+i\eta\bm{M}\cdot\nabla+U^{\prime},$
(115)

where $U^{\prime}\equiv U+\frac{1+i\eta}{2}(\nabla\cdot\bm{M})$.

It can be seen that when $\eta=i$, $\mathcal{H}=\tilde{\mathcal{O}}$. When
$\eta$ is real, $\mathcal{H}$ is Hermitian.

For the time being, we assume $\eta$ is real. The eigenfunctions of
$\mathcal{H}$ are $\psi^{\prime}_{n\bm{k}}$, satisfying

$\mathcal{H}(\eta)\psi^{\prime}_{n\bm{k}}(\eta)=E^{\prime}_{n\bm{k}}(\eta)\psi^{\prime}_{n\bm{k}}(\eta).$

They are of course Bloch wave functions and their periodic parts are
$w^{\prime}_{n\bm{k}}$’s, which satisfy

$\mathcal{H}^{\prime}(\eta)w^{\prime}_{n\bm{k}}(\eta)=E^{\prime}_{n\bm{k}}(\eta)w^{\prime}_{n\bm{k}}(\eta),$

where

$\mathcal{H}^{\prime}\equiv-(\nabla+i\bm{k})^{2}+i\eta\bm{M}\cdot(\nabla+i\bm{k})+U^{\prime}.$
(116)

Then we obtain the probability current

$\displaystyle\bm{J}^{\prime}_{\bm{k}}(\eta)=$
$\displaystyle-i(\langle\partial_{\bm{k}}w^{\prime}_{0\bm{k}}|\partial_{t}w^{\prime}_{0\bm{k}}\rangle-\langle\partial_{t}w^{\prime}_{0\bm{k}}|\partial_{\bm{k}}w^{\prime}_{0\bm{k}}\rangle)=-2(\bm{J}_{\bm{k}}^{\prime
h}+\bm{J}_{\bm{k}}^{\prime h*}),$ (117)

where

$\displaystyle\bm{J}_{\bm{k}}^{\prime h}(\eta)=$
$\displaystyle\frac{i}{2}\langle\partial_{\bm{k}}w^{\prime}_{0\bm{k}}|\partial_{t}w^{\prime}_{0\bm{k}}\rangle-\frac{i}{2}\langle\partial_{\bm{k}}w^{\prime}_{0\bm{k}}|w^{\prime}_{0\bm{k}}\rangle\langle
w^{\prime}_{0\bm{k}}|\partial_{t}w^{\prime}_{0\bm{k}}\rangle.$ (118)

Following a method of Niu and Thouless for the quantized adiabatic charge
transport [15], we first write $\bm{J}^{\prime}_{\bm{k}}$ in the form of Green
functions and then prove its insensitivity to $\bm{k}$ in the following.

From (118),

$\begin{split}\bm{J}_{\bm{k}}^{\prime
h}(\eta)=&\frac{i}{2}\langle\partial_{\bm{k}}w^{\prime}_{0\bm{k}}|\partial_{t}w^{\prime}_{0\bm{k}}\rangle-\frac{i}{2}\langle\partial_{\bm{k}}w^{\prime}_{0\bm{k}}|w^{\prime}_{0\bm{k}}\rangle\langle
w^{\prime}_{0\bm{k}}|\partial_{t}w^{\prime}_{0\bm{k}}\rangle\\\
=&\frac{i}{2}\sum_{n\neq
0}\langle\partial_{\bm{k}}w^{\prime}_{0\bm{k}}|w^{\prime}_{n\bm{k}}\rangle\langle
w^{\prime}_{n\bm{k}}|\partial_{t}w^{\prime}_{0\bm{k}}\rangle\\\
=&-\frac{i}{2}\sum_{n\neq 0}\langle
w^{\prime}_{0\bm{k}}|\partial_{\bm{k}}w^{\prime}_{n\bm{k}}\rangle\langle
w^{\prime}_{n\bm{k}}|\partial_{t}w^{\prime}_{0\bm{k}}\rangle.\end{split}$
(119)

Then calculate the derivatives of both sides of
$\mathcal{H}^{\prime}(\eta)w^{\prime}_{n\bm{k}}(\eta)=E^{\prime}_{n\bm{k}}(\eta)w^{\prime}_{n\bm{k}}(\eta)$
with respect to $\bm{k}$,

$\frac{\partial\mathcal{H}^{\prime}}{\partial\bm{k}}w^{\prime}_{n\bm{k}}+\mathcal{H}^{\prime}\frac{\partial
w^{\prime}_{n\bm{k}}}{\partial\bm{k}}=\frac{\partial
E^{\prime}_{n\bm{k}}}{\partial\bm{k}}w^{\prime}_{n\bm{k}}+E^{\prime}_{n\bm{k}}\frac{\partial
w^{\prime}_{n\bm{k}}}{\partial\bm{k}}.$ (120)

The inner product of both sides of the above equation with
$w^{\prime}_{0\bm{k}}$, for $n\neq 0$, leads to

$\begin{split}&\langle
w^{\prime}_{0\bm{k}}|\frac{\partial\mathcal{H}^{\prime}}{\partial\bm{k}}|w^{\prime}_{n\bm{k}}\rangle+E^{\prime}_{0\bm{k}}\langle
w^{\prime}_{0\bm{k}}|\frac{\partial
w^{\prime}_{n\bm{k}}}{\partial\bm{k}}\rangle=E^{\prime}_{n\bm{k}}\langle
w^{\prime}_{0\bm{k}}|\frac{\partial
w^{\prime}_{n\bm{k}}}{\partial\bm{k}}\rangle,\end{split}$ (121)

therefore,

$\begin{split}\langle w^{\prime}_{0\bm{k}}|\frac{\partial
w^{\prime}_{n\bm{k}}}{\partial\bm{k}}\rangle=-\frac{1}{E^{\prime}_{0\bm{k}}-E^{\prime}_{n\bm{k}}}\langle
w^{\prime}_{0\bm{k}}|\frac{\partial\mathcal{H}^{\prime}}{\partial\bm{k}}|w^{\prime}_{n\bm{k}}\rangle.\end{split}$
(122)

which is substituted into Eq. (119) to obtain

$\begin{split}\bm{J}_{\bm{k}}^{\prime h}(\eta)=&\frac{i}{2}\sum_{n\neq
0}\frac{\langle
w^{\prime}_{0\bm{k}}|\frac{\partial\mathcal{H}^{\prime}}{\partial\bm{k}}|w^{\prime}_{n\bm{k}}\rangle\langle
w^{\prime}_{n\bm{k}}|\partial_{t}w^{\prime}_{0\bm{k}}\rangle}{E^{\prime}_{0\bm{k}}-E^{\prime}_{n\bm{k}}}\\\
=&\frac{i}{2}\sum_{n\neq 0}\frac{\langle
w^{\prime}_{0\bm{k}}|(-2i)(\nabla+i\bm{k}-\frac{i\eta}{2}\bm{M})|w^{\prime}_{n\bm{k}}\rangle\langle
w^{\prime}_{n\bm{k}}|\partial_{t}w^{\prime}_{0\bm{k}}\rangle}{E^{\prime}_{0\bm{k}}-E^{\prime}_{n\bm{k}}}\\\
=&\sum_{n\neq
0}\frac{\langle\psi^{\prime}_{0\bm{k}}|(\nabla-\frac{i\eta}{2}\bm{M})|\psi^{\prime}_{n\bm{k}}\rangle\langle\psi^{\prime}_{n\bm{k}}|\partial_{t}\psi^{\prime}_{0\bm{k}}\rangle}{E^{\prime}_{0\bm{k}}-E^{\prime}_{n\bm{k}}},\end{split}$
(123)

where we have used

$\frac{\partial\mathcal{H}^{\prime}}{\partial\bm{k}}=-2i(\nabla+i\bm{k}-\frac{i\eta}{2}\bm{M}),$

and

$\psi_{n\bm{k}}=e^{i\bm{k}\bm{x}}w_{n\bm{k}}.$

Since
$\mathcal{H}\psi^{\prime}_{n\bm{k}}=E^{\prime}_{n\bm{k}}\psi^{\prime}_{n\bm{k}}$,
one obtains

$\begin{split}&\left(\nabla-\frac{i\eta}{2}\bm{M}\right)(\mathcal{H}\psi^{\prime}_{n\bm{k}})=\left(\nabla-\frac{i\eta}{2}\bm{M}\right)(E^{\prime}_{n\bm{k}}\psi^{\prime}_{n\bm{k}}),\end{split}$
(124)

Thus

$\begin{split}&\left[\nabla-\frac{i\eta}{2}\bm{M},\mathcal{H}\right]\psi^{\prime}_{n\bm{k}}+\mathcal{H}\left(\nabla-\frac{i\eta}{2}\bm{M}\right)\psi^{\prime}_{n\bm{k}}=E^{\prime}_{n\bm{k}}\left(\nabla-\frac{i\eta}{2}\bm{M}\right)\psi^{\prime}_{n\bm{k}},\end{split}$
(125)

Therefore

$\begin{split}&\langle\psi^{\prime}_{0\bm{k}}|\left[\nabla-\frac{i\eta}{2}\bm{M},\mathcal{H}\right]|\psi^{\prime}_{n\bm{k}}\rangle+\langle\psi^{\prime}_{0\bm{k}}|\mathcal{H}\left(\nabla-\frac{i\eta}{2}\bm{M}\right)|\psi^{\prime}_{n\bm{k}}\rangle=E^{\prime}_{n\bm{k}}\langle\psi^{\prime}_{0\bm{k}}|\left(\nabla-\frac{i\eta}{2}\bm{M}\right)|\psi^{\prime}_{n\bm{k}}\rangle,\end{split}$
(126)

Hence

$\begin{split}&\langle\psi^{\prime}_{0\bm{k}}|\left(\nabla-\frac{i\eta}{2}\bm{M}\right)|\psi^{\prime}_{n\bm{k}}\rangle=\frac{\langle\psi^{\prime}_{0\bm{k}}|\left[\mathcal{H},\nabla-\frac{i\eta}{2}\bm{M}\right]|\psi^{\prime}_{n\bm{k}}\rangle}{E^{\prime}_{0\bm{k}}-E^{\prime}_{n\bm{k}}}.\end{split}$
(127)

On the other hand,

$\begin{split}&\dot{\mathcal{H}}\psi^{\prime}_{0\bm{k}}+\mathcal{H}\partial_{t}\psi^{\prime}_{0\bm{k}}=\dot{E}^{\prime}_{0\bm{k}}\psi^{\prime}_{0\bm{k}}+E^{\prime}_{0\bm{k}}\partial_{t}\psi^{\prime}_{0\bm{k}}.\end{split}$
(128)

thus

$\begin{split}&\langle\psi^{\prime}_{n\bm{k}}|\dot{\mathcal{H}}|\psi^{\prime}_{0\bm{k}}\rangle+\langle\psi^{\prime}_{n\bm{k}}|\mathcal{H}|\partial_{t}\psi^{\prime}_{0\bm{k}}\rangle=\dot{E}^{\prime}_{0\bm{k}}\langle\psi^{\prime}_{n\bm{k}}|\psi^{\prime}_{0\bm{k}}\rangle+E^{\prime}_{0\bm{k}}\langle\psi^{\prime}_{n\bm{k}}|\partial_{t}\psi^{\prime}_{0\bm{k}}\rangle,\end{split}$
(129)

therefore

$\begin{split}&\langle\psi^{\prime}_{n\bm{k}}|\dot{\mathcal{H}}|\psi^{\prime}_{0\bm{k}}\rangle+E^{\prime}_{n\bm{k}}\langle\psi^{\prime}_{n\bm{k}}|\partial_{t}\psi^{\prime}_{0\bm{k}}\rangle=E^{\prime}_{0\bm{k}}\langle\psi^{\prime}_{n\bm{k}}|\partial_{t}\psi^{\prime}_{0\bm{k}}\rangle,\end{split}$
(130)

hence

$\begin{split}&\langle\psi^{\prime}_{n\bm{k}}|\partial_{t}\psi^{\prime}_{0\bm{k}}\rangle=\frac{\langle\psi^{\prime}_{n\bm{k}}|\dot{\mathcal{H}}|\psi^{\prime}_{0\bm{k}}\rangle}{E^{\prime}_{0\bm{k}}-E^{\prime}_{n\bm{k}}}.\end{split}$
(131)

Substituting (127) and (131) into Eq. (123), one obtains

$\bm{J}_{\bm{k}}^{\prime h}(\eta)=\sum_{n\neq
0}\frac{\langle\psi^{\prime}_{0\bm{k}}|[\mathcal{H},\nabla-\frac{i\eta}{2}\bm{M}]|\psi^{\prime}_{n\bm{k}}\rangle\langle\psi^{\prime}_{n\bm{k}}|\dot{\mathcal{H}}|\psi^{\prime}_{0\bm{k}}\rangle}{(E^{\prime}_{0\bm{k}}-E^{\prime}_{n\bm{k}})^{3}}.$
(132)

According to the residue theorem,

$\frac{1}{(E^{\prime}_{0k}-E^{\prime}_{nk})^{3}}=-\frac{1}{2}\ointctrclockwise_{\mathcal{C}}\frac{dz}{2\pi
i}\frac{1}{(z-E^{\prime}_{0q})^{2}}\frac{1}{(z-E^{\prime}_{nq})^{2}},$

where the path $\mathcal{C}$ encircles $E^{\prime}_{0q}$. Inserting this into
the above equation, we obtain

$\bm{J}_{\bm{k}}^{\prime
h}(\eta)=-\frac{1}{2}\ointctrclockwise_{\mathcal{C}}\frac{dz}{2\pi
i}Tr\left[P_{0}g\left[\mathcal{H},\nabla-\frac{i\eta}{2}\bm{M}\right]gP_{E}g\dot{\mathcal{H}}g\right],$
(133)

where $g\equiv\frac{1}{z-\mathcal{H}}$ and
$P_{0}\equiv|\psi^{\prime}_{0\bm{k}}\rangle\langle\psi^{\prime}_{0\bm{k}}|,P_{E}\equiv\sum_{n\neq
0}|\psi^{\prime}_{n\bm{k}}\rangle\langle\psi^{\prime}_{n\bm{k}}|$. In the same
way, we obtain

$\bm{J}_{\bm{k}}^{\prime
h*}(\eta)=-\frac{1}{2}\ointctrclockwise_{\mathcal{C}}\frac{dz}{2\pi
i}Tr\left[P_{E}g\left[\mathcal{H},\nabla-\frac{i\eta}{2}\bm{M}\right]gP_{0}g\dot{\mathcal{H}}g\right].$
(134)

It is also straightforward to prove that

$\ointctrclockwise_{\mathcal{C}}\frac{dz}{2\pi
i}\left(Tr\left[P_{0}g\left[\mathcal{H},\nabla-\frac{i\eta}{2}\bm{M}\right]gP_{0}g\dot{\mathcal{H}}g\right]+Tr\left[P_{E}g\left[\mathcal{H},\nabla-\frac{i\eta}{2}\bm{M}\right]gP_{E}g\dot{\mathcal{H}}g\right]\right)=0,$

by using the residue theorem. Therefore, we can obtain the probability current
for $\mathcal{H}$

$\begin{split}\bm{J}^{\prime}_{\bm{k}}(\eta)=&\ointctrclockwise_{\mathcal{C}}\frac{dz}{2\pi
i}Tr\left[g\left[\mathcal{H},\nabla-\frac{i\eta}{2}\bm{M}\right]gg\dot{\mathcal{H}}g\right]\\\
=&\ointctrclockwise_{\mathcal{C}}\frac{dz}{2\pi
i}Tr\left[\left[g,\nabla-\frac{i\eta}{2}\bm{M}\right]g\dot{\mathcal{H}}g\right].\end{split}$
(135)

Because $\nabla-\frac{i\eta}{2}\bm{M}=-\frac{1}{2}[\mathcal{H},\bm{x}]$, we
can further simplify $\bm{J}^{\prime}_{\bm{k}}(\eta)$. However, the inclusion
of the $\bm{x}$ operator in the integrand makes the integration not well
defined, since it diverges if the size of the system is infinite.
Consequently, we modify $\bm{x}$ so that it is periodic. Suppose the periods
along the $x$ and $y$ direction are both $L$. We now define the operator
$\bm{\xi}$, with its components satisfying

$\xi_{i}=x_{i}-L\theta(x_{i})+\frac{L}{2},$ (136)

where $x_{1}=x$, $x_{2}=y$, $\theta(x)$ is the standard Heaviside function.
Thus in the region $-\frac{L}{2}<\xi_{i}\leqslant\frac{L}{2}$. Then we can
obtain

$\displaystyle\nabla-\frac{i\eta}{2}\bm{M}=\frac{1}{2}[g^{-1},\bm{\xi}]+L\bm{j}(0),$
$\displaystyle\bm{j}(0)=\frac{1}{2}[\nabla\delta(\bm{x})+\delta(\bm{x})\nabla]-\frac{i\eta}{2}\bm{M}\delta(\bm{x}).$

Replace the specific term in (135) with the above expression,

$\bm{J}^{\prime}_{\bm{k}}(\eta)=\ointctrclockwise_{\mathcal{C}}\frac{dz}{2\pi
i}\left(\frac{\partial}{\partial
t}Tr[\bm{\xi}g]+\frac{1}{2}\frac{\partial}{\partial
z}Tr[g\\{\bm{\xi},\mathcal{H}\\}]+LTr[[g,j(0)]g\dot{\mathcal{H}}g]\right).$
(137)

The first term turns out to be zero after we take the average of it over time.
The second term is zero due to the periodicity of the path. The last term is
the only one that contributes. We write it in a more explicit form

$\begin{split}\bm{J}^{\prime}_{\bm{k}}(\eta)=&L\ointctrclockwise_{\mathcal{C}}\frac{dz}{2\pi
i}\iiint_{-L/2}^{L/2}d\bm{x}d\bm{x}^{\prime}d\bm{x}^{\prime\prime}[g(\bm{x}^{\prime},\bm{x})\bm{j}(0)g(\bm{x},\bm{x}^{\prime\prime})\dot{\mathcal{H}}(\bm{x}^{\prime\prime})g(\bm{x}^{\prime\prime},\bm{x}^{\prime})\\\
&-g(\bm{x}^{\prime},\bm{x})\bm{j}(0)g(\bm{x},\bm{x}^{\prime\prime})g(\bm{x}^{\prime\prime},\bm{x}^{\prime})\dot{\mathcal{H}}(\bm{x}^{\prime})].\end{split}$
(138)

The main analogy with Ref. [15] is that the single particle Hamiltonian
$h=-\frac{1}{2}(d/dx)^{2}+U(x,\tau)$ and the Hermitian operator
$\mathcal{H}=L_{H}-i\eta L_{A}=-\nabla^{2}+i\eta\bm{M}\cdot\nabla+U^{\prime}$
are similar. They both possess a kinetic term and a potential term. In each
case, in the deterministic limit, the potential term dominates, consequently
the Green functions $g(x,x^{\prime})$ decays exponentially if $x-x^{\prime}$
deviates from the peaks.

As the potential is periodic, the eigenfucntions of of $\mathcal{H}$ are Bloch
waves, which are superpositions of the Wanner functions, which are localized.
The eigenfunctions can be written as

$\begin{split}\psi_{\bm{r}_{0},\bm{k}}=&\sum_{m,n\in\mathds{Z}}e^{i(k_{x}mL_{x}+k_{y}nL_{y})}\Gamma_{x}(x-x_{0}-mL_{x})\Gamma_{y}(y-y_{0}-nL_{y}),\end{split}$
(139)

where $\bm{r}_{0}=(x_{0},y_{0})$, where $x_{0}\in[0,L_{x})$ and
$y_{0}\in[0,L_{y})$, can be regarded as the band index, $m$ and $n$ are
integers, $\Gamma_{x}(x-x_{0})$ is a localized function peaked at $x_{0}$,
$\Gamma_{y}(y-y_{0})$ is a localized function peaked at $y_{0}$.

Consequently, the Green functions can be calculated as follows

$\begin{split}g_{\alpha}(\bm{x},\bm{x^{\prime}})=&\langle\bm{x}|\frac{1}{z-\mathcal{H}}|\bm{x^{\prime}}\rangle\\\
=&\iint
d^{2}\bm{r}\langle\bm{x}|\frac{1}{z-\mathcal{H}}|\psi_{\bm{r},\bm{\alpha}}\rangle\langle\psi_{\bm{r},\bm{\alpha}}|\bm{x}^{\prime}\rangle.\end{split}$
(140)

In the deterministic limit or low-temperature limit, $\mathcal{H}\rightarrow
U^{\prime}$, therefore

$\begin{split}g_{\alpha}(\bm{x},\bm{x^{\prime}})\rightarrow&\iint
d^{2}\bm{r}\frac{\psi_{\bm{r},\bm{\alpha}}(\bm{x})\psi^{*}_{\bm{r},\bm{\alpha}}(\bm{x}^{\prime})}{z-U^{\prime}(\bm{r})}\\\
=&\sum_{m,n,m^{\prime},n^{\prime}}e^{ik_{x}(m-m^{\prime})L_{x}+ik_{y}(n-n^{\prime})L_{y}}\iint
d^{2}\bm{r}\frac{F(x,x^{\prime},y,y^{\prime})}{z-U^{\prime}(\bm{r})},\end{split}$
(141)

where
$F(x,x^{\prime},y,y^{\prime})\equiv\Gamma_{x}(x-x_{0}-mL_{x})\Gamma_{y}(y-y_{0}-nL_{y})\Gamma_{x}(x^{\prime}-x_{0}-m^{\prime}L_{x})\Gamma_{y}(y^{\prime}-y_{0}-n^{\prime}L_{y})\approx\Gamma_{x}(x-x^{\prime}-(m-m^{\prime})L_{x})\Gamma_{y}(y-y^{\prime}-(n-n^{\prime})L_{y})$,
therefore

$\begin{split}g_{\alpha}(\bm{x},\bm{x^{\prime}})\rightarrow&\frac{1}{z-U(\bm{x}^{\prime})}\sum_{m,n,m^{\prime},n^{\prime}}e^{ik_{x}(m-m^{\prime})L_{x}+ik_{y}(n-n^{\prime})L_{y}}\Gamma_{x}(x-x^{\prime}-(m-m^{\prime})L_{x})\Gamma_{y}(y-y^{\prime}-(n-n^{\prime})L_{y})\\\
=&\frac{\mathcal{N}}{z-U(\bm{x}^{\prime})}\sum_{s,t}e^{ik_{x}sL_{x}+ik_{y}tL_{y}}\Gamma_{x}(x-x^{\prime}-sL_{x})\Gamma_{y}(y-y^{\prime}-tL_{y}),\end{split}$
(142)

where $\mathcal{N}$ is the number of different values of $(m,n)$.

It is clear that in the deterministic or low-temperature limit, the Green
functions peak at points with $|x-x^{\prime}|=sL_{x}$ and
$|y-y^{\prime}|=tL_{y}$, and decay rapidly away from the the peaks.

On the other hand, $\bm{j}(0)$ in (138) contains Dirac delta functions centred
at $\bm{x}=0$, therefore the integrand is considerable only when
$\bm{x}=\bm{x}^{\prime}=\bm{x}^{\prime\prime}=0$, i.e. $s=t=0$. In this case,
in the Green Functions as given in (142), the $\mathbf{k}$-dependent terms,
only appearing as the exponents tend to vanish, consequently the Green
functions and thus $\bm{J}^{\prime}_{\bm{k}}(\eta)$ are insensitive to
$\bm{k}$ in the deterministic or low-temperature limit.

Finally we consider the analytical continuation of $\eta$ to $i$ [18]. Then
$\mathcal{H}=\tilde{\mathcal{O}}$ and
$\bm{J}^{\prime}_{\bm{k}}(i)=\bm{J}_{\bm{k}}$. Since
$\bm{J}^{\prime}_{\bm{k}}(\eta)$ is also insensitive to $\bm{k}$, so is
$\bm{J}^{\prime}_{\bm{k}}(i)$. Hence we can arrive at the conclusion that
$\bm{J}_{\bm{k}}$ is insensitive to $\bm{k}$.

On the other hand, our simulation results confirm topological quantization,
hence indirectly confirm the the insensitivity of $\bm{J}_{\bm{k}}$ to
$\bm{k}$, consistent with the validity of the analytic continuation.

##  Appendix C Simplification of $\bm{J}_{\bm{k}}^{h}$

After the introduction of the Bloch periodic function, $\bm{J}_{\bm{k}}^{h}$
can be transformed

$\displaystyle\bm{J}_{\bm{k}}^{h}=$ $\displaystyle\sum_{n\neq 0}\frac{\langle
e^{i\bm{k}\cdot\bm{r}}v_{0\bm{k}}|\left[\nabla+\frac{1}{2}\rho_{0}^{-1}(\hat{\bm{e}}_{x}\partial_{y}A-\hat{\bm{e}}_{y}\partial_{x}A)\right]|e^{i\bm{k}\cdot\bm{r}}w_{n\bm{k}}\rangle\langle
e^{i\bm{k}\cdot\bm{r}}v_{n\bm{k}}|e^{i\bm{k}\cdot\bm{r}}\dot{w}_{0\bm{k}}\rangle}{E_{0\bm{k}}-E_{n\bm{k}}}$
$\displaystyle=$ $\displaystyle\sum_{n\neq 0}\frac{\langle
v_{0\bm{k}}|\left[\nabla+i\bm{k}+\frac{1}{2}\rho_{0}^{-1}(\hat{\bm{e}}_{x}\partial_{y}A-\hat{\bm{e}}_{y}\partial_{x}A)\right]|w_{n\bm{k}}\rangle\langle
v_{n\bm{k}}|\dot{w}_{0\bm{k}}\rangle}{E_{0\bm{k}}-E_{n\bm{k}}}.$

In the Hilbert space of $w_{n\bm{k}}$, the transformed Fokker-Planck operator
must be transformed to

$\tilde{\mathcal{O}}^{\prime}=-(\nabla+i\bm{k})^{2}-\rho_{0}^{-1}(\hat{\bm{e}}_{x}\partial_{y}A-\hat{\bm{e}}_{y}\partial_{x}A)\cdot(\nabla+i\bm{k})+U,$
(143)

in order that
$\tilde{\mathcal{O}}^{\prime}w_{n\bm{k}}=E_{n\bm{k}}w_{n\bm{k}}$. Then
calculate the derivative of $\tilde{\mathcal{O}}^{\prime}$ versus $\bm{k}$

$\displaystyle\frac{\partial\tilde{\mathcal{O}}^{\prime}}{\partial\bm{k}}=-2i(\nabla+i\bm{k})-i\rho_{0}^{-1}(\hat{\bm{e}}_{x}\partial_{y}A-\hat{\bm{e}}_{y}\partial_{x}A)$
$\displaystyle\Rightarrow$
$\displaystyle\nabla+i\bm{k}+\frac{1}{2}\rho_{0}^{-1}(\hat{\bm{e}}_{x}\partial_{y}A-\hat{\bm{e}}_{y}\partial_{x}A)=\frac{i}{2}\frac{\partial\tilde{\mathcal{O}}^{\prime}}{\partial\bm{k}}.$

Consequently, $\bm{J}_{\bm{k}}^{h}$ can be simplified,

$\bm{J}_{\bm{k}}^{h}=\frac{i}{2}\sum_{n\neq 0}\frac{\langle
v_{0\bm{k}}|\frac{\partial\tilde{\mathcal{O}}^{\prime}}{\partial\bm{k}}|w_{n\bm{k}}\rangle}{E_{0\bm{k}}-E_{n\bm{k}}}\langle
v_{n\bm{k}}|\dot{w}_{0\bm{k}}\rangle.$ (144)

After that, calculate the derivative of
$\tilde{\mathcal{O}}^{\prime}w_{n\bm{k}}=E_{n\bm{k}}w_{n\bm{k}}$ versus
$\bm{k}$,

$\frac{\partial\tilde{\mathcal{O}}^{\prime}}{\partial\bm{k}}w_{n\bm{k}}+\tilde{\mathcal{O}}^{\prime}\frac{\partial
w_{n\bm{k}}}{\partial\bm{k}}=\frac{\partial
E_{n\bm{k}}}{\partial\bm{k}}w_{n\bm{k}}+E_{n\bm{k}}\frac{\partial
w_{n\bm{k}}}{\partial\bm{k}}.$

Then take the inner product of $v_{0\bm{k}}$ and the above equation,

$\left\langle
v_{0\bm{k}}\left|\frac{\partial\tilde{\mathcal{O}}^{\prime}}{\partial\bm{k}}\right|w_{n\bm{k}}\right\rangle+\left\langle\tilde{\mathcal{O}}^{\prime{\dagger}}v_{0\bm{k}}\left|\frac{\partial
w_{n\bm{k}}}{\partial\bm{k}}\right.\right\rangle=\frac{\partial
E_{n\bm{k}}}{\partial\bm{k}}\langle
v_{0\bm{k}}|w_{n\bm{k}}\rangle+E_{n\bm{k}}\left\langle
v_{0\bm{k}}\left|\frac{\partial
w_{n\bm{k}}}{\partial\bm{k}}\right.\right\rangle.$

It is straightforward to obtain
$\tilde{\mathcal{O}}^{\prime{\dagger}}v_{0\bm{k}}=E_{0\bm{k}}^{*}v_{0\bm{k}}$
and $\langle v_{0\bm{k}}|w_{n\bm{k}}\rangle=0$ for $n\neq 0$, which lead to

$\frac{\langle
v_{0\bm{k}}|\frac{\partial{\tilde{\mathcal{O}}^{\prime}}}{\partial\bm{k}}|w_{n\bm{k}}\rangle}{E_{0\bm{k}}-E_{n\bm{k}}}=-\left\langle
v_{0\bm{k}}\left|\frac{\partial
w_{n\bm{k}}}{\partial\bm{k}}\right.\right\rangle=\left\langle\left.\frac{\partial
v_{0\bm{k}}}{\partial\bm{k}}\right|w_{n\bm{k}}\right\rangle.$

By substituting the above equation into (144), one obtains

$\bm{J}_{\bm{k}}^{h}=\frac{i}{2}\sum_{n\neq 0}\left\langle\left.\frac{\partial
v_{0\bm{k}}}{\partial\bm{k}}\right|w_{n\bm{k}}\right\rangle\langle
v_{n\bm{k}}|\dot{w}_{0\bm{k}}\rangle.$ (145)

One has a completeness identity
$\mathbf{1}=\sum_{n}|\psi_{n}\rangle\langle\phi_{n}|=\sum_{n}|w_{n}\rangle\langle
v_{n}|$ , where the index $\bm{k}$ is omitted for simplicity. Its validity can
be justified by the calculation of the matrix elements,

$\displaystyle\langle\psi_{i}|\left[\sum_{n}|\psi_{n}\rangle\langle\phi_{n}|\right]|\psi_{j}\rangle$
$\displaystyle=\sum_{n}T_{in}\delta_{nj}=T_{ij}=\langle\psi_{i}|\psi_{j}\rangle,$
$\displaystyle\langle\phi_{i}|\left[\sum_{n}|\psi_{n}\rangle\langle\phi_{n}|\right]|\psi_{j}\rangle$
$\displaystyle=\sum_{n}\delta_{in}\delta_{nj}=\delta_{ij}=\langle\phi_{i}|\psi_{j}\rangle.$

Therefore, half of the probability current can be simplified further

$\begin{split}\bm{J}_{\bm{k}}^{h}=&\frac{i}{2}\langle\partial_{\bm{k}}v_{0\bm{k}}|\partial_{t}w_{0\bm{k}}\rangle-\frac{i}{2}\langle\partial_{\bm{k}}v_{0\bm{k}}|w_{0\bm{k}}\rangle\langle
v_{0\bm{k}}|\partial_{t}w_{0\bm{k}}\rangle.\end{split}$ (146)

## References

  * [1] U. K. Rößler, N. Bogdanov. and C. Pfleiderer,  Spontaneous skyrmion ground states in magnetic metals.  Nature(London) 442, 797 (2006).
  * [2] S. Mühlbauer, B. Binz, F. Jonietz, C. Pfleiderer, A. Rosch, A. Neubauer, R. Georgii, and P. Böni,  Skyrmion lattice in a chiral magnet.  Science 323, 915 (2009).
  * [3] X. Z. Yu, Y. Onose, N. Kanazawa, J. H. Park, J. H. Han, Y. Matsui, N. Nagaosa, and Y. Tokura,  Real-space observation of a two-dimensional skyrmion crystal.  Nature 465, 901 (2010).
  * [4] X. Z. Yu, N. Kanazawa, Y. Onose, K. Kimoto, W. Z. Zhang, S. Ishiwata, Y. Matsui, and Y. Tokura,  Near room-temperature formation of a skyrmion crystal in thin-films of the helimagnet FeGe.  Nat. Mater. 10, 106 (2011).
  * [5] S. Heinze, K. von Bergmann, M. Menzel, J. Brede, A. Kubetzka, R. Wiesendanger, G. Bihlmayer, and S. Blüegel,  Spontaneous atomic-scale magnetic skyrmion lattice in two dimensions.  Nat. Phys. 7, 713 (2011).
  * [6] N. Nagaosa and Y. Tokura,  Topological properties and dynamics of magnetic skyrmions.  Nat. Nanotech. 8, 899 (2013).
  * [7] P. Reimann, Brownian motors: noisy transport far from equilibrium, Phys. Rep. 361, 57 (2002).
  * [8] Y. Shi and Q. Niu,  Quantization and corrections of adiabatic particle transport in a periodic ratchet potential.  EPL 59, 324 (2002).
  * [9] S.-C. Tang and Y. Shi  Topological quantization of the flow of magnetic skyrmions driven by a ratchet-like potential under thermal fluctuations.  arXiv:1806.05354 (2018).
  * [10] M. Mochizuki, X. Z. Yu, S. Seki, N. Kanazawa, W. Koshibae, J. Zang, M. Mostovoy, Y. Tokura, and N. Nagaosa,  Thermally driven ratchet motion of a skyrmion microcrystal and topological magnon Hall effect.  Nat. Mater. 13, 241 (2014).
  * [11] C. Reichhardt, D. Ray, and C. J. Olson Reichhardt,  Magnus-induced ratchet effects for skyrmions interacting with asymmetric substrates.  New J. Phys. 17, 073034 (2015).
  * [12] X. Ma, C. J. Olson Reichhardt, and C. Reichhardt,  Reversible vector ratchets for skyrmion systems.  Phys. Rev. B 95, 104401 (2017).
  * [13] W. Chen, L. Liu, Y. Ji, and Y. Zheng,  Skyrmion ratchet effect driven by a biharmonic force.  Phys. Rev. B 99, 064431 (2019).
  * [14] D. J. Thouless,  Quantization of particle transport,  Phys. Rev. B 27, 6083 (1983)
  * [15] Q. Niu and D. J. Thouless,  Quantised adiabatic charge transport in the presence of substrate disorder and many-body interaction.  J. Phys. A: Math. Gen. 17, 2453 (1984).
  * [16] H. Ochoa and Y. Tserkovnyak,  Quantum skyrmionics,  Int. J. Mod. Phys. B 33, 1930005 (2019)
  * [17] R. E. Troncoso and Á. S. Núñez,  Thermally assisted current-driven skyrmion motion.  Phys. Rev. B 89, 224403 (2014).
  * [18] H. Risken,  The Fokker-Planck Equation (Springer, Berlin, 1989).
  * [19] J. Iwasaki, M. Mochizuki, and N. Nagaosa,  Universal current-velocity relation of skyrmion motion in chiral magnets.  Nat. Commun. 4, 1463 (2013).
  * [20] L. Kong and J. Zang,  Dynamics of an insulating skyrmion under a temperature gradient.  Phys. Rev. Lett. 111, 067203 (2013).
  * [21] S. Lin, C. D. Batista, C. Reichhardt, and A. Saxena,  AC current generation in chiral magnetic insulators and skyrmion motion induced by the spin seebeck effect.  Phys. Rev. Lett. 112, 187203 (2014).
  * [22] R. E. Troncoso and Á. S. Núñez,  Brownian motion of massive skyrmions in magnetic thin films.  Annals of Physics 351, 850 (2014).
  * [23] J. L. García-Palacios and F. J. Lázaro,  Langevin-dynamics study of the dynamical properties of small magnetic particles.  Phys. Rev. B 58, 14937 (1998).
  * [24] I. Dzyaloshinsky,  A thermodynamic theory of ”weak” ferromagnetism of antiferromagnetics.  Journal of Physics and Chemistry of Solids 4, 241 (1958).
  * [25] T. Moriya,  Anisotropic superexchange interaction and weak ferromagnetism.  Phys. Rev. 120, 91 (1960).
  * [26] R. Bartussek, P. Hänggi, and J. G. Kissner,  Periodically rocked thermal ratchets.  EPL 28, 459 (1994).