# Microscopic theory of current-induced skyrmion transport and its application
in disordered spin textures

Emil Östberg<EMAIL_ADDRESS>Department of Physics, Division of
Mathematical Physics, Lund University, 22100 Lund, Sweden  Emil Viñas Boström
<EMAIL_ADDRESS>Max Planck Institute for the Structure and Dynamics
of Matter, Luruper Chaussee 149, 22761 Hamburg, Germany  Nano-Bio Spectroscopy
Group, Departamento de Fisica de Materiales, Universidad del Pais Vasco, 20018
San Sebastian, Spain  Claudio Verdozzi<EMAIL_ADDRESS>Department of Physics, Division of Mathematical Physics and ETSF, Lund
University, 22100 Lund, Sweden

###### Abstract

Magnetic skyrmions hold great promise for realizing compact and stable memory
devices that can be manipulated at very low energy costs via electronic
current densities. In this work, we extend a recently introduced method to
describe classical skyrmion textures coupled to dynamical itinerant electrons.
In this scheme, the electron dynamics is described via nonequilibrium Green’s
functions (NEGF) within the generalized Kadanoff-Baym ansatz, and the
classical spins are treated via the Landau-Lifshitz-Gilbert equation. The
framework is here extended to open systems, by the introduction of a non-
interacting approximation to the collision integral of NEGF. This, in turn,
allows us to perform computations of the real-time response of skyrmions to
electronic currents in large quantum systems coupled to electronic reservoirs,
which exhibit a linear scaling in the number of time steps. We use this
approach to investigate how electronic spin currents and dilute spin disorder
affects skyrmion transport and the skyrmion Hall drift. Our results show that
the skyrmion dynamics is sensitive to the specific form of spin disorder, such
that different disorder configurations leads to qualitatively different
skyrmion trajectories for the same applied bias. This sensitivity arises from
the local spin dynamics around the magnetic impurities, a feature that is
expected not to be well captured by phenomenological or spin-only
descriptions. At the same time, our findings illustrate the potential of
engineering microscopic impurity patterns to steer skyrmion trajectories.

##  I Introduction

Technological progress often stems from discovery and exploitation of new
materials and forms of energy. While self-evident, this paradigm has recently
undergone criticism and revision, due to mounting awareness of the negative
impact that indiscriminate technological development has on environment and
climate. This is also true for electronics: it has become clear that
production, use and casual disposal of electronic devices can lead to a sharp
increase in energy consumption, waste, and greenhouse effects Waste . Thus,
together with a steady increase in the use of high performance technology,
there is a need for novel electronics with reduced dimensionality, large
integration and low energy consumption Green .

Pursuing on equal footing these two directives is the core aim of spintronics
Spintronics1 : magnetic excitations allow for less energy-intensive ways of
storing and processing digital information, and thus devices based on magnetic
materials and phenomena offer an attractive alternative to conventional
electronics. For a long time, fundamental and applied research in magnetism
was largely concerned with macroscopic samples, and primarily with simple
magnetic orders, such as ferro- and antiferromagnets. However, more recently,
it has been possible to experimentally realize magnetic systems with non-
trivial magnetic textures, creating unprecedented possibilities for spintronic
applications Spintronics2 ; Parkin2015 .

A notable example in this respect is provided by magnetic skyrmions Mertig2021
. These are topologically nontrivial spin textures stabilized by a competition
of exchange, Dzyaloshinskii-Moriya interactions (DMIs) and magnetic
anisotropies Skyrmionics . With their compact size, topological protection,
and non-intensive energy requirements for manipulation, skyrmions are of great
potential interest to realize race track memories Sampaio2013 ; Nagaosa ;
Romming2015  and quantum computation devices Yang16 ; Chauwin19 ; Mirlin19 ;
Psaroudaki2021 . This, however, requires efficient ways of writing, deleting,
and manipulating skyrmions on short time scales and with high spatial
precision, via electronic spin currents.

Figure 1: A sketch of the composite spin-electron system. A central region
consisting of itinerant electrons interacting with localized spins is in
contact with two-dimensional leads. By applying a spin dependent bias in the
leads an electronic current is generated in the central system giving rise to
skyrmion motion.

Here we describe a scheme to perform large scale simulations of the
intertwined dynamics of interacting and open spin-electron systems (Fig. 1).
Our method is an extension of the approach introduced in Refs. LeadsPaper ;
LaserPaper , and amounts to propagating the equation of motion for the
electronic spin-dependent single-particle density matrix together with the
Landau-Lifshitz-Gilbert equation for the classical spins. The former equation
can be derived from the general theory of non-equilibrium Green’s functions
using the so-called generalized Kadanoff-Baym ansatz (GKBA) Lipavsky86 ;
Hermanns12 ; Hermanns14 ; Latini14 ; Bostrom18 ; Hopjan18 ; KarlssonHopjan18 ;
Kalvova18 ; Joost20 ; Karlsson21 ; Pavlyukh22a ; Pavlyukh22b ; BonitzWBL23 ;
Tuovinen23 , and therefore allows to systematically introduce the effects of
electron-electron interactions via diagrammatic many-body perturbation theory.
In what follows, we apply this scheme to study current-induced skyrmion
motion, fully accounting for the dynamics of the itinerant electrons resulting
from an applied bias. In agreement with phenomenological theories Iwasaki13 ,
we find that for a clean sample skyrmions are pinned below a critical spin
current density $I_{0}$, after which the velocity is a linear function of the
current $I-I_{0}$. The situation is found to be qualitatively different in
presence of (dilute) spin disorder, where the skyrmion motion is strongly
dependent on the location, size and form of the disorder configuration (for
previous work on the role of disorder, see e.g.Thorwart2021 ). This provides a
clear indication that treatments based on the standard Thiele or Landau-
Lifshitz-Gilbert (LLG) equation are not always adequate, and that the dynamics
of the electrons must be explicitly taken into account.

The manuscript is structured as follows: In Sec. II we briefly review previous
approaches to coupled spin-electron dynamics. In Sec. III we introduce the
system and Hamiltonian to be considered, and discuss the coupling between the
itinerant electrons and external reservoirs. The spin and electron equations
of motion are presented in Sec. IV and Sec. V, and Sec. VI introduces an
approximate wide band limit (AWBL) as a numerically efficient way to propagate
the equations of motion in presence of large central regions connected to
external reservoirs. In Sec. VII we present some observables used to
interrogate the skyrmion content of the spin configuration. In Sec. VIII we
use the AWBL to investigate skyrmion motion induced by current densities in
the itinerant electron system, and in Sec. IX we consider skyrmion motion in
presence of magnetic disorder. Finally, in Sec. X we conclude with a
discussion of experimental signatures and possible material platforms for
which our results are of relevance.

##  II Review of previous approaches

Magnetic skyrmions are made up of localized magnetic moments, typically
arising from a large Hund’s coupling $J$ favoring a high-spin state the $d$\-
or $f$-orbitals of the magnetic ion. Therefore, from a microscopic
perspective, it is natural to expect that a quantum description of the
magnetic structure should be necessary 111For current work in this direction,
see e.g. Haller22 . However, since magnetic moments in typical skyrmion
materials (consisting of transition metal ions such as Fe, Co and Mn) are of
large magnitude, spin fluctuations are suppressed and a classical
approximation usually works well.

In most cases skyrmion textures are large compared to the underlying lattice
constant. More precisely, the ratio of the electronic hopping $t$ and the
strength $\alpha$ of the spin-orbit interaction is typically of the order
$\alpha/t\sim 0.01-0.1$, resulting in a spin spiral wavelength $\lambda\sim
10-100$ nm Nagaosa ; Experimental_review . In certain cases however, such as
at an interface between metallic thin films and a material with large spin-
orbit coupling, the effective Dzyaloshinksii-Moriya interaction can be
significantly enhanced and lead to skyrmions with radii on the order of a few
nm Heinze11 . For skyrmions of large sizes it is common to take a continuum
limit of the microscopic spin Hamiltonian, resulting in a magnetic energy
functional that can be minimized with micromagnetic methods. As a result,
calculations for lattice skyrmions of realistic size (up to about $10$ nm in
radius in two-dimensional systems) are usually based on an atomistic
description of classical spins, with external fields such electromagnetic
radiation or electronic current densities included as non-dynamics variables.
For example, the motion of skyrmions in response to an external current
density can been described via an effective equation for the skyrmion center
of mass, the so-called Thiele equation Thiele73 , assuming that the form of
the skyrmion is rigid. In the Thiele equation, electrons enter only via the
external current density, taken from a static solution of the macroscopic
Maxwell equations. A more detailed description often considered in the
literature is to obtain the individual spin dynamics from the Landau-Lifshitz-
Gilbert (LLG) equation Iwasaki13 , and to include the effects of electrons and
external fields via a generalized force Nikolic1 ; Nikolic2 . The LLG
equation, with itinerant electrons included implicitly, has provided important
insights into skyrmion behavior in a large range of materials, and its use is
widespread. Furthermore, other approaches have been introduced, besides the
LLG, that go beyond the semiclassical Thiele’s description to study the
dynamics of skyrmions PhysRevLett.127.047203 ; PhysRevB.107.224418 ;
PhysRevB.107.174418 ; PhysRevB.108.064427 .

Indeed, there are many cases where neglecting the explicit dynamics of the
electrons severely hinders a more detailed understanding, and possibly even
the development of novel physical ideas and technological opportunities in
skyrmionics.

