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

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 [1].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 [2].
Pursuing on equal footing these two directives is the core aim of spintronics [3]: 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 FIG. 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.magnetic textures, creating unprecedented possibilities for spintronic applications [4,5].
A notable example in this respect is provided by magnetic skyrmions [6].These are topologically nontrivial spin textures stabilized by a competition of exchange, Dzyaloshinskii-Moriya interactions (DMIs) and magnetic anisotropies [7].With their compact size, topological protection, and non-intensive energy requirements for manipulation, skyrmions are of great potential interest to realize race track memories [8][9][10] and quantum computation devices [11][12][13][14].This, however, requires efficient ways of writing, deleting, and manipulating skyrmions on short time scales and with high spatial precision, via electronic spin currents.
Here we describe a scheme to perform large scale simulations of the intertwined dynamics of interacting and arXiv:2312.12201v1[cond-mat.mes-hall]19 Dec 2023 open spin-electron systems (Fig. 1).Our method is an extension of the approach introduced in Refs.[15,16], 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) [17][18][19][20][21][22][23][24][25][26][27][28][29][30], 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 [31], 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.[32]).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 spinelectron 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 necessary1 .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 α of the spin-orbit interaction is typically of the order α/t ∼ 0.01 − 0.1, resulting in a spin spiral wavelength λ ∼ 10 − 100 nm [9,34].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 [35].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 [36], 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 [31], and to include the effects of electrons and external fields via a generalized force [37,38].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 [39][40][41][42].
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. [15], 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 [16,37,38,[43][44][45].For example, it was recently demonstrated [16] 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 photogeneration 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 × 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 [15].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 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.

A. 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 Here c † iσ creates an itinerant electron at site i with spin projection σ, the hopping amplitude between nearestneighbor sites i and j (denoted by ⟨ij⟩) is given by t ij , α ij accounts for spin-orbit interactions, and B = Bẑ is an external magnetic field along the z-axis.The electronic spin operator at site i is defined by ŝi = σσ ′ c † iσ τ σσ ′ c iσ ′ , where τ denotes the vector of Pauli matrices.The parameters J ij = J ji and D ij = −D ji give the exchange and Dzyaloshinskii-Moriya interaction (DMI) [46,47] between spins Ŝi and Ŝj , and K quantifies an easyaxis 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 [35], 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 D ij can be of different form.Assuming the overall magnitude of the interaction is fixed at |D ij | = D, the spatial dependence typically takes one of the following two forms Here d ij is the vector between lattice sites i and j, and ẑ 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 Blochtype DMI that commonly arises in systems with noncentrosymmetric crystal structures, gives rise to Blochtype (spiral) skyrmions.Even so, many bulk materials that support Bloch-type skyrmions can be fabricated as thin film samples, while still hosting skyrmions [34].
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 where the parameter α 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 D ij ∝ α ij [47].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ϕij whose phase encodes the interaction with an external electromagnetic field [16].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.

B. The reservoirs and their coupling to the central region
The leads are taken to be two-dimensional and semiinfinite (see Fig. 1), and are described by the Hamiltonian Here α ∈ {L, R} denotes the left (L) and right (R) lead, σ labels the spin projection, and the parameter t Rα is the hopping amplitude within lead α.The hopping amplitude is related to the nominal band width by W α = 4t Rα , and v α σ (t) is a time-dependent bias measured from the chemical potential µ.The leads interact with the central region via the Hamiltonian H CR = ασ H CR α,σ , where 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 ∼ 100 nm [48], 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 Ŝi → ⟨ Ŝi ⟩ ≡ S i = Sn i and |n i | = 1.This limit is exact for S → ∞ [49,50], and considered to be a suitable approximation already when S > 1 [35].Taking the classical limit leads to a semiclassical approximation for the coupled spin and electron subsystems [15,16], where the quantum electronic Hamiltonian depends parametrically on the classical variables 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 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 → J ij and similarly for 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 αn i × (∂n i /∂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 [51][52][53][54][55], which describes the realtime 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 Here the brackets ⟨•⟩ 0 denote an ensemble average with respect to the thermal density operator ρ = e −βH /Z, where β = 1/k B T is the inverse temperature, H the equilibrium Hamiltonian and Z = tr e −βH the canonical partition function.In the low-temperature limit β → ∞, the ensemble average reduces to a ground state expectation value.The operator c H iσ (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 [52], and all operators inside the brackets are time-ordered on the contour by the operator T γ .As before, the indices i and j run over all lattice sites, while σ and σ ′ 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 , but for notational simplicity we keep the dependence on 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 ′ ).The equations of motion for the one-particle Green's function, the so-called Kadanoff-Baym equations, can be written as [51] Here h(z) is the time-dependent mean-field Hamiltonian, and the self-energy Σ(z, z ′ ) 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 [56].Denoting the correlation and embedding self-energies respectively by Σ c and Σ emb , the total self-energy can be written as Σ = Σ c + Σ emb .We note that in the present case the electronic Hamiltonian contains no interaction terms, and therefore Σ 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 2 Rα /t α 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 [20,30,57,58] where the superscripts denote the so-called lesser (<) and retarded (R) components of the self-energy [54].Here s(t) is a smooth function used to equilibrate the central system in presence of the leads, V α (t) is the time-dependent bias in lead α, and f (ϵ) is the Fermi-Dirac distribution.
A. Generalized Kadanoff-Baym ansatz Provided that we have full knowledge of the self-energy Σ, the Kadanoff-Baym equations are an exact reformulation of the many-body problem.In practice however, Σ 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 [9], 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 timestep scaling can be achieved by employing the so-called generalized Kadanoff-Baym ansatz (GKBA) [17], 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 [17][18][19][20][21][22][23][24][25][26][27][28][29][30].Within NEGF-GKBA and its recent reformulation as a time-linear scheme [25,26], it is currently possible to simulate the long-time dynamics of electronic systems with a basis size on the order of 100 orbitals [28].
The GKBA for electronic degrees of freedom is achieved by the following factorization [17], where the superscript A denotes the advanced component of the Green's function.Using transformation rules introduced by Langreth [59], we can project Eq. 13 from the Keldysh contour to the physical time axis, and obtain a dynamical equation for G < (t, t ′ ).When taking t ′ → t and using the GKBA, the equation for G < (t, t ′ ) reduces equation of motion for the one-particle density matrix ρ(t) = iG < (t, t) of the form.
To close the equation for ρ we further assume that , with the interaction between the leads and the central region adiabatically turned on before the bias is applied.Here Γ = α Γ α accounts for the presence of the leads [20].

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 ρ 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 [29,30,57,58], 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 With this approximation the advanced Green's function becomes a function of the time difference, G A (t, t ′ ) = G A (t − t ′ ).As a consequence, the entire collision integral becomes a function of t−t ′ , 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 S : T 2 → S 2 can be defined.On a lattice, it can be shown that a suitable definition of this topological charge is given by [60] where the sum is over the triangulated lattice and ∆ is the set that contains the indexes for the lattice sites for every triangle in the lattice.The solid angle Ω jkl is defined by and the function η jkl = sgn[(S j • (S k × 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π of the unit sphere gives the integer winding number (for appropriate boundary conditions).The solid angle Ω jkl serves as an indicator of the extent to which the spins twist.Specifically, within the core of a skyrmion, Ω 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 The expressions for M and ρ ix,iy are respectively given by ρ ix,iy = Ω (ix,iy),(ix+1,iy),(ix+1,iy+1) + Ω (ix,iy),(ix,iy+1),(ix+1,iy+1) .
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 A.

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 < τ using the contact function s(t) = sin 2 (πt/2τ ).After time t = τ a bias is applied, here considered to be of the form 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 ℏ = 1.
A. 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 × 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 ∈ [0, 800] measured in units of ℏ/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 leaddevice connection γ.
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 α = 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.
The purpose of using a large α 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 B and afterwards, where we consider skyrmion dynamics in central regions of 20 × 40 sites), the damping is set small (α = 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 R cm .Accordingly, in the results in the next sections, R cm is calculated using Eq.20, i.e without including the contribution of the peripheral sites of the 20 × 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 γ, 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 γ 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.

B. 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 × 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 γ = 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 [9].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 previ-FIG.3. Benchmarks for the approximate wide band limit.(a) Initial relaxed configuration of a 4 × 4 spin texture interacting with a corresponding layer of electrons.(b − d) Spin current I R,↑ and electronic spin density ρ 11↑ of the lattice site encircled in green in (a), for different values of γ 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 α.The parameters used are B = 0.05, K = 0.6, α = 3, αR = 0, J = 0.5, D = 0.4, τ = 800, µ = 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 γ are given in each panel.
ous studies [61,62] In particular, it was recently proposed that the magnitude of the skyrmion mass is affected by the interactions between localized moments and free electrons [37,38], 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 → 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.
A. 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 xdirection, 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 ∼ 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 ∼ 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.

B. 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 × 4 spins, a rotation of the xycomponents of the spins away from the ferromagnetic background can be seen in an area extending up to 12×12 sites.This indicates that the effective skyrmion size is much larger than the central 4 × 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.

C. Case 3: A larger impurity configuration
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 ∼ 2100, the skyrmion shrinks and expands in a quasi steady state.

D. 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 [63].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 wideband 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 [37,38].
Finally, we note that the methods presented here open the door for simulations of the intertwined spinelectron 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 [16,64].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.
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.

FIG. 2 .
FIG. 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 µ, and are connected to a central region of 20×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↑ = −V R↑ and V L↓ = V L↓ = 0) is applied in the leads.