Phenomenology of current-induced skyrmion motion in antiferromagnets

We study current-driven skyrmion motion in uniaxial thin film antiferromagnets in the presence of the Dzyaloshinskii-Moriya interactions and in an external magnetic field. We phenomenologically include relaxation and current-induced torques due to both spin-orbit coupling and spatially inhomogeneous magnetic textures in the equation for the N\'eel vector of the antiferromagnet. Using the collective coordinate approach we apply the theory to a two-dimensional antiferromagnetic skyrmion and estimate the skyrmion velocity under an applied DC electric current.


Introduction
For decades ferromagnets have been the main components of spintronic devices [1]. Antiferromagnets (magnetically ordered materials with compensated magnetization) have long remained in their shadows [2,3], even though they possess properties that make them appealing as potential alternatives of ferromagnets and as next generation data storage devices. The insensitivity of antiferromagnets to external magnetic fields makes them more robust against magnetic perturbations, they operate on faster timescales possibly enabling ultrafast information processing, and antiferromagnetic metals, alloys and semiconductors are not limited to combinations of only a few elements like Fe, Co, Ni [4]. On the other hand, the insensitivity to external magnetic fields makes antiferromagnets also much harder to manipulate and control [3]. In recent years, several breakthroughs were achieved in overcoming this obstacle. The anisotropic magnetoresistance (AMR) effect was proposed [5] and utilized [6] to electrically detect antiferromagnetically ordered states. Another important step towards antiferromagnetic spintronics was the prediction [7] and subsequent observation [8] of Néel spin-orbit torques in a certain class of antiferromagnets that can electrically manipulate the antiferromagnetic Néel vector.
Independently of the developments with antiferromagnets, skyrmions have been gaining momentum in the field of data storage [9]. Magnetic skyrmions, for example, are topological windings of the magnetization on the nanoscale that appear in noncentrosymmetric materials. In ferromagnets, skyrmions are promising candidates as information carriers for future information-processing devices and have been intensely studied in recent years [9,10,11,12,13]. One of the key properties that makes them attractive is the very small current densities that are needed to set them in motion [10].
In contrast, not much is known about skyrmions in antiferromagnets. While their existence has been predicted by Bogdanov and colleagues [14] (see also [15]), there are no experimental observations. The differences and possible advantages of antiferromagnets over ferromagnets lead to the question: how will skyrmions interact with currents in antiferromagnets?
To utilize skyrmions, it is necessary to know how to create, manipulate and detect them in magnetic thin-film nanostructures. Initial numerical analyses have been performed on two-dimensional antiferromagnetic films [16,17] focusing on the manipulation of skyrmions by electric currents. In both treatments the skyrmion dynamics were studied by solving the equations of motion for the two sublattices numerically where damping and current-induced torques were implemented in the same spirit as for ferromagnets. In this article, we use a phenomenological and analytical approach to gain further insights.
The latter analytical approach has been applied in recent works to study the magnetization dynamics of antiferromagnets under the influence of electric currents [18,19,20,21]. However, so far the effect of the inversion symmetry breaking inducing the Dzyaloshinskii-Moriya (DM) interactions has not been taken into account.
where the second term in the last equation is neglected because it is quadratic in the driving forces [14].
Here we apply the results for the general magnetization dynamics to a rigid antiferromagnetic skyrmion to analyze its current-induced motion. We use the ansatz n(r, t) = n sk (r R(t)), where n sk (r) is the static profile of an isolated skyrmion and R(t) is the skyrmion position. As collective coordinates we take {⇠ i } = {R x , R y }. After multiplying (8) by @n/@⇠ ↵ for ↵ = x, y and integrating over space, the equation of motion transforms to The coe cients read where represents the friction term and h e↵ , h denote the transverse terms that originate from the external magnetic field (similar to a Lorentz force). The constants I j are dimensionless integrals determined by the skyrmion profile and we discuss them later. The appearance of the e↵ective mass m e↵ is the main di↵erence compared to ferromagnetic skyrmion motion resulting from the di↵erent nature of the magnetization dynamics in antiferromagnets.

