Technological progress often stems from the 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 the environment and climate. This is also true for electronics: it has become clear that the 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 these two directives on equal footing 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 nontrivial 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 racetrack memories [8–10] and quantum computation devices [11–14]. However, this 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 open spin–electron systems (Figure 1). Our method is an extension of the approach introduced in [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 (LLG) equation for the classical spins. The former equation can be derived from the general theory of nonequilibrium Green’s functions using the so-called generalized Kadanoff–Baym ansatz (GKBA) [17–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 ₀, after which the velocity is found to be a linear function of the current I − I ₀. The situation is found to be qualitatively different in the 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 LLG equation are not always adequate and that the dynamics of the electrons must be explicitly taken into account.

This paper is structured as follows: Section 2 briefly reviews previous approaches to coupled spin–electron dynamics, and Section 3 introduces the system and Hamiltonian to be considered and discusses the coupling between the itinerant electrons and external reservoirs. The spin and electron equations of motion are presented in Section 4 and Section 5, and Section 6 introduces an approximate wide-band limit (AWBL) as a numerically efficient way to propagate the equations of motion in the presence of large central regions connected to external reservoirs. Section 7 presents some observables used to interrogate the skyrmion content of the spin configuration. In Section 8, we use the AWBL to investigate the skyrmion motion induced by current densities in the itinerant electron system, and in Section 9, we consider skyrmion motion in the presence of magnetic disorder. Finally, in Section 10, we conclude with a discussion of experimental signatures and possible material platforms for which our results are of relevance.

Magnetic skyrmions are made up of localized magnetic moments, typically arising from a large Hund’s coupling J favoring a high-spin state of 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 ¹ . 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 to 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]. However, in certain cases, such as at an interface between metallic thin films and a material with large spin–orbit coupling, the effective Dzyaloshinskii–Moriya interaction can be significantly enhanced and lead to skyrmions with radii on the order of a few nanometers [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 approximately 10 nm in radius in two-dimensional systems) are usually based on an atomistic description of classical spins, with external fields such as electromagnetic radiation or electronic current densities included as non-dynamical variables. For example, the motion of skyrmions in response to an external current density can be described via an effective equation for the skyrmion’s center of mass, the so-called Thiele equation [36], assuming that the form of the skyrmion is rigid. Based on 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 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 equation, which go beyond the semiclassical Thiele’s description to study the dynamics of skyrmions [39–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 [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 description of skyrmion dynamics (as well as more general spin textures) has been addressed in several contexts [16, 37, 38, 43–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 photo-generation of skyrmions resulting from laser excitation can occur on much shorter time scales than previously thought.

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

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 Figure 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 are coupled 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 the density set by the chemical potential. Furthermore, 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 H t = H C t + H R t + H C R , (1) where a possible time dependence is explicitly indicated both for the central region and for the reservoirs. This is necessary to initiate the dynamics of the system. We now discuss each contribution to H separately.

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 S ̂ i → ⟨ S ̂ i ⟩ ≡ S i = S n 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 semi-classical 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 d n i d t = − n i × ∑ j S J i j n j + D i j × n j + B + S K n i z z ̂ + g ⟨ s ̂ i ⟩ , (11) where the last term provides 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 amounts to adding an extra term α n _i × (∂ n _i /∂t) to Eq. 11. The equations of motion for the itinerant electrons are discussed in detail in the next section.

To describe the dynamics of the electrons, we use a nonequilibrium Green’s functions (NEGFs) description within the GKBA. NEGFs are a general and powerful approach to nonequilibrium phenomena [51–55], to describe 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 G i j σ σ ′ z , z ′ = − i T γ c i σ H z c j σ ′ H z ′ † 0 . (12) 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 is the equilibrium Hamiltonian, and Z = tr e ^−βH is the canonical partition function. In the low-temperature limit β → ∞, the ensemble average reduces to a ground-state expectation value. The operator c i σ 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), where 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 mentioned 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 . Therefore, we should write G i j σ σ ′ ( z , z ′ ) = G i j σ σ ′ ( z , z ′ ; { S i } ) , 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 Green’s function, with the implicit understanding that all quantities are matrices in the 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] i d d z − h z G z , z ′ = δ z , z ′ + ∫ γ ⅆ z ̄ Σ z , z ̄ G z ̄ , z ′ . (13) Here, h(z) is the time-dependent mean-field Hamiltonian and the self-energy Σ(z, z′) contains 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 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 by Σ _c and Σ_emb, respectively, 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 a wide-band limit (WBL) approximation to the embedding self-energy, obtained when the hopping amplitude inside the leads tends to infinity, while the ratio t R α 2 / 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 to the assumption that the density of states of the leads is constant over the bandwidth of the central system. Performing the WBL approximation results in the following expression for the embedding self-energy [20, 30, 57, 58]: Σ emb < t , t ′ = i ∑ α Γ α s t s t ′ e − i ∫ t t ′ V α t ̄ d t ̄ ∫ d ϵ 2 π f ϵ − μ e − i ϵ t − t ′ Σ emb R t , t ′ = − i 2 δ t − t ′ s t s t ′ ∑ α Γ α , 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 the presence of the leads, V _α (t) is the time-dependent bias in lead α, and f(ϵ) is the Fermi–Dirac distribution.

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 an 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 G A t , t ′ = i θ t ′ − t T ̄ e − i ∫ t ′ t d t ̄ H e + H s − e t ̄ + i Γ / 2 ≈ i θ t ′ − t e − i t − t ′ H e + i Γ / 2 . (16) 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 Supplementary Material S1). 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.

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 ² and the central region can be approximately compactified to the torus T ² (assuming a ferromagnetic ordering along the edges), a topological charge measuring the winding number of the map S: T ²↦S ² can be defined. On a lattice, it can be shown that a suitable definition of this topological charge is given by [60] Q = 1 4 π ∑ j k l ∈ Δ Ω jkl , (17) 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 exp i Ω jkl / 2 = ρ jkl − 1 1 + S j ⋅ S k + S k ⋅ S l + S l ⋅ S j + i η jkl S j ⋅ S k × S l , ρ jkl = 2 1 + S j ⋅ S k 1 + S k ⋅ S l 1 + S l ⋅ S j , (18) and the function η jkl = s g n [ ( 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 R c m = 1 M ∑ i x = 2 N x − 1 ∑ i y = 2 N y − 1 ρ i x , i y i x i y . (19) The expressions for M and ρ i x , i y are respectively given by M = ∑ i x = 2 N x − 1 ∑ i y = 2 N y − 1 ρ i x , i y , (20) ρ i x , i y = | Ω i x , i y , i x + 1 , i y , i x + 1 , i y + 1 + Ω i x , i y , i x , i y + 1 , i x + 1 , i y + 1 | . (21) In Eqs 19, 20, 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’s center of mass. The rationale behind these oscillations and the specific way we perform the site exclusion are further detailed in Section 8.1.

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²(πt/2τ). After time t = τ, a bias is applied, here considered to be of the form v α t = θ t − τ V if  α = L  and  σ = ↑ − V if  α = R  and  σ = ↑ 0 otherwise . (22) This bias generates a spin current through the central system. In the simulations below, we fix the energy unit by setting the hopping t _ij = t _s = −1 and take ℏ = 1.

Knowing the trajectory of the skyrmion in the base case, we can now investigate the effects of placing impurities in its path. In the panels of Figures 5–7 showing the spin configuration, the positions of the impurities are marked in green. Here, we restrict ourselves 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 = 1,400 in order not to disturb the system before the bias is applied at t = 1,500.

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 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 NEGFs, coupled to classical spins governed by an LLG equation.

To describe the skyrmion motion in large systems, we introduced an approximate approach to include the effect of the electronic reservoir coupled to the central system. This so-called AWBL shows a satisfactory level of agreement with the full 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 skyrmions. 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 to 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 competing with the 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 validating methods that, still in a spin-only 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 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 [16, 63]. Apart from its importance for skyrmionic systems, the presented method will, therefore, be of large relevance to describe the coupled dynamics of spins, electrons, and lattice vibrations in spintronic devices and for material engineering, thereby helping facilitate a microscopic understanding of driven magnetic materials.