The importance of explicitly accounting for itinerant electrons in the
description of skyrmion dynamics was originally pointed out in Ref. LeadsPaper
, where a small two-dimensional spin texture containing a single skyrmion was
made to interact with a nanowire carrying a time-dependent current. Since
then, the significance of including electronic degrees of freedom in the
dynamical description of skyrmions (as well as more general spin textures),
has been addressed in several contexts LaserPaper ; bluegel1 ; Sinova ;
Bluegel2 ; Nikolic1 ; Nikolic2 . For example, it was recently demonstrated
LaserPaper  that explicitly accounting for the dynamics of itinerant electrons
can be of crucial importance in capturing the interaction between skyrmions
and laser light. In particular, such a treatment shows that photo-generation
of skyrmions resulting from laser excitation can occur on much shorter
timescales than previously thought.

Similarly, electrons and currents are expected to play a key role in skyrmion
transport phenomena both for clean samples and when disorder is present.
However, to investigate such dynamics requires to overcome an additional
hurdle on the methodological side, namely that it is necessary to
simultaneously treat large system sizes and extended time scales within an
open quantum system framework. In the remainder of this manuscript, we show
how such a framework can be constructed.

##  III System and Hamiltonian

We consider a system at zero temperature and of finite size, referred to as
the central region $C$, in contact with two macroscopic reservoirs (see Fig.
2). The central region consists of a finite square lattice with $N=N_{x}\times
N_{y}$ sites, and at each lattice site there is one electronic orbital and one
localized spin. The electronic orbitals are populated by spinful itinerant
electrons that can tunnel in and out of the reservoirs, also referred to as
leads. The leads couple only to the electronic degrees of freedom, i.e. there
are no localized spins in the reservoirs. Inside the leads the electrons are
assumed to be non-interacting and spin polarized, with a density set by the
chemical potential. Also, in the central region the electrons are assumed to
be non-interacting, although electron-electron interactions are
straightforward (but numerically demanding) to include LeadsPaper . However,
the itinerant electrons of the central region interact with the localized
magnetic moments via a local exchange coupling, and the magnetic moments in
turn interact among themselves via various magnetic interactions. The
Hamiltonian of the composite system is

$\displaystyle H(t)=H_{C}(t)+H_{R}(t)+H_{CR},$ (1)

where a possible time-dependence is explicitly indicated both for the central
region and the reservoirs. This is necessary to initiate the dynamics of the
system. We now discuss each contribution to $H$ separately.

###  III.1 The central region

The Hamiltonian of the central region is $H_{C}(t)=H_{e}+H_{s}+H_{s-e}$, where
$H_{e}$ describes the itinerant electrons, $H_{s}$ the localized spins, and
$H_{s-e}$ the spin-electron coupling. These terms are respectively given by

$\displaystyle H_{e}$
$\displaystyle=\sum_{\expectationvalue{ij}\sigma\sigma^{\prime}}c_{i\sigma}^{\dagger}(-t_{ij}\mathbf{1}+\bm{\alpha}_{ij}\cdot\bm{\tau})_{\sigma\sigma^{\prime}}c_{j\sigma^{\prime}}-B\sum_{i}\hat{s}_{i}^{z}$
(2) $\displaystyle H_{s}$
$\displaystyle=-\frac{1}{2}\sum_{\expectationvalue{ij}}J_{ij}\hat{\mathbf{S}}_{i}\cdot\hat{\mathbf{S}}_{j}-\frac{1}{2}\sum_{\expectationvalue{ij}}\mathbf{D}_{ij}\cdot\quantity(\hat{\mathbf{S}}_{i}\crossproduct\hat{\mathbf{S}}_{j})-B\sum_{i}\hat{S}_{i}^{z}-\frac{K}{2}\sum_{i}(\hat{S}_{i}^{z})^{2}$
(3) $\displaystyle H_{s-e}$
$\displaystyle=-g\sum_{i}\hat{\mathbf{S}}_{i}\cdot\hat{\mathbf{s}}_{i}.$ (4)

Here $c^{\dagger}_{i\sigma}$ creates an itinerant electron at site $i$ with
spin projection $\sigma$, the hopping amplitude between nearest-neighbor sites
$i$ and $j$ (denoted by $\expectationvalue{ij}$) is given by $t_{ij}$,
$\bm{\alpha}_{ij}$ accounts for spin-orbit interactions, and ${\bf
B}=B\hat{\bf z}$ is an external magnetic field along the $z$-axis. The
electronic spin operator at site $i$ is defined by $\hat{\bf
s}_{i}=\sum_{\sigma\sigma^{\prime}}c_{i\sigma}^{\dagger}\bm{\tau}_{\sigma\sigma^{\prime}}c_{i\sigma^{\prime}}$,
where $\bm{\tau}$ denotes the vector of Pauli matrices. The parameters
$J_{ij}=J_{ji}$ and ${\bf D}_{ij}=-{\bf D}_{ji}$ give the exchange and
Dzyaloshinskii-Moriya interaction (DMI) Dzyaloshinskii58 ; Moriya60  between
spins $\hat{\bf S}_{i}$ and $\hat{\bf S}_{j}$, and $K$ quantifies an easy-axis
single-ion anisotropy. The itinerant electrons interact with localized spins
via a local exchange coupling of strength $g$. We note that in magnetic thin
films, skyrmions are often stabilized by the easy-axis anisotropy Heinze11 ,
and an external magnetic field is not strictly necessary. However, including a
Zeeman term in the Hamiltonian provides an additional physical means to
control the equilibrium magnetic state, and in particular to tune the system
between a ferromagnetic state and the skyrmion crystal phase.

Depending on the symmetries of the system, the DMI vector ${\bf D}_{ij}$ can
be of different form. Assuming the overall magnitude of the interaction is
fixed at $|{\bf D}_{ij}|=D$, the spatial dependence typically takes one of the
following two forms