Making sure the constraints are fulfilled
We need to make sure that n 2 = 1 and n · M = 0.
Obtaining the skyrmion velocity Collective coordinate approach, as done previously The idea is to parametrize the Néel order parameter by collective coordinate ⇠ i At the same time, however, the Néel vector depends on the spherical angles n(t) = n ⇣ ✓(t), ⌘ Which one exactly of the spherical angles contains the time-dependence? Theṅ Subsequently multiply the whole equation by @n @⇠j · . . . and integrate over space. We choose the skyrmion center coordinates Word of caution regarding the derivatives @n @⇠i In all the subsequent calculations, I have equated for i = x (⇠ x = R x )and the equivalent for i = y (⇠ y = R y ). However this is only true for the isolated skyrmion, that is centered at the origin n sk (r). Once the skyrmion is set in motion, then the actual Neel vector becomes n = n(r R(t)). Consequently, the partial derivatives need to be replaced by for i = x and the equivalent for i = y. Thus, all the partial derivatives of the skyrmion Neel vector need to be replaced by a negative sign in the subsequent section.

11
Making sure the constraints are fulfilled We need to make sure that n 2 = 1 and n · M = 0. Obtaining the skyrmion velocity Collective coordinate approach, as done previously The idea is to parametrize the Néel order parameter by collective coordinate ⇠ i At the same time, however, the Néel vector depends on the spherical angles n(t) = n ⇣ ✓(t), ⌘ Which one exactly of the spherical angles contains the time-dependence? Theṅ Subsequently multiply the whole equation by @n @⇠j · . . . and integrate over space. We choose the skyrmion center coordinates ⇠ i = R x (t), R y (t).
Word of caution regarding the derivatives @n @⇠i In all the subsequent calculations, I have equated for i = x (⇠ x = R x )and the equivalent for i = y (⇠ y = R y ). However this is only true for the isolated skyrmion, that is centered at the origin n sk (r). Once the skyrmion is set in motion, then the actual Neel vector becomes n = n(r R(t)). Consequently, the partial derivatives need to be replaced by for i = x and the equivalent for i = y. Thus, all the partial derivatives of the skyrmion Neel vector need to be replaced by a negative sign in the subsequent section.
ime, however, the Néel vector depends on the spherical angles n(t) = ich one exactly of the spherical angles contains the time-dependence? Theṅ ultiply the whole equation by @n @⇠j · . . . and integrate over space. We choose enter coordinates ⇠ i = R x (t), R y (t).
tion regarding the derivatives @n @⇠i quent calculations, I have equated R x )and the equivalent for i = y (⇠ y = R y ). However this is only true for yrmion, that is centered at the origin n sk (r). Once the skyrmion is set in e actual Neel vector becomes n = n(r R(t)). Consequently, the partial d to be replaced by he equivalent for i = y. Thus, all the partial derivatives of the skyrmion d to be replaced by a negative sign in the subsequent section. 11 so that, in total,

Correction to the equation of motion
Thus, the equation for the Néel vector needs to be extended to satisfy the normalization constraints Obtaining the skyrmion velocity Collective coordinate approach, as done previously The idea is to parametrize the Néel order parameter by collective coordinate ⇠ i Figure 1. Schematic representation of the skyrmion motion in an antiferromagnet driven by an electric current j x and an external magnetic field H ẑ, as described by (1). The friction force is denoted by ΓṘ and the longitudinal current-induced force by ∆j. The combination of these forces leads to a skyrmion motion with the velocitẏ R sk . For simplicity, a static skyrmion shape is depicted, without taking into account the shape changes that the magnetic field and the electric current induce.
Here, starting from general principles and symmetry considerations, we construct the equations describing the macroscopic antiferromagnetic magnetization dynamics in the presence of a DC current and an external magnetic field. In particular, we i) take into account the effect of the DM interactions corresponding to a particular symmetry class of magnetic materials, ii) incorporate phenomenologically the spin-orbit torques of an antiferromagnetic system and iii) study the role of the external magnetic field and its effect on the skyrmion motion. We demonstrate that the latter has, in fact, no effect on the skyrmion motion. It does, however, affect its shape. Our main result is an equation of motion for the position R of the skyrmion in the thin film given by Here, m eff is the effective skyrmion mass, and Γ a friction constant. The external magnetic field does not appear in the equation. The dissipative spin-orbit and spintransfer torques (described by ∆) lead to a longitudinal current-induced force on the skyrmion. In figure 1 we illustrate the resulting skyrmion motion. We emphasize that we included both homogeneous and inhomogeneous current-induced torques in our analysis, thus obtaining a more general form of the equation of motion than considered in [16,17]. The remainder of the paper is structured as follows. In section 2, we discuss the phenomenological model and the resulting equations describing the current-induced antiferromagnetic dynamics. We discuss the current-induced torques and we derive a closed equation of motion for the antiferromagnetic order parameter. In section 3 we reformulate the result into an equation for the position of a skyrmion by using the collective coordinate approach. Finally, we provide an estimate for the skyrmion velocity.

Phenomenological model
In this section we analyze the general magnetization dynamics of the antiferromagnets of interest. We first discuss the magnetic energy of the model and then derive the equation of motion for the Néel vector.

Energy functional
We consider a two-sublattice antiferromagnet within the exchange approximation with the sublattice magnetizations M 1 −M 2 and |M 1 | = |M 2 | = M s , where M s is the saturation magnetization. Further, we consider a noncentrosymmetric lattice of the crystallographic class C nv [14]. This class encompasses also two-dimensional films which have structural inversion asymmetry along theẑ-direction, for example, due to the presence of an interface.
It is most convenient to formulate the theory in terms of the antiferromagnetic order parameter (also called the Néel vector) n = (M 1 −M 2 )/2M s and the total magnetization m = (M 1 + M 2 )/2M s . Our phenomenological approach is based on the exchange approximation [22], which requires rotational invariance of the magnetization vectors and invariance of the theory with respect to an exchange of the two sublattices, i.e., under the transformations n → −n and m → m. The magnetic energy that follows from these considerations is ‡ [14,22] The first and the second terms describe the inhomogeneous and homogenous exchange interaction with the constants A and H exc , respectively. The uniaxial anisotropy is parameterized by the constant H an > 0 and the external magnetic field is denoted by H. The remaining terms describe the DM interactions, where d ẑ and D represent the homogeneous and inhomogeneous parts, respectively. The latter part of the DM interactions, also called a Lifshitz invariant, is the main ingredient needed to stabilize a skyrmion in this system and inhomogeneous textures in general [14].
The antiferromagnetic vectors obey the constraints n 2 = 1 and n · m = 0. Throughout this work we make the assumption that the homogeneous exchange interaction H exc is the dominant energy scale, so that H exc d, H an . For typical values of the external magnetic field that do not destroy the antiferromagnetic order the exchange constants dominates the field too, H exc H [6].

Landau-Lifshitz-Gilbert equations
In this section we derive the equations of motion for the Néel order parameter for a timeindependent magnetic field and electric current. The Landau-Lifshitz-Gilbert equations to leading order in H exc are given by [18,19,23] n − γ 2M s n × δf δm + τ n , ‡ The energy functional seems to satisfy the C ∞v symmetry group, however, the discrete symmetry origin of the antiferromagnetic vectors needs to be respected, which lowers the symmetry to C nv .
where γ is the gyromagnetic ratio, α G § is a phenomenological Gilbert damping coefficient and τ n,m represent the respective current-induced torques to be discussed in section 2.3. The functional derivatives of the energy density f are 1 2M s The static magnetization of this system without electric currents is found by taking the cross product of the first line of (3) with n and identifying the time-independent components under the constraints n 2 = 1 and n · m = 0, which yields [14,15] The presence of the homogeneous DM interactions term d shows that even in the absence of an external magnetic field the total magnetization is non-zero. However, this term does not contribute to the motion of an antiferromagnetic texture within the approximations we consider and we neglect it for the remainder of the paper.