$\displaystyle\mathbf{D}_{ij}$ $\displaystyle=D\,\mathbf{d}_{ij}\hskip
48.36958pt\text{(Bloch)}$ (5) $\displaystyle\mathbf{D}_{ij}$
$\displaystyle=D\,\mathbf{d}_{ij}\times\hat{\mathbf{z}}.\hskip
28.45274pt\text{(N\'{e}el)}$

Here $\mathbf{d}_{ij}$ is the vector between lattice sites $i$ and $j$, and
$\hat{\mathbf{z}}$ is a normal vector to the two-dimensional lattice plane.
The Néel-type DMI gives rise to Néel-type (hedgehog) skyrmions, and is
typically generated by the inversion symmetry breaking induced at a surface.
The Néel type DMI is therefore most common in single layer or few layer
substrates, since thicker materials tend to restore the bulk inversion
symmetry. On the other hand, the Bloch-type DMI that commonly arises in
systems with non-centrosymmetric crystal structures, gives rise to Bloch-type
(spiral) skyrmions. Even so, many bulk materials that support Bloch-type
skyrmions can be fabricated as thin film samples, while still hosting
skyrmions Experimental_review .

Figure 2: Schematic illustration of the transport setup and the relevant
parameters of the system-reservoir Hamiltonian. The leads have a common
bandwidth $W$ and chemical potential $\mu$, and are connected to a central
region of $20\times 40$ sites comprising both electrons and localised spins
via a nearest-neighbour tunneling amplitude $t_{L/R}$ at the left/right edge.
To initiate the skyrmion dynamics, a spin-dependent bias
($V_{L\uparrow}=-V_{R\uparrow}$ and $V_{L\downarrow}=V_{L\downarrow}=0$) is
applied in the leads.

In this work, we will limit ourselves to a Bloch-type DMI, while considering
for computational simplicity a monolayer geometry for the localized spins.
Similarly we will focus on a Rashba-type spin-orbit coupling of the form

$\displaystyle\bm{\alpha}_{ij}=i\alpha_{R}\,\mathbf{d}_{ij}\times\hat{\mathbf{z}},$
(6)

where the parameter $\alpha_{R}$ sets the overall magnitude of the
interaction. This reflects the fact that the forms of the spin-orbit coupling
and Dzyaloshinskii-Moriya interaction are intimately related, such that to
lowest order in a strong coupling expansion ${\bf
D}_{ij}\propto\bm{\alpha}_{ij}$ Moriya60 . We note that the Hamiltonian
$H_{C}$ can be generalized further by including both a local on-site potential
$v_{i}(t)$ for the itinerant electrons, as well as a complex hopping
$t_{ij}e^{i\phi_{ij}}$ whose phase encodes the interaction with an external
electromagnetic field LaserPaper . It is similarly straightforward to include
a direct coupling between the localized spins and external electromagnetic
fields e.g. via the inverse Faraday effect. For the present work such terms
are not of relevance, and are therefore omitted.

###  III.2 The reservoirs and their coupling to the central region

The leads are taken to be two-dimensional and semi-infinite (see Fig. 1), and
are described by the Hamiltonian

$\displaystyle H_{R}(t)$
$\displaystyle=\sum_{\alpha\sigma}H^{R}_{\alpha,\sigma}(t)$ (7) $\displaystyle
H^{R}_{\alpha\sigma}(t)$
$\displaystyle=t_{R\alpha}\sum_{\expectationvalue{ij}\sigma}a_{i\sigma,\alpha}^{\dagger}a_{j\sigma,\alpha}+\sum_{i\sigma}v^{\alpha}_{\sigma}(t)a_{i\sigma,\alpha}^{\dagger}a_{i\sigma,\alpha}.$
(8)

Here $\alpha\in\\{L,R\\}$ denotes the left (L) and right (R) lead, $\sigma$
labels the spin projection, and the parameter $t_{R\alpha}$ is the hopping
amplitude within lead $\alpha$. The hopping amplitude is related to the
nominal band width by $W_{\alpha}=4t_{R\alpha}$, and $v^{\alpha}_{\sigma}(t)$
is a time-dependent bias measured from the chemical potential $\mu$. The leads
interact with the central region via the Hamiltonian
$H_{CR}=\sum_{\alpha\sigma}H^{CR}_{\alpha,\sigma}$, where

$\displaystyle
H^{CR}_{L\sigma}=t_{L}\sum_{i_{y}=1}^{N_{y}}\big{[}a^{\dagger}_{(1,i_{y})\sigma,L}c_{(1,i_{y})\sigma}+H.c.\big{]}$
(9) $\displaystyle
H^{CR}_{R\sigma}=t_{R}\sum_{i_{y}=1}^{N_{y}}\big{[}a^{\dagger}_{(1,i_{y})\sigma,L}c_{(N_{x},i_{y})\sigma}+H.c.\big{]}.$
(10)

Here $i_{y}$ denotes the $y$-component of the site index $i$, and runs over
the transverse dimension $N_{y}$ of the system. For the left (right) lead, the
values of the $x$-component $i_{x}$ run backward (forward) along the $x$-axis.
The geometry of the system-reservoir system is is pictorially illustrated in
Fig. 1, with semi-infinite and spin-polarized leads coupled to the left and
right edges of the central system, supporting a spin-polarized current that
interacts with the magnetic texture. Furthermore, a schematic of the
parameters relevant to the transport setup, and to support charge and spin
currents, is reported in Fig. 2.

##  IV Spin equations of motion

We now consider the system’s nonequilibrium dynamics, which is initiated by
applying a voltage bias to the reservoirs. Since the radius of a typical
skyrmion is $\sim 100$ nm Wang2018 , a full quantum mechanical description of
the spin texture is extremely challenging. A commonly adopted strategy is then
to resort to a classical description of the spins, where $\hat{\bf
S}_{i}\rightarrow\langle\hat{\bf S}_{i}\rangle\equiv{\bf S}_{i}=S{\bf n}_{i}$
and $|{\bf n}_{i}|=1$. This limit is exact for $S\to\infty$ Lieb73 ; Fradkin13
, and considered to be a suitable approximation already when $S>1$ Heinze11 .
Taking the classical limit leads to a semi-classical approximation for the
coupled spin and electron subsystems LeadsPaper ; LaserPaper , where the
quantum electronic Hamiltonian depends parametrically on the classical
variables ${\bf n}_{i}$. In this approximation the system’s time evolution is
therefore governed by two coupled differential equations, one for the
classical spins and one for the quantum mechanical electrons.

To obtain a dynamical equation for the localized spins, we start from the
Heisenberg equations of motion for the spin operators and then take the
classical limit. The result is a Landau-Lifshitz equation of the form

$\displaystyle\derivative{\mathbf{n}_{i}}{t}$
$\displaystyle=-\mathbf{n}_{i}\crossproduct\Big{[}\sum_{\expectationvalue{j}}S(J_{ij}\mathbf{n}_{j}+\mathbf{D}_{ij}\crossproduct\mathbf{n}_{j})+(B+SKn_{iz})\hat{\bf
z}+g\langle\hat{\mathbf{s}}_{i}\rangle\Big{]},$ (11)

where the last term gives the coupling between the classical spins and the
instantaneous quantum average of the itinerant electron spins. In the
following we absorb the factors of $S$ into the couplings, i.e. we define
$SJ_{ij}\rightarrow J_{ij}$ and similarly for $\mathbf{D}_{ij}$ and $K$, and
for simplicity we assume $J_{ij}=J$. As commonly done in spintronics
simulations, we add a small Gilbert damping to the equations of motion to
stabilize the dynamics, which amount to adding an extra term
$\alpha\mathbf{n}_{i}\crossproduct(\partial\mathbf{n}_{i}/\partial t)$ in Eq.
(11). The equations of motion for the itinerant electrons is discussed in
detail in the next Section.

##  V Nonequilibrium Green’s functions theory

To describe the dynamics of the electrons, we use the theory of non-
equilibrium Green’s functions (NEGF) within the generalized Kadanoff-Baym
ansatz (GKBA). NEGF are a general and powerful approach to nonequilibrium
phenomena KadanoffBaym ; Keldysh ; balzer2012book ; Stefanucci2013 ;
Hopjan2014 , which describes the real-time dynamics of a system by exactly
including external perturbations with an arbitrary temporal dependence. Within
this theory the time-dependent expectation values of any single-particle
observable, such as currents, densities and magnetization, can be obtained
from the one-particle, two-time Green’s function

$\displaystyle
G_{ij}^{\sigma\sigma^{\prime}}(z,z^{\prime})=-i\expectationvalue{\mathcal{T}_{\gamma}\quantity{c_{i\sigma}^{H}(z)[c^{H}_{j\sigma^{\prime}}(z^{\prime})]^{\dagger}}}_{0}.$
(12)

Here the brackets $\expectationvalue{\cdot}_{0}$ denote an ensemble average
with respect to the thermal density operator $\rho=e^{-\beta H}/Z$, where
$\beta=1/k_{B}T$ is the inverse temperature, $H$ the equilibrium Hamiltonian
and $Z=\tr e^{-\beta H}$ the canonical partition function. In the low-
temperature limit $\beta\to\infty$, the ensemble average reduces to a ground
state expectation value. The operator $c_{i\sigma}^{H}(z)$ is the annihilation
operator of an itinerant electron written in the Heisenberg picture with
respect to the full time-dependent Hamiltonian $H(z)$, $z$ is a complex time
argument living on the Keldysh contour Keldysh , and all operators inside the
brackets are time-ordered on the contour by the operator
$\mathcal{T}_{\gamma}$. As before, the indices $i$ and $j$ run over all
lattice sites, while $\sigma$ and $\sigma^{\prime}$ denote spin projections.

We note that within the semi-classical scheme discussed above, the one-
particle Green’s function depends parametrically on the classical spin
variables via the spin-electron coupling $H_{s-e}$. We should therefore write
$G_{ij}^{\sigma\sigma^{\prime}}(z,z^{\prime})=G_{ij}^{\sigma\sigma^{\prime}}(z,z^{\prime};\\{{\bf
S}_{i}\\})$, but for notational simplicity we keep the dependence on ${\bf
S}_{i}$ implicit in the following. In addition, we suppress the site and spin
indexes of the Green’s function, with the implicit understanding that all
quantities are matrices in site and spin space, and only explicitly indicate
the time variables of $G(z,z^{\prime})$. The equations of motion for the one-
particle Green’s function, the so-called Kadanoff-Baym equations, can be
written as KadanoffBaym

$\displaystyle\quantity[i\derivative{z}-h(z)]G(z,z^{\prime})=\delta(z,z^{\prime})+\int_{\gamma}\differential\bar{z}\,\Sigma(z,\bar{z})G(\bar{z},z^{\prime}).$
(13)

Here $h(z)$ is the time-dependent mean-field Hamiltonian, and the self-energy
$\Sigma(z,z^{\prime})$ contain all correlation effects beyond the Hartree
approximation. A strength of the NEGF formalism is that non-interacting leads
can be exactly incorporated in the equation of motion for the Green’s function
of the central system, through the introduction of a so-called embedding self-
energy Myohanen_2008 . Denoting the correlation and embedding self-energies
respectively by $\Sigma_{c}$ and $\Sigma_{\rm emb}$, the total self-energy can
be written as $\Sigma=\Sigma_{c}+\Sigma_{\rm emb}$. We note that in the
present case the electronic Hamiltonian contains no interaction terms, and
therefore $\Sigma_{c}=0$.

In the following we consider the wide-band limit (WBL) approximation to the
embedding self-energy, obtained when the hopping amplitude inside the leads
tend to infinity while the ratio $t_{R\alpha}^{2}/t_{\alpha}$ remains fixed.
In the formal treatment, this approximation is only performed in the extended
direction of the leads (the one perpendicular to the system edge). In physical
terms, it amounts the assumption that the density of states of the leads is
constant over the band width of the central system. Performing the WBL
approximation results in the following expression for the embedding self-
energy Latini14 ; Tuovinen14 ; Ridley17 ; Tuovinen23

$\displaystyle\Sigma_{\rm emb}^{<}(t,t^{\prime})$
$\displaystyle=i\sum_{\alpha}\Gamma_{\alpha}s(t)s(t^{\prime})e^{-i\int_{t}^{t^{\prime}}V_{\alpha}(\bar{t})\differential\bar{t}}\int\frac{\differential\epsilon}{2\pi}f(\epsilon-\mu)e^{-i\epsilon(t-t^{\prime})}$
$\displaystyle\Sigma_{\rm emb}^{R}(t,t^{\prime})$
$\displaystyle=-\frac{i}{2}\delta(t-t^{\prime})s(t)s(t^{\prime})\sum_{\alpha}\Gamma_{\alpha},$

where the superscripts denote the so-called lesser ($<$) and retarded ($R$)
components of the self-energy Stefanucci2013 . Here $s(t)$ is a smooth
function used to equilibrate the central system in presence of the leads,
$V_{\alpha}(t)$ is the time-dependent bias in lead $\alpha$, and $f(\epsilon)$
is the Fermi-Dirac distribution.

###  V.1 Generalized Kadanoff-Baym ansatz

Provided that we have full knowledge of the self-energy $\Sigma$, the
Kadanoff-Baym equations are an exact reformulation of the many-body problem.
In practice however, $\Sigma$ needs to be approximated, which is commonly done
using diagrammatic many-body perturbation theory. A computational difficulty
met with the full Kadanoff-Baym formalism is that the numerical solution of
Eq. 13 scales cubically with the number of time steps, which is a consequence
of the memory integral (the right hand side of Eq. 13) accounting for the full
history of the dynamical evolution. Calculations with double-time Green’s
functions, as in Eq. (12), are highly expensive and scale unfavorably with
basis size and simulation time. Since skyrmions typically occur on large
lattice distances of $100$ nm Nagaosa , and move on typical time-scales of $1$
ps, the simulation of interacting spin-electron systems through a
straightforward use of the Kadanoff-Baym equations is prohibitive. A
significant improvement in the time-step scaling can be achieved by employing
the so-called generalized Kadanoff-Baym ansatz (GKBA) Lipavsky86 , where only
the time-diagonal of the Green’s function needs to be time evolved. This
approximation was originally derived for equilibrium systems and in the weak
scattering limit, but has been found to work well also out of equilibrium
Lipavsky86 ; Hermanns12 ; Hermanns14 ; Latini14 ; Bostrom18 ; Hopjan18 ;
KarlssonHopjan18 ; Kalvova18 ; Joost20 ; Karlsson21 ; Pavlyukh22a ;
Pavlyukh22b ; BonitzWBL23 ; Tuovinen23 . Within NEGF-GKBA and its recent
reformulation as a time-linear scheme Joost20 ; Karlsson21 , it is currently
possible to simulate the long-time dynamics of electronic systems with a basis
size on the order of $100$ orbitals Pavlyukh22b .

The GKBA for electronic degrees of freedom is achieved by the following
factorization Lipavsky86 ,

$\displaystyle
G^{<}(t,t^{\prime})=iG^{R}(t,t^{\prime})G^{<}(t^{\prime},t^{\prime})-iG^{<}(t,t)G^{A}(t,t^{\prime}),$
(14)

where the superscript $A$ denotes the advanced component of the Green’s
function. Using transformation rules introduced by Langreth Langreth1976 , we
can project Eq. 13 from the Keldysh contour to the physical time axis, and
obtain a dynamical equation for $G^{<}(t,t^{\prime})$. When taking
$t^{\prime}\to t$ and using the GKBA, the equation for $G^{<}(t,t^{\prime})$
reduces equation of motion for the one-particle density matrix
$\rho(t)=iG^{<}(t,t)$ of the form.

$\frac{\partial\rho(t)}{\partial
t}+i\commutator{h(t)}{\rho(t)}=-\Big{[}\int_{t_{0}}^{t}\differential\bar{t}\,\Sigma^{<}(t,\bar{t})G^{A}(\bar{t},t)+\int_{t_{0}}^{t}\differential\bar{t}\Sigma^{R}(t,\bar{t})\rho(\bar{t})G^{A}(\bar{t},t)\Big{]}+H.c.$
(15)

To close the equation for $\rho$ we further assume that
$G^{R}(t,t^{\prime})=-i\theta(t-t^{\prime})\mathcal{T}e^{-i\int_{t^{\prime}}^{t}d\bar{t}\,[h(\bar{t})+i\Gamma/2]}$,
with the interaction between the leads and the central region adiabatically
turned on before the bias is applied. Here
$\Gamma=\sum_{\alpha}\Gamma_{\alpha}$ accounts for the presence of the leads
Latini14 .

##  VI Approximate wide band limit (AWBL)

Introducing the GKBA improves the time-scaling of the Kadanoff-Baym equations,
but solving the equation of motion for $\rho$ still scales quadratically with
the number of time steps. Since we are interested in large systems with both
fast and slow degrees of freedom, it is useful to make further approximations
to reduce the time scaling. To this end, and inspired by other work on open
systems Tuovinen14 ; Ridley17 ; BonitzWBL23 ; Tuovinen23 , we introduce here
an approximate wide band limit (AWBL) based on an approximation to the
collision integral (the right-hand side of Eq. 15). As demonstrated below,
this prescription offers a good trade-off between computation time and
accuracy, and allows us to simulate large systems as required to describe
skyrmion textures. The proposed AWBL amounts to neglecting the spin-electron
coupling in the advanced Green’s function $G^{A}$, such that

$\displaystyle G^{A}(t,t^{\prime})$
$\displaystyle=i\theta(t^{\prime}-t){\bar{\mathcal{T}}}e^{-i\int_{t^{\prime}}^{t}d\bar{t}\,[H_{e}+H_{s-e}(\bar{t})+i\Gamma/2]}$
(16) $\displaystyle\approx
i\theta(t^{\prime}-t)e^{-i(t-t^{\prime})(H_{e}+i\Gamma/2)}.$ (17)

With this approximation the advanced Green’s function becomes a function of
the time difference, $G^{A}(t,t^{\prime})=G^{A}(t-t^{\prime})$. As a
consequence, the entire collision integral becomes a function of
$t-t^{\prime}$, only a single evaluation needs to be performed at each time
step, and the time-step scaling becomes linear (further technical details can
be found in the Supplemental Material). In general the AWBL is a quite drastic
approximation, but as demonstrated below it works rather well for the present
system. One reason for this is that the spin-electron coupling only
constitutes a higher-order correction to the embedding self-energy, since the
spins do not couple directly to the leads. A second reason is that for most of
the time-evolution, and in the dominant region of the central system, the spin
texture is constant both in time and space.

##  VII Skyrmion indicators

Before discussing the response of magnetic skyrmions to electronic currents,
we define the observables used to determine the presence of a magnetic
skyrmion in a classical spin texture. Since the spins live on the unit sphere
$S^{2}$, and the central region can be approximately compactified to the torus
$T^{2}$ (assuming a ferromagnetic ordering along the edges), a topological
charge measuring the winding number of the map ${\bf S}:T^{2}\mapsto S^{2}$
can be defined. On a lattice, it can be shown that a suitable definition of
this topological charge is given by (BERG1981, )

$\displaystyle Q$
$\displaystyle=\frac{1}{4\pi}\sum_{\\{jkl\\}\in\Delta}\Omega_{jkl},$ (18)

where the sum is over the triangulated lattice and $\Delta$ is the set that
contains the indexes for the lattice sites for every triangle in the lattice.
The solid angle $\Omega_{jkl}$ is defined by

$\displaystyle\exp(i\Omega_{jkl}/2)$
$\displaystyle=\rho_{jkl}^{-1}\big{(}1+\mathbf{S}_{j}\cdot\mathbf{S}_{k}+\mathbf{S}_{k}\cdot\mathbf{S}_{l}+\mathbf{S}_{l}\cdot\mathbf{S}_{j}+i\eta_{jkl}\mathbf{S}_{j}\cdot(\mathbf{S}_{k}\crossproduct\mathbf{S}_{l})\big{)}$
$\displaystyle\rho_{jkl}=\sqrt{2(1+\mathbf{S}_{j}\cdot\mathbf{S}_{k})(1+\mathbf{S}_{k}\cdot\mathbf{S}_{l})(1+\mathbf{S}_{l}\cdot\mathbf{S}_{j})},$
(19)

and the function $\eta_{jkl}={\rm
sgn}[(\mathbf{S}_{j}\cdot(\mathbf{S}_{k}\crossproduct\mathbf{S}_{l})]$ ensures
that the last term is positive. Geometrically the topological charge sums up
the solid angles spanned by all spin triplets $\\{jkl\\}$, which when divided
by the surface $4\pi$ of the unit sphere gives the integer winding number (for
appropriate boundary conditions). The solid angle $\Omega_{jkl}$ serves as an
indicator of the extent to which the spins twist. Specifically, within the
core of a skyrmion, $\Omega_{jkl}$ has a large magnitude, gradually
diminishing radially from the core. The extent of twisting at a specific
location can be conceptualized as the skyrmion density at that particular
site, with this density being invariably distributed across multiple sites.
Furthermore, in the presence of a skyrmion within the system we define the
center of mass as

$\displaystyle\mathbf{R}_{cm}$
$\displaystyle=\frac{1}{M}\sum_{i_{x}=2}^{N_{x}-1}\sum_{i_{y}=2}^{N_{y}-1}\rho_{i_{x},i_{y}}\begin{pmatrix}i_{x}\\\
i_{y}\end{pmatrix},$ (20)

The expressions for $M$ and $\rho_{i_{x},i_{y}}$ are respectively given by

$\displaystyle M$
$\displaystyle=\sum_{i_{x}=2}^{N_{x}-1}\sum_{i_{y}=2}^{N_{y}-1}\rho_{i_{x},i_{y}},$
(21) $\displaystyle\rho_{i_{x},i_{y}}$
$\displaystyle=\bigg{|}\Omega_{(i_{x},i_{y}),(i_{x}+1,i_{y}),(i_{x}+1,i_{y}+1)}$
(22)
$\displaystyle+\Omega_{(i_{x},i_{y}),(i_{x},i_{y}+1),(i_{x}+1,i_{y}+1)}\bigg{|}.$
(23)

In Eqs. (20, 21), the index $i_{x}$ ($i_{y}$) ranges between 2 and $N_{x}-1$
(2 and $N_{y}-1$) rather than, as intuitively expected, between 1 and $N_{x}$
(1 and $N_{y}$). Observations of spin oscillations at the periphery of the
central region motivate the exclusion of these edges in the determination of
the skyrmion center of mass. The rationale behind these oscillations and the
specific way we perform the site exclusion are further detailed in Section
VIII.1.

##  VIII Current driven skyrmion motion

We now consider the motion of skyrmions generated by a current density in the
itinerant electron system, in turn driven by an external bias. However, first
we need to prepare the system in a state featuring a skyrmion, a state that
may not necessarily correspond to the ground state of the semi-classical
system. This entails the selection of a suitable initial spin configuration,
subsequently coupled in a self-consistent manner to the electron system. The
interaction with the leads is slowly turned on in the time interval $t<\tau$
using the contact function $s(t)=\sin^{2}(\pi t/2\tau)$. After time $t=\tau$ a
bias is applied, here considered to be of the form

$\displaystyle v_{\alpha}(t)=\theta(t-\tau)\begin{cases}\phantom{-}V&\text{if
}\alpha=L\text{ and }\sigma=\,\uparrow\\\ -V&\text{if }\alpha=R\text{ and
}\sigma=\,\uparrow\\\ \phantom{-}0&\text{otherwise}\end{cases}.$ (24)

This bias generates a spin current through the central system. In the
simulations below we fix the the energy unit by setting the hopping
$t_{ij}=t_{s}=-1$, and take $\hbar=1$.

###  VIII.1 Benchmarks for the approximate wide band limit, and the role of
damping

As a primer for our discussion of disorder effects on skyrmion dynamics, we
briefly investigate the performance of the AWBL in the absence of disorder.
For this purpose, we compare the dynamics obtained within the WBL and AWBL of
a small central region denoted by $C_{16}$, consisting of a $4\times 4$
cluster as illustrated in Fig. 3a.

We start from the ground state of the isolated region $C_{16}$, and gradually
connect the system to the leads in the time interval $t\in[0,800]$ measured in
units of $\hbar/t_{s}$. This prepares the stationary initial state of the
fully connected system. At time $t=800$, a symmetric and spin-polarized bias
of the form shown in Eq. 24 is then applied. The subsequent part of the
simulation, i.e. for times $t>800$, is the one relevant for the skyrmion
dynamics to be discussed later. To compare the performance of the WBL and
AWBL, we choose as indicators the spin current and the spin-up density at a
given site. Representative results are shown in Fig. 3b-d, and indicate that
the agreement between the WBL and AWBL varies noticeably with the strength of
the spin-electron coupling $g$ and the lead-device connection $\gamma$.

Figure 3: Benchmarks for the approximate wide band limit. $(a)$ Initial
relaxed configuration of a $4\times 4$ spin texture interacting with a
corresponding layer of electrons. $(b-d)$ Spin current $I_{R,\uparrow}$ and
electronic spin density $\rho_{11\uparrow}$ of the lattice site encircled in
green in $(a)$, for different values of $\gamma$ and $g$, as obtained within
the wide band limit (WBL) and approximate wide band limit (AWBL). The spin
texture is approximately the same for all cases and throughout the time
evolution, due to a large value of the spin damping $\alpha$. The parameters
used are $B=0.05$, $K=0.6$, $\alpha=3$, $\alpha_{R}=0$, $J=0.5$, $D=0.4$,
$\tau=800$, $\mu=0$ and $V=3$. The color on the spins indicates their
$z$-component, while the arrows show their direction in the $xy$-plane. The
values of $g$ and $\gamma$ are given in each panel.

It is useful at this point to make a couple of remarks about the damping term
in the spin dynamics. The results in Fig. 3b-d are obtained with a large
Gilbert damping of $\alpha=3.0$, which results in a spin texture that remains
essentially stationary during the time evolution. Therefore, these comparisons
mainly concern the regime of electron dynamics in presence of a static
magnetic background. Using the initial spin of Fig. 3a) in $G^{A}$ instead of
neglecting $H_{s-e}$ would give a very good agreement between the AWBL and the
WBL. Strictly speaking, one can modify the AWBL to include a fixed spin
texture in $G^{A}$ instead of neglecting it completely and still retain time-
linear scaling. However, for the time-dependent disorder simulations to be
discussed later, the spin texture changes significantly in time due to the
skyrmion motion. We have verified numerically that, replacing in $G^{A}$ the
exact contribution from the time evolving texture with a static one (chosen at
any time during the time evolution) is not better than using a $G^{A}$ with
zero spin-electron coupling. This is why no spin configuration is included in
$G^{A}$ in Eq.(16), or later on in the dynamics in presence of disorder.

Figure 4: Skyrmion motion in response to electronic spin currents. (a) Center
of mass motion of the skyrmion and (b) initial spin configuration at $t=200$.
The highlighted region corresponds to the area shown in Figs. 5 to 7. The
parameters used are $g=1$, $B=0.05$, $K=0.6$, $\alpha=0.1$, $\alpha_{R}=0$,
$J=0.5$, $D=0.4$ (corresponding to Bloch-type DMI), $\gamma=0.2$, $\tau=1500$,
$\mu=0$ and $V=2$. The color of the spins indicates their $z$-component, while
the arrows display their direction in the $xy$-plane.

The purpose of using a large $\alpha$ in Fig. 3 is to attain a stationary
initial state before the bias is applied. In Fig. 3, such large damping is
maintained throughout the whole time evolution (i.e. during the lead
attachment and afterwards, when the bias is applied), for consistency and to
focus on the electronic behavior. However, for the later results in the paper
(Sect. VIII.2 and afterwards, where we consider skyrmion dynamics in central
regions of $20\times 40$ sites), the damping is set small ($\alpha=0.1$), in
order not to influence the system’s intrinsic dynamics, but still facilitate
the adiabatic preparation of the initial state. In this respect, we have found
that, even when evolving the system with no bias, and using very small or no
damping after fully connecting the leads, there is a very slow build-up of
oscillations in the current and the spin densities. Such oscillations
primarily occur in the sites at the edges of the central region (and
especially at the sites in contact with the leads). This holds both for the
small cluster in Fig. 3a and for the larger systems studied later, and
indicates that when the bias is applied the system has has not yet reached in
full the ground state in the presence of the leads, due to the very complex
energy and fluctuating landscape provided by the interaction of the electrons
with the classical spins. Though, in the spirit of having a microscopic
current inducing the skyrmion motion, it is still meaningful to apply a bias
to this configuration, and interpret the oscillations as physical in
character. At the same time, these oscillations affect in a rather artificial
way the estimate of the skyrmion center of mass $\mathbf{R}_{cm}$.
Accordingly, in the results in the next sections, $\mathbf{R}_{cm}$ is
calculated using Eq. 20, i.e without including the contribution of the
peripheral sites of the $20\times 40$ region.

Coming back to the results in Fig. 3b-d, an interesting feature is that the
agreement between the AWBL and WBL spin currents, for small $\gamma$, is
generally better than agreement between the spin densities. While we do not
have a clear explanation for this, we note that this trend is confirmed in
additional simulations (not shown), where $\gamma$ is varied at constant $g$
and viceversa. To summarize, our results suggest that the AWBL produces fairly
accurate electronic and spin currents, overall in good agreement with the full
WBL. It thereby constitutes a microscopic and semi-quantitative method to
investigate skyrmion dynamics at low computational cost, appropriate for both
transient and steady state regimes.

Figure 5: Skyrmion motion in presence of a small impurity on the skyrmion
path. (a) Center of mass movement of the skyrmion and (b) some snapshots of
the $z$-component of the spins in the highlighted area of Fig. 4. The
parameters used are $g=1$, $B=0.05$, $K=0.6$, $\alpha=0.1$, $\alpha_{R}=0$,
$J=0.5$, $D=0.4$ (corresponding to Bloch-type DMI), $\gamma=0.2$, $\tau=1500$,
$\mu=0$ and $V=2$. The color of the spins indicates their $z$-component.

###  VIII.2 Current induced motion

Having shown that the AWBL is in good quantitative agreement with the WBL for
small systems, we now employ it to describe skyrmion motion in systems beyond
the scope of the WBL. These results provide a useful benchmark for the
discussion in later sections, the effect of magnetic impurities on skyrmion
motion is analyzed. To obtain the reference results, we consider a central
region with $40\times 20$ sites, and investigate the dynamics in response to
an applied bias. These results are henceforth referred to as ”base case”. In
the base case as well as in the impurity studies below, a potential of
strength $V=2$ is applied in the leads, and the coupling between the system
and leads is of strength $\gamma=0.2$.

The system starts off with a relaxed skyrmion texture in the center (see Fig.
4b), and in response to the bias starts to drift towards the right with a
slight additional downwards movement (see Fig. 4a). As seen from the figure,
the velocity in the $x$-direction is not constant but rather ramps up, then
slows down, and finally starts to ramp up again. While the dominant rightwards
motion is due to the spin current, the downwards force is due to the so-called
skyrmion Hall effect Nagaosa . Additionally, there are some small-scale
oscillations, which are especially visible in the $y$-direction. These
oscillations could be due to the discrete nature of the underlying lattice,
since the skyrmion has a preferred equilibrium. When it moves between lattice
sites, the skyrmion form gets distorted, but returns to its original
configuration when it arrives at a new site. We finally note that it takes
some time for the skyrmion to pick up speed as the current is going through
the system, which is a general trend in our simulations. This is indicative of
the presence of a finite skyrmion mass, which has been discussed in previous
studies PhysRevB.90.174434 ; PhysRevLett.109.217201  In particular, it was
recently proposed that the magnitude of the skyrmion mass is affected by the
interactions between localized moments and free electrons Nikolic1 ; Nikolic2
, consistent with the present findings.

##  IX Disorder effects

Knowing the trajectory of the skyrmion in the base case, we can now
investigate the effects of putting impurities in its path. In the panels of
Figs. 5 to 7 showing the spin configuration, the positions of the impurities
are marked in green. Here we restrict to a simple form of impurity,
corresponding to a reduction in the symmetric exchange interaction by $J\to
J/2$ between the impurity site and its nearest neighbors. Therefore, it is
expected that the skyrmion will be hindered by a smaller interaction. This
change is effected at $t=1400$ in order not to disturb the system before the
bias is applied at $t=1500$.

Figure 6: Skyrmion motion in presence of a small impurity near the skyrmion
path. (a) Center of mass movement of the skyrmion and (b) some snapshots of
the $z$-component of the spins in the highlighted area of Fig. 4. The
parameters used are $g=1$, $B=0.05$, $K=0.6$, $\alpha=0.1$, $\alpha_{R}=0$,
$J=0.5$, $D=0.4$ (corresponding to Bloch-type DMI), $\gamma=0.2$, $\tau=1500$,
$\mu=0$ and $V=2$. The color of the spins indicates their $z$-component.

###  IX.1 Case 1: A small impurity on the skyrmion path

As a typical example of a small magnetic impurity, we consider an impurity
cluster consisting of a square of four lattice sites (see Fig. 5). One
possible realization of this form of disorder is the presence of adsorbate
atoms in the center of these four lattice sites. The impurity is placed right
in the path of the skyrmion, slightly to the right and below the skyrmion’s
initial position.

The center of mass movement of the skyrmion in response to a spin current is
shown in Fig. 5. In the $x$-direction, it seems like the skyrmion does not
feel the effect of the impurity until it is on top of it, since the movement
of the center of mass is quite similar to the base case before $t\sim 1700$.
However, a careful comparison reveals that it moves slightly slower. In
contrast, the skyrmion starts to pick up speed in the $y$-direction as it
approaches the impurity, indicating that the velocity in the $x$-direction is
transferred to a downwards motion. After the skyrmion has entered the
impurity, it gets pinned and struggles to escape. More specifically, between
$t=1700$ and $t=2000$ the skyrmion is stuck to the impurity, and can only move
around as long as some part of it is still on the impurity. A few times during
the center of mass dynamics (see Fig. 5), there are bumps, which indicate that
the skyrmion gathered velocity to try to escape but was pulled back again by
the impurity. Finally, at $t\sim 2000$, the skyrmion manages to escape and
move past the impurity. The dynamics presented in Fig. 5 clearly show the
advantage of a microscopic description, where the motion of the skyrmion can
be tracked atom by atom. During the dynamics the shape of the skyrmion is
changing and the direction of motion changes several times, effects that are
not possible to capture with a Thiele-type description. Furthermore, it is
quite reasonable to expect that the observed dynamics are a direct result of
the motion of the itinerant electrons, especially at the skyrmion boundary,
where the local spin changes must take into account the redistribution of the
electronic spin density during the time evolution.

###  IX.2 Case 2: A small impurity near the skyrmion path

Next, we shift the impurity position two lattice sites upwards compared to the
previous section (Fig. 6), such that the impurity is now slightly above the
skyrmion’s path. At the beginning of the dynamics, the center of mass movement
seems very similar to the previous case, with the velocity along the
$x$-direction slightly reduced compared to the base case, and a significant
downwards movement. However, in the present case the movement does not force
the skyrmion on top of the impurity, and the impurity pinning is avoided. In
the center of the mass movement of Fig. 6, as well as in the snapshots at
$t=1700$ and $t=1800$ of Fig. 6, the skyrmion is seen to move beneath the
impurity. It is likely that, if the system was extended further in the
$y$-direction, the skyrmion would continue its motion downwards. In the
present case however, the edge of the system prevents this, since our
simulations suggest that the repulsive force induced by the edges makes the
skyrmion to recoil upwards in the $y$-direction. This illustrates, on the
microscopic level, the importance of considering the track width when
designing efficient racetrack memories. If the width of the track is too
small, the skyrmion might not be able to get around impurities, while if the
width is too large skyrmions can be deflected by impurities and substantially
deviate from their path.

Although the skyrmion’s core, shown in blue in Fig. 4, is comprised of only
$4\times 4$ spins, a rotation of the $xy$-components of the spins away from
the ferromagnetic background can be seen in an area extending up to $12\times
12$ sites. This indicates that the effective skyrmion size is much larger than
the central $4\times 4$ region, which in turn implies that the skyrmion feels
the presence of the edge much earlier that would be expected from only
considering its core.

###  IX.3 Case 3: A larger impurity configuration

Figure 7: Skyrmion motion in presence of an extended impurity. (a) Center of
mass movement of the skyrmion and (b) some snapshots of the $z$-component of
the spins in the highlighted area of Fig. 4. The parameters used are $g=1$,
$B=0.05$, $K=0.6$, $\alpha=0.1$, $\alpha_{R}=0$, $J=0.5$, $D=0.4$
(corresponding to Bloch-type DMI), $\gamma=0.2$, $\tau=1500$, $\mu=0$ and
$V=2$. The color of the spins indicates their $z$-component.

As an example of dynamics in presence of a more extended magnetic impurity,
Fig. 7 shows the center of mass motion of the skyrmion in a system containing
a columnar impurity that covers almost the entire lattice in the
$y$-direction. Again, the impurity is only two lattice sites thick, but this
time with only four sites at the top and the bottom left unperturbed. The
characteristic downwards bump in the center of mass motion along the
y-direction is also seen here, but it is twice as big as in the previous
cases. As seen in Fig. 7, the velocity in the $y$-direction is large before
the skyrmion is suddenly halted. This detail, together with the observation
that at $t=1900$ the skyrmion is slightly smaller, suggests that a small
skyrmion could be annihilated even by a small impurity. After $t=1900$, the
skyrmion recoils upwards but the collision with the top does not seem as
severe as in the case of Fig. 5. Likely, this is because the velocity is lower
and the skyrmion is already stuck.

The observed motion along the vertical direction is most likely aided by the
lower wall. However, the impurity likely also plays an important role, since
we have seen that it is hard for the skyrmion to escape once it moves on top
of the impurity. In addition to the vertical motion, Fig. 7 also shows that
the skyrmion keeps wiggling along the $x$-direction while moving. In
particular, the snapshot at $t=1900$ strongly indicates that this movement is
due to the skyrmion expanding to the right. Subsequently, for times after
$t\sim 2100$, the skyrmion shrinks and expands in a quasi steady state.

###  IX.4 Summary of the results

To summarize this section, the results of Figs. 5 to 7 show that by placing
impurities in the selected places along a skyrmion’s path, it is possible to
engineer the skyrmion trajectory and to steer skyrmions in the desired
direction in logic devices PhysRevApplied.12.064053 . At the same time it is
also clear that, in order to accurately engineer the desired skyrmion dynamics
in the presence of impurities and provide some conceptual guidance for the
experimental realization of skyrmion architectures, a microscopic and explicit
equal-footing account of the behavior of itinerant electrons and localized
spins is of high relevance.

##  X Discussion and outlook

In this work we have presented a recently introduced method to deal with
magnetic skyrmions, where a skyrmion magnetic texture made of classical spins
is embedded in a background of quantum electrons. We considered the skyrmion
dynamics in a quantum transport geometry, where a central region of coupled
electrons and spins is connected to electron reservoirs. Our treatment is
based on a two-component description where the electrons are treated via
nonequilibrium Green’s functions (NEGF), coupled to classical spins governed
by a Landau-Lifshitz-Gilbert (LLG) equation.

To describe the skyrmion motion in large systems, we introduced an approximate
approach to include the effects electronic reservoirs coupled to the centeral
system. This so-called approximate wide-band limit (AWBL) shows a satisfactory
level of agreement with the full wide-band limit (WBL) approach. Similar to
other schemes in the literature it is a time-linear scaling method, and in
addition has a very advantageous scaling prefactor of great convenience to
deal with systems of large size.

This approach was used to investigate the effects of magnetic impurities on
the current-induced motion of magnetic skyrmion. Our results show that it is
possible to characterize at the atomic level processes such as skyrmion
scattering, recoil, drift, and trapping in disordered samples. This is a
subject of high relevance, since understanding how tailor the disorder level
of a given sample is crucial in order to control skyrmion dynamics.

A natural extension of this work is to take into account the additional
correlation effects arising from electron-electron interaction. This would
bring the simulation closer to realistic systems, where correlations between
the electrons could e.g. lead to destabilization of a skyrmion, or introduce
effects a competing with disorder. In a different direction, the results
obtained here can be used to benchmark other methods that deal with electron
and localized spins. For example, one can imagine to validate methods that,
still in a framework, go beyond the standard LLG equation Nikolic1 ; Nikolic2
.

Finally, we note that the methods presented here open the door for simulations
of the intertwined spin-electron dynamics of small but realistic nanoscale
systems. This provides a new approach to investigate the ultrafast dynamics of
magnetic systems with general noncollinear orders, as initiated by ultrashort
laser or current pulses LaserPaper ; Truc2023 . Apart from its importance for
skyrmionic systems, the presented method will therefore be of large relevance
to describe the intertwined dynamics of spins, electrons and lattice
vibrations in spintronics devices and for material engineering, thereby
helping to facilitate a microscopic understanding of driven magnetic
materials.

## Conflict of Interest Statement

The authors declare no competing financial or non-financial interests.

## Author Contributions

E.Ö. performed all calculations and interpretation of results under the
supervision of E.V.B. and C.V. E.Ö also contributed to the extension of pre-
existing code. The project was conceived by E.V.B. and C.V. The overall
supervision of the project was by C.V. All authors collaborated in writing the
paper.

## Funding

EVB acknowledges funding from the European Union’s Horizon Europe research and
innovation programme under the Marie Skłodowska-Curie grant agreement No
101106809. CV acknowledges funding from the Swedish Research Council (grant No
VR 2022 04486)

## Data Availability Statement

The data supporting the findings of the present manuscript, as well as the
files used to generate the figures, are available from the authors upon
reasonable request.

## References

  * (1) Cucchiella, F., D’Adamo, I., Lenny Koh, S. & Rosa, P.  Recycling of weees: An economic assessment of present and future e-waste streams.  _Renewable and Sustainable Energy Reviews_ 51, 263–272 (2015).  URL https://www.sciencedirect.com/science/article/pii/S1364032115005808.
  * (2) Irimia-Vladu, M.  “green” electronics: biodegradable and biocompatible materials and devices for sustainable future.  _Chemical Society Reviews_ 43, 588–610 (2014).
  * (3) Hirohata, A. _et al._ Review on spintronics: Principles and device applications.  _Journal of Magnetism and Magnetic Materials_ 509 (2020).
  * (4) Coronado, E.  Molecular magnetism: from chemical design to spin control in molecules, materials and devices.  _Nature Reviews Materials_ 5, 87–104 (2020).
  * (5) Parkin, S. & Yang, S.-H.  Memory on the racetrack.  _Nature Nanotechnology_ 10, 195–198 (2015).  URL https://doi.org/10.1038/nnano.2015.41.
  * (6) Göbel, B., Mertig, I. & Tretiakov, O. A.  Beyond skyrmions: Review and perspectives of alternative magnetic quasiparticles.  _Physics Reports_ 895, 1–28 (2021).  URL https://www.sciencedirect.com/science/article/pii/S0370157320303525.
  * (7) Back, C. _et al._ The 2020 skyrmionics roadmap.  _Journal of Physics D: Applied Physics_ 53, 363001 (2020).  URL https://dx.doi.org/10.1088/1361-6463/ab8418.
  * (8) Sampaio, J., Cros, V., Rohart, S., Thiaville, A. & Fert, A.  Nucleation, stability and current-induced motion of isolated magnetic skyrmions in nanostructures.  _Nature Nanotechnology_ 8, 839–844 (2013).  URL https://doi.org/10.1038/nnano.2013.210.
  * (9) Nagaosa, N. & Tokura, Y.  Topological properties and dynamics of magnetic skyrmions.  _Nature Nanotechnology_ 8, 899–911 (2013).  URL https://www.nature.com/articles/nnano.2013.243.
  * (10) Romming, N., Kubetzka, A., Hanneken, C., von Bergmann, K. & Wiesendanger, R.  Field-dependent size and shape of single magnetic skyrmions.  _Phys. Rev. Lett._ 114, 177203 (2015).  URL https://link.aps.org/doi/10.1103/PhysRevLett.114.177203.
  * (11) Yang, G., Stano, P., Klinovaja, J. & Loss, D.  Majorana bound states in magnetic skyrmions.  _Phys. Rev. B_ 93, 224505 (2016).  URL https://link.aps.org/doi/10.1103/PhysRevB.93.224505.
  * (12) Chauwin, M. _et al._ Skyrmion logic system for large-scale reversible computation.  _Phys. Rev. Appl._ 12, 064053 (2019).  URL https://link.aps.org/doi/10.1103/PhysRevApplied.12.064053.
  * (13) Rex, S., Gornyi, I. V. & Mirlin, A. D.  Majorana bound states in magnetic skyrmions imposed onto a superconductor.  _Phys. Rev. B_ 100, 064504 (2019).  URL https://link.aps.org/doi/10.1103/PhysRevB.100.064504.
  * (14) Psaroudaki, C. & Panagopoulos, C.  Skyrmion qubits: A new class of quantum logic elements based on nanoscale magnetization.  _Phys. Rev. Lett._ 127, 067201 (2021).  URL https://link.aps.org/doi/10.1103/PhysRevLett.127.067201.
  * (15) Viñas Boström, E. & Verdozzi, C.  Steering magnetic skyrmions with currents: A nonequilibrium green’s functions approach.  _physica status solidi (b)_ 256, 1800590 (2019).  URL https://onlinelibrary.wiley.com/doi/abs/10.1002/pssb.201800590.
  * (16) Viñas Boström, E., Rubio, A. & Verdozzi, C.  Microscopic theory of light-induced ultrafast skyrmion excitation in transition metal films.  _npj Comput Mater_ 8 (2022).  URL https://doi.org/10.1038/s41524-022-00735-5.
  * (17) Lipavský, P., Špička, V. & Velický, B.  Generalized kadanoff-baym ansatz for deriving quantum transport equations.  _Phys. Rev. B_ 34, 6933–6942 (1986).  URL https://link.aps.org/doi/10.1103/PhysRevB.34.6933.
  * (18) Hermanns, S., Balzer, K. & Bonitz, M.  The non-equilibrium green function approach to inhomogeneous quantum many-body systems using the generalized kadanoff–baym ansatz.  _Physica Scripta_ 2012, 014036 (2012).  URL https://dx.doi.org/10.1088/0031-8949/2012/T151/014036.
  * (19) Hermanns, S., Schlünzen, N. & Bonitz, M.  Hubbard nanoclusters far from equilibrium.  _Phys. Rev. B_ 90, 125111 (2014).  URL https://link.aps.org/doi/10.1103/PhysRevB.90.125111.
  * (20) Latini, S., Perfetto, E., Uimonen, A.-M., van Leeuwen, R. & Stefanucci, G.  Charge dynamics in molecular junctions: Nonequilibrium green’s function approach made fast.  _Phys. Rev. B_ 89, 075306 (2014).  URL https://link.aps.org/doi/10.1103/PhysRevB.89.075306.
  * (21) Boström, E. V., Mikkelsen, A., Verdozzi, C., Perfetto, E. & Stefanucci, G.  Charge separation in donor–c60 complexes with real-time green functions: The importance of nonlocal correlations.  _Nano Letters_ 18, 785–792 (2018).  URL https://doi.org/10.1021/acs.nanolett.7b03995.
  * (22) Hopjan, M., Stefanucci, G., Perfetto, E. & Verdozzi, C.  Molecular junctions and molecular motors: Including coulomb repulsion in electronic friction using nonequilibrium green’s functions.  _Phys. Rev. B_ 98, 041405 (2018).  URL https://link.aps.org/doi/10.1103/PhysRevB.98.041405.
  * (23) Karlsson, D., Hopjan, M. & Verdozzi, C.  Disorder and interactions in systems out of equilibrium: The exact independent-particle picture from density functional theory.  _Phys. Rev. B_ 97, 125151 (2018).  URL https://link.aps.org/doi/10.1103/PhysRevB.97.125151.
  * (24) Kalvová, A., Velický, B. & Špička, V.  Generalized master equation for a molecular bridge improved by vertex correction to the generalized kadanoff-baym ansatz.  _Europhysics Letters_ 121, 67002 (2018).  URL https://dx.doi.org/10.1209/0295-5075/121/67002.
  * (25) Schlünzen, N., Joost, J.-P. & Bonitz, M.  Achieving the scaling limit for nonequilibrium green functions simulations.  _Phys. Rev. Lett._ 124, 076601 (2020).  URL https://link.aps.org/doi/10.1103/PhysRevLett.124.076601.
  * (26) Karlsson, D., van Leeuwen, R., Pavlyukh, Y., Perfetto, E. & Stefanucci, G.  Fast green’s function method for ultrafast electron-boson dynamics.  _Phys. Rev. Lett._ 127, 036402 (2021).  URL https://link.aps.org/doi/10.1103/PhysRevLett.127.036402.
  * (27) Pavlyukh, Y., Perfetto, E., Karlsson, D., van Leeuwen, R. & Stefanucci, G.  Time-linear scaling nonequilibrium green’s function methods for real-time simulations of interacting electrons and bosons. i. formalism.  _Phys. Rev. B_ 105, 125134 (2022).  URL https://link.aps.org/doi/10.1103/PhysRevB.105.125134.
  * (28) Pavlyukh, Y., Perfetto, E., Karlsson, D., van Leeuwen, R. & Stefanucci, G.  Time-linear scaling nonequilibrium green’s function method for real-time simulations of interacting electrons and bosons. ii. dynamics of polarons and doublons.  _Phys. Rev. B_ 105, 125135 (2022).  URL https://link.aps.org/doi/10.1103/PhysRevB.105.125135.
  * (29) Balzer, K., Schlünzen, N., Ohldag, H., Joost, J.-P. & Bonitz, M.  Accelerating nonequilibrium green function simulations with embedding self-energies.  _Phys. Rev. B_ 107, 155141 (2023).  URL https://link.aps.org/doi/10.1103/PhysRevB.107.155141.
  * (30) Tuovinen, R., Pavlyukh, Y., Perfetto, E. & Stefanucci, G.  Time-linear quantum transport simulations with correlated nonequilibrium green’s functions.  _Phys. Rev. Lett._ 130, 246301 (2023).  URL https://link.aps.org/doi/10.1103/PhysRevLett.130.246301.
  * (31) Iwasaki, J., Mochizuki, M. & Nagaosa, N.  Universal current-velocity relation of skyrmion motion in chiral magnets.  _Nature Communications_ 4 (2013).
  * (32) Stier, M., Strobel, R., Krause, S., Häusler, W. & Thorwart, M.  Role of impurity clusters for the current-driven motion of magnetic skyrmions.  _Phys. Rev. B_ 103, 054420 (2021).  URL https://link.aps.org/doi/10.1103/PhysRevB.103.054420.
  * (33) Haller, A., Groenendijk, S., Habibi, A., Michels, A. & Schmidt, T. L.  Quantum skyrmion lattices in heisenberg ferromagnets.  _Phys. Rev. Res._ 4, 043113 (2022).  URL https://link.aps.org/doi/10.1103/PhysRevResearch.4.043113.
  * (34) Everschor-Sitte, K., Masell, J., Reeve, R. M. & Kläui, M.  Perspective: Magnetic skyrmions—overview of recent progress in an active research field.  _Journal of Applied Physics_ 124, 240901 (2018).  URL https://doi.org/10.1063/1.5048972.
  * (35) Heinze, S. _et al._ Spontaneous atomic-scale magnetic skyrmion lattice in two dimensions.  _Nat. Phys._ 7, 713–718 (2011).  URL https://doi.org/10.1038/nphys2045.
  * (36) Thiele, A. A.  Steady-state motion of magnetic domains.  _Phys. Rev. Lett._ 30, 230–233 (1973).  URL https://link.aps.org/doi/10.1103/PhysRevLett.30.230.
  * (37) Reyes-Osorio, F. & Nikolic, B. K.  Anisotropic skyrmion mass induced by surrounding conduction electrons: A Schwinger-Keldysh field theory approach.  _arXiv.2302.04220_ (2023).
  * (38) Reyes-Osorio, F. & Nikolic, B. K.  Gilbert damping in metallic ferromagnets from Schwinger-Keldysh field theory: Nonlocality, nonuniformity, and anisotropy in the presence of spin-orbit coupling.  _arXiv.2302.04220_ (2023).
  * (39) Weißenhofer, M., Rózsa, L. & Nowak, U.  Skyrmion dynamics at finite temperatures: Beyond Thiele’s equation.  _Phys. Rev. Lett._ 127, 047203 (2021).  URL https://link.aps.org/doi/10.1103/PhysRevLett.127.047203.
  * (40) Huang, L., Burnell, G. & Marrows, C. H.  Transient retrograde motion of spin wave driven skyrmions in magnetic nanotracks.  _Phys. Rev. B_ 107, 224418 (2023).  URL https://link.aps.org/doi/10.1103/PhysRevB.107.224418.
  * (41) Baláž, P. _et al._ Theory of superdiffusive spin transport in noncollinear magnetic multilayers.  _Phys. Rev. B_ 107, 174418 (2023).  URL https://link.aps.org/doi/10.1103/PhysRevB.107.174418.
  * (42) Gupta, R. _et al._ Element-resolved evidence of superdiffusive spin current arising from ultrafast demagnetization process.  _Phys. Rev. B_ 108, 064427 (2023).  URL https://link.aps.org/doi/10.1103/PhysRevB.108.064427.
  * (43) Ghosh, S., Freimuth, F., Gomonay, O., Blügel, S. & Mokrousov, Y.  Driving spin chirality by electron dynamics in laser-excited antiferromagnets.  _Communications Physics_ 5, 69 (2022).  URL https://doi.org/10.1038/s42005-022-00840-3.
  * (44) Zarzuela, R. & Sinova, J.  Effective spin-charge transport theory and spin-transfer physics in frustrated magnets within the slave-boson approach.  _Phys. Rev. B_ 105, 024423 (2022).  URL https://link.aps.org/doi/10.1103/PhysRevB.105.024423.
  * (45) Ghosh, S., Blügel, S. & Mokrousov, Y.  Ultrafast optical generation of antiferromagnetic meron-antimeron pairs with conservation of topological charge.  _Phys. Rev. Res._ 5, L022007 (2023).  URL https://link.aps.org/doi/10.1103/PhysRevResearch.5.L022007.
  * (46) Dzyaloshinsky, I.  A thermodynamic theory of ”weak” ferromagnetism of antiferromagnetics.  _J. Phys. Chem. Solids_ 4, 241 (1958).  URL http://www.sciencedirect.com/science/article/pii/0022369758900763.
  * (47) Moriya, T.  Anisotropic superexchange interaction and weak ferromagnetism.  _Phys. Rev._ 120, 91 (1960).  URL https://link.aps.org/doi/10.1103/PhysRev.120.91.
  * (48) Wang, X. S., Yuan, H. Y. & Wang, X. R.  A theory on skyrmion size.  _Communications Physics_ 1, 31 (2018).  URL https://doi.org/10.1038/s42005-018-0029-0.
  * (49) Lieb, E. H.  The classical limit of quantum spin systems.  _Commun. Math. Phys._ 31, 327–340 (1973).  URL https://doi.org/10.1007/bf01646493.
  * (50) Fradkin, E.  _Field Theories of Condensed Matter Physics_ (Cambridge University Press, 2013).
  * (51) Kadanoff, L. P. & Baym, G.  _Quantum statistical mechanics : Green’s function methods in equilibrium and nonequilibrium problems_ (W.A. Benjamin, 1962).
  * (52) Keldysh, L. V.  Diagram technique for nonequilibrium processes.  _Soviet Phys. JETP_ 20, 1018–26 (1965).
  * (53) Balzer, K. & Bonitz, M.  _Nonequilibrium Green’s Functions Approach to Inhomogeneous Systems_.  Lecture Notes in Physics (Springer Berlin Heidelberg, 2012).  URL https://books.google.se/books?id=2-XXMgEACAAJ.
  * (54) Stefanucci, G. & van Leeuwen, R.  _Nonequilibrium Many-Body Theory of Quantum Systems: A Modern Introduction_ (Cambridge University Press, 2013).
  * (55) Hopjan, M. & Verdozzi, C.  _Probing Strongly Correlated Materials in Non-equilibrium: Basic Concepts and Possible Future Trends in First Principle Approaches_ (Springer Berlin Heidelberg, Berlin, Heidelberg, 2014).  URL https://doi.org/10.1007/128_2014_528.
  * (56) Myöhänen, P., Stan, A., Stefanucci, G. & van Leeuwen, R.  A many-body approach to quantum transport dynamics: Initial correlations and memory effects.  _Europhysics Letters_ 84, 67001 (2008).  URL https://dx.doi.org/10.1209/0295-5075/84/67001.
  * (57) Tuovinen, R., Perfetto, E., Stefanucci, G. & van Leeuwen, R.  Time-dependent landauer-büttiker formula: Application to transient dynamics in graphene nanoribbons.  _Phys. Rev. B_ 89, 085131 (2014).  URL https://link.aps.org/doi/10.1103/PhysRevB.89.085131.
  * (58) Ridley, M., MacKinnon, A. & Kantorovich, L.  Partition-free theory of time-dependent current correlations in nanojunctions in response to an arbitrary time-dependent bias.  _Phys. Rev. B_ 95, 165440 (2017).  URL https://link.aps.org/doi/10.1103/PhysRevB.95.165440.
  * (59) Langreth, D. C.  Linear and nonlinear response theory with applications.  In Devreese, J. T. & van Doren, V. E. (eds.) _Linear and Nonlinear Electron Transport in Solids_ , 3–32 (Springer US, Boston, MA, 1976).  URL https://doi.org/10.1007/978-1-4757-0875-2_1.
  * (60) Berg, B. & Lüscher, M.  Definition and statistical distributions of a topological number in the lattice o(3) $\sigma$-model.  _Nuclear Phys. B_ 190, 412–424 (1981).  URL https://doi.org/10.1016/0550-3213(81)90568-X.
  * (61) Schütte, C., Iwasaki, J., Rosch, A. & Nagaosa, N.  Inertia, diffusion, and dynamics of a driven skyrmion.  _Phys. Rev. B_ 90, 174434 (2014).  URL https://link.aps.org/doi/10.1103/PhysRevB.90.174434.
  * (62) Makhfudz, I., Krüger, B. & Tchernyshyov, O.  Inertia and chiral edge modes of a skyrmion magnetic bubble.  _Phys. Rev. Lett._ 109, 217201 (2012).  URL https://link.aps.org/doi/10.1103/PhysRevLett.109.217201.
  * (63) Chauwin, M. _et al._ Skyrmion logic system for large-scale reversible computation.  _Phys. Rev. Appl._ 12, 064053 (2019).  URL https://link.aps.org/doi/10.1103/PhysRevApplied.12.064053.
  * (64) Truc, B. _et al._ Light-induced metastable hidden skyrmion phase in the mott insulator cu2OSeO3.  _Advanced Materials_ 35 (2023).  URL https://doi.org/10.1002/adma.202304197.