Current-induced torques
We follow a phenomenological approach to derive the current-induced torques. It is based on the Onsager reciprocity relations, which, in the case under consideration, relate the process of inducing charge currents by a time-varying magnetic texture to the effect that charge currents have on the magnetization dynamics. Following [19], to lowest order in the spatial gradients and the magnetization m and zeroth order in spin-orbit coupling we find three contributions towards the magnetically pumped charge density j pump /σ that obey the symmetries of the system (rotation and n → −n): ηM s /γ n · (ṁ × ∂ i n), βM s /γṅ · ∂ i n and ζM s /γ n · (ṅ × ∂ i m) [19]. After applying the Onsager relations these are transformed and result in the torques Here, the terms parameterized by the coefficients η, ζ describe reactive spin-transfer torques, whereas the term with β describes a dissipative spin-transfer torque. Furthermore, the inversion symmetry breaking gives rise to another set of torques even in homogeneous systems, which do not involve gradients of the antiferromagnetic vectors. Following the same approach as above [24], we find the pumped charge currents § In the notation of [19], this is the G 2 damping constant.
Antiferromagnetic skyrmion motion 6 that are lowest order in the spin-orbit coupling: C 1 M s /γẑ ×ṁ, C 2 M s /γẑ × (n ×ṅ) and C 3 M s /γẑ × (m ×ṁ), which after applying the Onsager relations lead to Here, the terms proportional to C 1 are field-like (reactive) spin-orbit torques and the terms proportional to C 2 , C 3 are the anti-damping (dissipative) spin-orbit torques.
In ferromagnetic systems similar torques have been discussed in [25,26] and in antiferromagnets a microscopic analysis has been performed in [7]. We emphasize that the spin-orbit torques and the spin-transfer torques differ in their nature. Whereas in the latter the free electrons are polarized by the local magnetization while moving through the texture and interact with it after being polarized, in the spinorbit torques the polarization is due to the spin-orbit coupling in the system and not due to the magnetization. In that sense, the spin-transfer torques are a result of a nonlocal interaction between the electrons and the magnetic moments, while the spin-orbit torques are local.
In the later steps of the calculation, presented below, the form of both the spintransfer and the spin-orbit torques will be simplified by retaining only the leading order terms in H exc .

Equation of motion
Writing out the torques explicitly, the Landau-Lifshitz-Gilbert equations (3) becomė Here, only terms to leading order in H exc have been kept, apart from the η and C 1 ones in the second line, which we retained in order to keep the constraints n 2 = 1 and n · m = 0 fulfilled. The term containing the functional derivative of the energy density with respect to m is of subleading order in H exc , however, it is of the same order as the left hand side after substituting (9) below and needs to be kept as well.
An expression for the total magnetization can be obtained from the equation forṅ: where m 0 now does not contain the homogeneous DM interactions contribution d, as discussed in section 2.2. With this, we are able to write a closed equation for the Néel Antiferromagnetic skyrmion motion 7 order parameter n ×n = H shape + H forces .
Here, we have grouped the right hand side into terms that determine the shape of the antiferromagnetic texture and terms that induce or affect its motion. The fields are given by and In deriving the equation for the Néel order parameter only terms up to linear order inṅ and j have been kept. An example of a term of higher order that has been omitted is n × (j · ∇)ṅ.
Equation (10) is an important result of this work and describes the magnetization dynamics of a uniaxial antiferromagnet with inversion symmetry broken along theẑdirection under the influence of an external time-independent magnetic field and DC electric current.

Skyrmion motion
In the previous section we analyzed the magnetization dynamics of a uniaxial antiferromagnet in the presence of electric currents and a time-independent external magnetic field. Here, we focus on the translational motion of a magnetic skyrmion, rewrite the equations of motion by using the collective coordinate approach, and obtain an estimate for the skyrmion velocity as a result of the electric currents.

Collective coordinates
Previous work has shown that the dynamics of magnetic textures in ferromagnets and antiferromagnets can often be described by only a few variables [15,20]. The approach necessitates the choice of a finite set of collective coordinates ξ i (t) which are used to specify the time evolution of the Néel order parameter n(r, t) = n(r, {ξ i (t)}). In particular, we usė where the second term in the last equation is neglected because it is quadratic in the velocities, which are assumed to be small [20].
Here, we apply the approach to the translational motion of an antiferromagnetic skyrmion to analyze its current-induced dynamics. We assume that the skyrmion profile is composed of a static, cylindrical and rigid component n sk and motion-and currentinduced corrections δn that break the cylindrical symmetry (see section 3.2). For the time evolution we use the ansatz n(r, t) = n(r − R(t)), where n = n sk + δn and R(t) is the skyrmion position. As collective coordinates we take {ξ i } = {R x , R y }. After multiplying (10) by n × ∂n/∂x α for α = x, y and integrating over space, the equation of motion for the skyrmion position to leading order in the electric currents becomes The coefficients read where Γ represents a friction term and ∆ characterizes the effect of the dissipative current-induced torques. The dimensionless constant I is determined by the skyrmion profile, which we discuss later. The characteristic lengthscale x 0 is the domain wall width of the system and is given in the following section. The dependence of the effective mass m eff on the exchange constant H exc is the main difference compared to ferromagnetic skyrmion motion and results from the different nature of the magnetization dynamics in antiferromagnets.
In deriving (13) we have considered both homogeneous and inhomogeneous currentinduced torques, as well as an external magnetic field applied in theẑ-direction. Our results thus paint a richer picture of the current-induced antiferromagnetic skyrmion motion than discussed recently. References [16,17] predicted longitudinal currentinduced forces on the skyrmion position, as opposed to the ferromagnetic skyrmion motion, where the nonzero magnetization always leads to a transverse force. Reference [16] dealt with homogeneous torques only, whereas in [17] only gradient torques have been considered. In both references no magnetic field has been included. We find, however, that even in the presence of an applied field the skyrmion motion remains longitudinal. This is further substantiated by the findings of [15,27], where it is shown that gyroscopic forces (which are linear in the magnetic field) are not present in antiferromagnets for objects that exhibit the topology of skyrmions.
For an electric current applied in thex-direction, an expression for the longitudinal velocity v sk can be readily obtained by the steady-state solution of (13), which yields The choice of this factor comes from general considerations for the conserved quantities in the system. For the translational motion, the relevant quantity is the momentum [15]. Note that this also agrees with the approach taken in [20]. This is the velocity corresponding to the zero-field scenario considered in [16,17]. Before we proceed with its estimate, we need to analyze the skyrmion shape in more detail.

Skyrmion profile
To calculate the constant I, we need to determine the profile of the antiferromagnetic skyrmion. The model in (2) allows for skyrmion solutions, as long as the external field is applied along theẑ-direction or is zero [14,15]. We assume the external magnetic field to be the dominant contribution towards H eff , so that the current-induced contribution towards the effective field does not destroy the skyrmion (that is, H C 1 M s j/2γ). The skyrmion profile n = n sk + δn is determined from the steady state of (10) in the absence of dissipative and damping terms n × δf δn Here, the static component n sk solves the equation with both j,ṅ = 0, whereas the corrections δn originate from a non-zero dynamic termṅ and the current-induced effective field C 1 M s /2γ(ẑ × j). We are considering slow skyrmion dynamics, therefore the corrections can be assumed small, |δn| |n sk |. The profile equation needs to be solved numerically [14]. For the purposes of the present work, we restrict the further analysis to the static component n sk (the corrections δn will not lead to a qualitative change in the skyrmion velocity estimate). It is convenient to rewrite (16) into spherical coordinates for the Néel vector, n sk = (sin θ cos ψ, sin θ sin ψ, cos θ), and cylindrical coordinates for the spatial variables, r = ρ(cos φ, sin φ). For the crystallographic class under consideration, a skyrmion solution exists when the Néel order parameter has the same azimuthal direction as the cylindrical spatial vector (that is ψ = φ, see figure 2) [14]. The corresponding equation becomes Here x 0 = A/ |H an | is the characteristic lengthscale (domain wall width) of the system, H 0 = |H exc H an | the spin-flop field, D 0 = 4/π A |H an | the threshold value of the inhomogeneous DM interactions constant that stabilizes modulated magnetic structures at zero magnetic field andρ = ρ/x 0 the rescaled radial coordinate. Figure 3 shows the skyrmion profile for the particular choice of D/D 0 = 0.9 and H/H 0 = 0.3. The boundary conditions used are θ(ρ = 0) = π and θ(ρ → ∞) = 0, where the latter condition is implemented by a shooting method. With this, we are in a position to calculate numerically the integral I. The result is given in Appendix B.

Skyrmion velocity
Now, we are in a position to give an estimate for the magnitude of the longitudinal skyrmion velocity (15).
We estimate the spin-orbit torque coefficient to be C 2 3.4 × 10 −3 m 2 A −1 s −1 (see Appendix C and [7]). The Gilbert damping parameter is α G 0.01 [19]. The characteristic length of the system is of the order of the skyrmion size, which is typically x 0 10 −8 m [14,19,20]. Typical experimentally used current densities in ferromagnets are of the order of j 10 11 Am −2 [25]. With this, we estimate v sk,C 2 255 ms −1 .
In contrast, [16] predicts a skyrmion velocity of ∼ 1700 m/s. In that work, homogeneous torques of a similar form, but of different origin have been considered. The homogeneous torques there arise from a spin-polarized current injected vertically into the system, whereas the homogeneous torques in the present work are due to the spin-orbit coupling in the antiferromagnet. Using the values that they provide (see table 1), we arrive at v sk,C 2 500 ms −1 , which has the same order of magnitude as the result of [16]. In (8) the dissipative spin-transfer torque coefficient β has dimensions of m 3 A −1 s −1 . Its dimensionless counterpartβ is obtained byβ = βne, where n is the electron density and e the electron charge. Typically, in ferromagnetic systems this value is taken to be of the order of the Gilbert damping,β α G [28]. Typical metallic electron densities are of the order of n 10 29 m −3 and e = 1.6 × 10 −19 As, so that for these parameters v sk,β 5 ms −1 .
This component of the velocity corresponds to the skyrmion velocity in [17]. Our estimate agrees with the findings of that work for the same choice of parameters (see Table 1. Estimated skyrmion velocities using the parameters in the present work and in [16,17]. In all cases the velocities v sk,C2 and v sk,β are estimated according to the expressions in (15). The numerical results for v are taken from [16] and [17], respectively. This value has been calculated from the expression that the authors provide in [16] for their Slonczewski-like spin-transfer torque coefficient β = |h/(µ 0 e)| P/(2dM s ), multiplied by the gyromagnetic ratio |γ| = 2.211 × 10 5 mA −1 . Here,h is the reduced Planck constant, µ 0 the vacuum permeability, e the electron charge, P = 0.4 is the polarization rate of the spin-polarized current, d = 0.4 nm the film thickness and M s = 290 kA m −1 is the saturation magnetization; For comparison, we focus only on one set of values for α G , β and v of the range provided in [17];

Parameters
The authors use a value of j = 200m/s for the drift velocity of the electrons. We calculate the corresponding current density by taking the electron density to be n 10 29 m −3 and the electron charge e = 1.6 × 10 −19 As.

Conclusion and Outlook
In summary, we extended the phenomenological theory of a uniaxial antiferromagnet with DM interactions to incorporate the current-induced spin-orbit torques together with the already studied spin-transfer torques. We used this theory to analyze the translational skyrmion motion in the presence of a time-independent external magnetic field and a DC electric current. We find that the magnetic field merely modifies the shape of the antiferromagnetic skyrmion and does not contribute towards the skyrmion motion. Further, our results show that the skyrmion moves in a straight line, along the direction of the applied electric current. This agrees with the numerical results of [16,17], which were obtained for skyrmions in the absence of a magnetic field. Depending on the choice of parameters, we find skyrmion velocities that are in the range of 1 − 1000 ms −1 , in agreement with the numerical results of [16,17].
Numerical simulations need to be performed in the presence of an external magnetic field to verify our analytical results. Another direction to proceed would be to allow for a variable skyrmion radius within the collective coordinate approach. The latter could uncover additional interesting physics.