Phenomenological description of the nonlocal magnetization relaxation in magnonics, spintronics, and domain-wall dynamics

Phenomenological description of the nonlocal magnetization relaxation in magnonics, spintronics, and domain-wall dynamics


A phenomenological equation called Landau-Lifshitz-Baryakhtar (LLBar) equation, which could be viewed as the combination of Landau-Lifshitz (LL) equation and an extra “exchange damping” term, was derived by Baryakhtar using Onsager’s relations. We interpret the origin of this “exchange damping” as nonlocal damping by linking it to the spin current pumping. The LLBar equation is investigated numerically and analytically for the spin wave decay and domain wall motion. Our results show that the lifetime and propagation length of short-wavelength magnons in the presence of nonlocal damping could be much smaller than those given by LL equation. Furthermore, we find that both the domain wall mobility and the Walker breakdown field are strongly influenced by the nonlocal damping.

75.78.Cd, 76.50.+g, 75.60.Ch

I Introduction

The genuine complexity of magnetic and spintronic phenomena occurring in magnetic samples and devices imposes both fundamental and technical limits on the applicability of quantum-mechanical and atomistic theories to their modeling. To a certain degree, this challenge can be circumvented by exploiting phenomenological theories based on the continuous medium approximation. The theories operate with the magnetization (i.e. the magnetic moment density) and the effective magnetic field as generalized coordinates and forces respectively Landau and Lifshitz (1935); Akhiezer et al. (1968). The effective magnetic field is defined in terms of various magnetic material parameters, which are determined by fitting theoretical results to experimental data, and at least in principle, can be calculated using the quantum-mechanical or atomistic methods. However, solving the phenomenological models analytically is still a formidable task in the majority of practically important cases. The difficulty is primarily due to the presence of the long range magneto-dipole interaction and associated non-uniformity of the ground state configurations of both the magnetization and effective magnetic field. Hence, the phenomenological models are solved instead numerically, using either finite-difference or finite-element methods realized in a number of micromagnetic solvers Donahue and Porter (1999); Scholz et al. (2003); Fischbacher et al. (2007); Berkov and Gorn (2008); Vansteenkiste et al. (2014).

Traditionally, the software for such numerical micromagnetic simulations of magnetization dynamics is based on solving the Landau-Lifshitz equation Landau and Lifshitz (1935) with a transverse magnetic relaxation term, either in the original (Landau) Landau and Lifshitz (1935) or ”Gilbert” Gilbert (2004) form. Over time, dictated by the experimental and technological needs, the solvers have been modified to include finite temperature effects García-Palacios and Lázaro (1998) and additional contributions to the magnetic energy (and therefore to effective magnetic field) Shu et al. (2004). The recent advances in spintronics and magnonics have led to the implementation of various spin transfer torque terms Zhang and Li (2004); Thiaville et al. (2005) and periodic boundary conditions Lebecki et al. (2008); Wang et al. (2010); Krueger and Selke (2013). Furthermore, the progress in experimental investigations of ultrafast magnetization dynamics Kirilyuk et al. (2010) has exposed the need to account for the variation of the length of the magnetization vector in response to excitation by femtosecond optical pulses, leading to inclusion of the longitudinal relaxation of the magnetization within the formalism of numerical micromagnetics Au et al. (2013). Provided that a good agreement between the simulated and measured results is achieved, a microscopic (i.e. quantum-mechanical or atomistic) interpretation of the experiments can then be developed.

The described strategy relies on the functional completeness of the phenomenological model. For instance, a forceful use of incomplete equations to describe phenomena originating from terms missing from the model may result in false predictions and erroneous values of fitted parameters, and eventually in incorrect conclusions. The nature of the magnetic relaxation term and associated damping constants in the Landau-Lifshitz equation is of paramount importance both fundamentally and technically. It is this term that is responsible for establishment of equilibrium both within the magnetic sub-system and with its environment (e.g. electron and phonon sub-systems), following perturbation by magnetic fields, spin currents, and/or optical pulses Kirilyuk et al. (2010). Moreover, it is the same term that will eventually determine the energy efficiency of any emerging nano-magnetic devices, including both those for data storage Terris and Thomson (2005) and manipulation Kruglyak et al. (2010).

In this report, we demonstrate how the phenomenological magnetic relaxation term derived by Baryakhtar to explain the discrepancy between magnetic damping constants obtained from ferromagnetic resonance (FMR) and magnetic domain wall velocity measurements in dielectrics Baryakhtar (1998); Bar’yakhtar et al. (1986); Baryakhtar and Danilevich (2013) can be applied to magnetic metallic samples. We show that the Landau-Lifshitz equation with Baryakhtar relaxation term (Landau-Lifshitz-Baryakhtar or simply LLBar equation) contains the Landau-Lifshitz-Gilbert (LLG) equation as a special case, while also naturally including the contribution from the nonlocal damping in the tensor form of Zhang and Zhang Zhang and Zhang (2009) and De Angeli De Angeli et al. (2009). The effects of the longitudinal relaxation and the anisotropic transverse relaxation on the magnetization dynamics excited by optical and magnetic field pulses, respectively, in continuous films and magnetic elements were discussed e.g. in Ref. [17; 25; 26]. So, here we focus primarily on the manifestations of the Baryakhtar relaxation in problems specific for magnonics Kruglyak et al. (2010) and domain wall dynamics Hillebrands and Thiaville (2006); Weindler et al. (2014). This is achieved by incorporating the LLBar equation within the code of the Object Oriented Micromagnetic Framework (OOMMF) Donahue and Porter (1999), probably the most popular micromagnetic solver currently available, and by comparing the results of simulations with those from simple analytical models. Specifically, we demonstrate that the Baryakhtar relaxation leads to increased damping of short wavelength spin waves and to modification of the domain wall mobility, the latter being also affected by the longitudinal relaxation strength.

The paper is organized as follows. In Sec. II, we review and interpret the Baryakhtar relaxation term. In Sec. III, we calculate and analyze the spin wave decay in a thin magnetic nanowire. In Sec. IV, we simulate the the suppression of standing spin waves in thin film. In Sec. V, we analyze the domain wall motion driven by the external field and compare the relative strength of contributions from the longitudinal and nonlocal damping. We conclude the discussion in Sec. VI.

Ii Basic equations

In the most general case, the LLBar equation can be written as Baryakhtar (1998); Dvornik et al. (2013a)


where is the gyromagnetic ratio and the relaxation term is


Here and in the rest of the paper, the summation is automatically assumed for repeated indices. The two relaxation tensors and describe relativistic and exchange contributions, respectively, as originally introduced in Ref. [21].

To facilitate comparison with the LLB equation as written in Ref. [29], the magnetic interaction energy of the sample is defined as


where is the equilibrium magnitude of the magnetization vector at a given temperature and zero micromagnetic effective field, i.e. the effective field derived from the micromagnetic energy density , as used in standard simulations at constant temperature under condition (i.e. with only the transverse relaxation included). The second term in right-hand side of Eq. (3) describes the energy density induced by the small deviations of the magnetization length from its equilibrium value at the given temperature, i.e., , and is the longitudinal magnetic susceptibility. Therefore, the associated effective magnetic field is


where , is the effective magnetic field associated to . Hereafter we assume that our system is in contact with the heat bath, so that the equilibrium temperature and associated value of and remain constant irrespective of the magnetization dynamics.

In accordance with the standard practice of both micromagnetic simulations and analytical calculations, to solve LLBar equations (1-4), one first needs to the corresponding static equations obtained by setting the time derivatives to zero and thereby to derive the spatial distribution of the magnetization in terms of both its length and direction. We note that, in general (e.g. as in the case of a domain wall), the resulting distribution of the longitudinal effective field and therefore also of the equilibrium magnetization length is nonuniform, so that the length is not generally equal to . With the static solution at hands, the dynamical problem is solved so as to find the temporal evolution of the magnetization length and direction following some sort of a perturbation. Crudely speaking, the effect of the relaxation terms is that, at each moment of time, the magnetization direction relaxes towards the instantaneous direction of the effective magnetic field, while the magnetization length relaxes towards the value prescribed by the instantaneous longitudinal effective magnetic field. The effective field itself varies with time, which makes the problem rather complex. However, this is the same kind of complexity as the one that has always been inherent to micromagnetics. The account of the longitudinal susceptibility within the LLBar equation only brings another degree of freedom (the length of the magnetization) into the discussion. One should note however that the longitudinal susceptibility has a rather small value at low temperature and so its account is only required at temperatures of the order of the Curie temperature.

We neglect throughout the paper any effects due to the anisotropy of relaxation, which could be associated e.g. with the crystalline structure of the magnetic material Baryakhtar (1998); Dvornik et al. (2013a). This approximation is justified for polycrystalline and amorphous soft magnetic metals, as has been confirmed by simulations presented in Ref. [25]. Hence, we represent the relaxation tensors as and where parameters and are the relativistic and exchange relaxation damping constants and is the unit tensor. Then, Eq. (1) is reduced to


We separate the equations describing the dynamics and relaxation of the length and direction of the magnetization vector. Representing the latter as a product of its magnitude and directional unit vector , we can write


We multiply this equation by to obtain,


Then, subtracting the product of equation (7) and from equation (6), we obtain


where . In the rest of the paper, we will use to represent the component of the vector that is perpendicular (transverse) to vector . Note that only the perpendicular component of the torque contributes to . For given temperature, is constant and we can define . In the limiting case of , and thus is recognized as the Gilbert damping constant from the LLG equation. Let us now consider the case of . The corresponding contribution to the relaxation term, which we denote here as , can be written as


where and the quantity has the form of some magnetization current density (magnetization flux).

For the following, it is useful to split the effective field into its perpendicular (relative to ) part (, “perpendicular field”) and parallel part (, “parallel field”), i.e., , and then to consider the associated magnetic fluxes and torques separately. The magnetic flux of and then the contribution of the associated torque onto is


The perpendicular field can be represented as


So, we can write for the magnetization flux associated with the perpendicular field


The right-hand side of Eq. (12) could be regarded as the torque generated by spin current pumping since can be considered as the exchange spin current Tserkovnyak et al. (2009), and then for the associated perpendicular torque , we obtain,


where we have introduced variable . We show that the torque could be written as (see Appendix A for details)


where is a tensor Zhang and Zhang (2009); Fähnle and Zhang (2013),


In the limit of , we assume and therefore obtain


At the same time, Eq. (8) can then be written as




and is the transverse component of the effective field. The first term in Eq. (18) is kept as since . In practice, we use Eq. (17) rather than Eq. (II) for numerical implementation. As shown in Eq. (II) the damping terms contain both the form Hankiewicz et al. (2008); Tserkovnyak et al. (2009) and tensor form  Zhang and Zhang (2009). Hence, we conclude that the exchange damping can be explained as the nonlocal damping, and Eq. (17) is the phenomenological equation to describe the nonlocal damping.

The intrinsic Gilbert damping is generally considered to have the relativistic origin Landau and Lifshitz (1935); Hickey and Moodera (2009). Phenomenologically, the Gilbert damping is local and the damping due to the nonuniform magnetization dynamics being ignored Gilbert (2004). The exchange relaxation term in the LLBar equation describes the nonlocal damping due to the nonuniform effective field. Despite the complexity of various damping mechanisms, the spin current in conducting ferromagnets can be calculated, e.g. using the time-dependent Pauli equation within the s-d model. The spin current is then given by , where is the conductivity Zhang and Zhang (2009), and thus the nonlocal damping of the tensor form can be obtained Zhang and Zhang (2009); Fähnle and Zhang (2013). As we can see from Appendix A, this spin current densities and have the same form, and therefore, we can establish that . The spin current component (see Appendix A) gives the term  Tserkovnyak et al. (2009), and the value of can be therefore interpreted as , where is the conduction electron density, the effective mass and is the transverse spin scattering time Nembach et al. (2013).

It is of interest to compare Eq. (5) with Landau-Lifshitz-Bloch (LLB) equation Atxitia and Chubykalo-Fesenko (2011), which could be written as


where is the reduced magnetization and is the equilibrium magnetization value at temperature . The effective field contains the usual micromagnetic contributions as well as the contribution from the temperature,


where and with Atxitia and Chubykalo-Fesenko (2011). By substituting Eq. (20) into Eq. (19), one arrives at


Meanwhile, if we neglect the term in Eq. (5) and insert the effective field Eq. (3) into Eq. (5), we obtain


As we can see, Eq. (22) is a special case of LLB equation with the assumption that and . However, the LLB equation does not contain the -term (nonlocal damping term) which is the main focus in this work.

Iii Spin wave decay

To perform the micromagnetic simulation for the spin wave decay, we have implemented Eq. (17) as an extension for the finite difference micromagnetics package OOMMF. A new variable for the exchange damping is introduced with , where is a coefficient to scale to the same order as . In practice, was chosen to be .

The simulation geometry has dimensions , and , and the cell size is . The magnetization aligns along the direction for the equilibrium state and the parameters are typical of Permalloy: the exchange constant , the saturation magnetization and the Gilbert damping damping coefficient . The spin waves are excited locally in the region , and to prevent the spin wave reflection the damping coefficient is increased linearly Han et al. (2009) from 0.01 at to 0.5 at .

Figure 1: The spin wave amplitude decay along the rod, for a spin wave was excited locally by applying a microwave of frequency and amplitude in the region . The data were fitted using Eq. (23) with and .

Figure 1 illustrates the spin wave amplitude decay along the rod. The component of magnetization unit vector data for nm were fitted using (23) to extract the wave vector and the decay constant , and good agreement is observed due to the effective absence of spin wave reflection. We use data after having computed the time development of the magnetization for 4 ns to reach a steady state. The injected spin wave energy is absorbed efficiently enough within the right 200 nm of the rod due to the increased damping.

To analyze the simulation data, we exploit the uniform plane wave assumption with its exponential amplitude decay due to energy dissipation, i.e. magnetization with the form , where is the characteristic parameter of the spin wave damping. For a small amplitude spin wave propagation we have Seo et al. (2009)


where , and the effective field of the long rod can be expressed as


where the ‘easy axis’ anisotropy field originates from the demagnetizing field, and the constant measures the strength of the exchange field,


To test the spin wave decay for this system, a sinusoidal field was applied to the rod in the region to generate spin waves.

Figure 2: The spin wave decay constant–wave vector product as a function of the frequency for different values. The slateblue line was drawn using Eq. (29) for the case .

Figure 2 shows the product of spin wave decay constant and wave vector as a function of the frequency. The dependence is linear for the case, which is agreement with the zero adiabatic spin torque case Seo et al. (2009). The addition of a nonzero term leads to a nonlinear relation, and the amplitude of the spin wave decay constant that is significantly larger than that given by the linear dependence. We also performed the simulation for case by using Eq. (5) which shows that term is the leading factor for this nonlinearity (the relative error is less than 1% for ). To analyze the nonlinear dependence, we introduce the complex wave vector ,


By linearizing Eq. (17) and setting the determinant of the matrix to zero we obtained (see Appendix B for details):


where and . The second term of the Eq. (27) is expected to be equal to zero, i.e., . There are two scenarios to consider: First is the case, could be extracted by taking the imaginary part of at Eq. (26):


The linear dependence of as a function of frequency matches the data plotted in Fig. 2. For the case, solving Eq. (27) yields in the linear with respect to the damping constants approximation,


where the dispersion relation for the rod is . Eq. (29) shows there is an extra term associated with the exchange damping term besides the linear dependence between and . The slateblue line in Fig. 2 is plotted using Eq. (29) with and , which shows a good approximation for the simulation data. Besides, this exchange damping could be important in determining the nonadiabatic spin torque. We could establish the value of using the existing experimental data, such as the transverse spin current data Nembach et al. (2013) gives which hints the lifetime and propagation length of short-wavelength magnons could be much shorter than those given by the LLG equation Baryakhtar (2014).

Iv Suppression of standing spin waves

In the presence of nonlocal damping, the high frequency standing spin waves in the thin films are suppressed Baryakhtar (2014). If the magnetization at the surfaces are pinned, the spin wave resonance can be excited by a uniform alternating magnetic field Kittel (1958). With given out-of-plane external field in the -direction, the frequencies of the excited spin waves of the film are given by Stancil and Prabhakar (2009),


where is the film thickness, , and . The excited spin wave modes are labeled by the integer , and the odd has a nonvanishing interaction with the given uniform alternating magnetic field Kittel (1958).

Figure 3: Imaginary parts of the dynamical susceptibility of the film for (a) thickness nm, and (b) thickness nm.

To reduce the simulation time, we consider a system with cross-sectional area nm in -plane and apply the two-dimensional periodic boundary conditions Wang et al. (2010) to the system. We use the Permalloy as the simulation material with external field A/m and the cell size is nm. Instead of applying microwaves to the system, we calculate the magnetic absorption spectrum of the film by applying a sinc-function field pulse to the system Dvornik et al. (2013b). With the collected average magnetization data, the dynamic susceptibility is computed using Fourier transformation Liu et al. (2008). For example, the component is computed using when the pulse is parallel to the -axis.

Figure 3(a) shows the imaginary part of the dynamic susceptibility for a film with nm. As we can see, the spin wave of modes are excited, and the influence of the “exchange damping” is small. However, the presence of the “exchange damping” suppresses the spin wave excitation ( mode) significantly for the film with thickness nm, as shown in Fig. 3(b). The reason is because the damping of the standing spin waves is proportional to in the presence of exchange damping Baryakhtar (2014).

V Domain wall motion

We implemented Eq. (5) in a finite element based micromagnetic framework to study the effect of parallel relaxation process on domain wall motion. The simulated system for the domain wall motion is a one-dimensional (1D) mesh with length 20000 nm and discretization size , a head-to-head domain wall is initialized with its center near nm. In this section, the demagnetizing fields are simplified as and the demagnetizing factors are chosen to be , and , respectively. The domain wall moves under the applied field for ns and the domain wall velocities at different external field strengths are computed. Figure 4 shows the simulation results of domain wall motion under external fields for different susceptibilities without consideration of exchange damping, i.e., . For Nickel and Permalloy, the longitudinal susceptibility is around at room temperature and increases with the temperature up until the Curie point Atxitia and Chubykalo-Fesenko (2011). We find that the longitudinal susceptibilities have no influence on the maximum velocity but change the Walker breakdown field significantly. The domain wall velocity in the limit is almost the same with the case of , which could be explained by the relation that the difference proportional to the ratio of in Eq. (53).

Figure 4: Simulations results of domain wall velocities for various susceptibilities. The parameters used are: , , and . The vertical dash lines are the breakdown fields computed using Eq. (53).

To investigate the effect of longitudinal magnetic susceptibility and exchange relaxation damping on the domain wall motion, we use the remainder of this section for analytical studies. We start from the constant saturation magnetization of one-dimensional domain wall model, such as the 1D head-to-head wall Nakatani et al. (2005). The static 1D domain wall profile can be expressed as


where is the perpendicular component of the unit magnetization vector, is the wall width parameter and is the position of the domain wall center.

We consider the case that the system is characterized by two anisotropies, easy uniaxial anisotropy and hard plane anisotropy , which originate from demagnetization. The aim is to analyze the impact of the longitudinal magnetic susceptibility under the 1D domain wall model, the demagnetization energy density could be written as


where and . In the limit case case, the effective anisotropy energy density can be rewritten as


where is used and is the ratio between hard plane anisotropy strength and easy uniaxial anisotropy strength.

The dynamics of the domain wall with 1D profile can be described using 3 parameters O’Brien et al. (2012): the domain width , the domain wall position and the domain wall tilt angle . In this domain wall model, one can assume that is only a function of time. Thus, the magnetization profile for the head-to-head domain wall is given by


Using the magnetization unit vector to calculate the exchange energy is a good approximation for the case , thus, the total energy density can be rewritten as




Within the 1D domain wall profile, , the longitudinal component of the effective field is obtained:


where is defined as


As we can see, is a function of the tilt angle and the domain wall width . At the static state, should equal zero, i.e., , which gives


Eq. (39) shows that the difference between magnetization length and reaches its maximum at the center of the domain wall due to the effect of the exchange field, which also peaks in the centre of the domain wall. According to Eq. (39), we can estimate that the magnetization length difference is for case. Figure 5 shows the magnetization length differences of a 1D domain wall for various , it can be seen that this approximation for agrees very well with the simulation results.

Figure 5: Simulation results of the magnetization length difference for a 1D domain wall located at nm with and J/m. The demagnetizing factors are selected to be and .

In the dynamic case, is not equal to zero. If we wrote Eq. (37) as , we can find that the nontrivial term that contributes to is


As an approximation for , we expect  Sobolev et al. (1997), which gives


In this approximation, we have ignored the terms containing and thus the amplitude of is influenced by the domain wall velocity only. We employ the Lagrangian equation combined with dissipation function to compute the domain wall dynamics Hillebrands and Thiaville (2006). The Lagrange equations are


where refers to , and . The dissipation function is defined by where


is the dissipation density function.

v.1 Parallel relaxation

We neglect the exchange damping term with assumption that and arrive at


where and are the perpendicular and parallel components of the effective field. If we also assume that , can be approximated by Eq. (11),


Substituting Eq. (41) and Eq. (45) into Eq. (44) and integrating over space, we obtain


where we have ignored the term. This term leads to the optimal domain wall width Hillebrands and Thiaville (2006):


and for the optimal domain wall width reduces to . In what follows, the domain wall width parameter is approximated by the optimal wall width. The parameter is then given by


and it is straightforward to find its minimum , which corresponds to .

The introduced paramter in Eq. (46) is given by and its value is determined by the ratio of and , which could be although we assume . Following the treatment of Ref. [27], the integrated Lagrangian action is given by


where is the Berry phase term Shibata et al. (2011), is a result of the varying magnetization that introduced a pinning potential. However, the potential is fairly small and therefore is negligible since . By substituting Eq. (49) and Eq. (44) into Eq. (42),


where . The domain wall dynamics is governed by Eq. (50), by eliminating we obtain an equation about ,


where is the Walker breakdown field. From Eq. (51) we can find that the critical value of is approximately equal to if , which leads to the maximum value of to be . There exists an equilibrium state such that if ,


which means the Walker breakdown field for case is increased to , i.e.,


where is the maximum vlaue of . For this steady-state wall motion, the domain wall velocity is




Therefore, in the limit case , and the domain wall mobility is given by


where is the minimum value of . In Fig. 4 the corresponding Walker breakdown fields are plotted in vertical dash lines, which gives a good approximation for and cases. The simulation results show that the Walker breakdown field could be changed significantly if the longitudinal susceptibility is comparable to the damping constant.

v.2 Nonlocal damping

Figure 6: Simulations results of domain wall velocities for the limit case that with various exchange dampings. The parameters used are: , and . The vertical dash lines are the breakdown fields computed with Eq. (60).

In this part we consider the domain wall motion influenced by exchange damping for the case that . The dissipation density function (43) thus becomes


where and are the two components of the effective field, and is computed using Eq. (45). After calculation we obtain


We take the same Lagrangian action (49) for and arrive at


Similarly, the corresponding Walker breakdown field changes to


where . The domain wall mobility is given by


As we can see, the nonlocal damping term influences the domain wall motion as well, and we can establish that . Therefore, for the scenarios that , the contributions from the Gilbert and nonlocal damping are of the order of magnitude for both the domain wall mobility and Walker breakdown field.

Figure 6 shows the domain wall velocities for domain wall motion driven by external fields in the limiting case of . The simulation results are based on a one-dimensional mesh with length nm with cell size of nm. The damping is set to and the demagnetizing factors are chosen to be , and . As predicted by Eq. (60), the nonlocal damping leads to an increment of the Walker breakdown field, and Eq. (60) fits the simulation results very well.

Vi Summary

We explain the “exchange damping” in the Landau-Lifshitz-Baryakhtar (LLBar) equation as nonlocal damping by linking it to the spin current pumping, and therefore the LLBar (17) can be considered as a phenomenological equation to describe the nonlocal damping. In the presence of nonlocal damping, the lifetime and propagation length of short-wavelength magnons could be much shorter than those given by the LLG equation. Our simulation results show that the spin wave amplitude decays much faster in the presence of nonlocal damping when spin waves propagate along a single rod. The analytical result shows that there is extra nonlinear dependence scaling with between (the product of spin wave decay constant and wave vector ) and frequency due to the nonlocal damping. Using the micromagnetic simulation based on the LLBar equation, we show that the difference between magnetization length and reaches its maximum at the center of the domain wall. For the cases that where is the longitudinal magnetic susceptibility and is the Gilbert damping, the Walker breakdown field will increase significantly. By using a 1D domain wall model, we also show that both the domain wall mobility and the Walker breakdown field are strongly influenced by the nonlocal damping as well.

We acknowledge the financial support from EPSRC’s DTC grant EP/G03690X/1. W.W. thanks the China Scholarship Council for financial assistance. The research leading to these results has received funding from the European Community’s Seventh Framework Programme (FP7/2007-2013) under Grant Agreement n247556 (NoWaPhen) and from the European Union’s Horizon 2020 research and innovation programme under the Marie Skłodowska-Curie grant agreement No 644348 (MagIC).

Appendix A Derivation of equation (Ii)

We split the perpendicular spin current into two components,


where we write as ,


The torque generated by spin current is given by , i.e.,


where we have used the identities and . Meanwhile, the corresponding torque can be computed by , which gives


Note that can be changed into the tensor form,




Therefore, we obtain for ,


where is a tensor,


Appendix B Derivation of equation (27)

We introduce a new variable to represent the second term in the (23), i.e., , so we have


where . Considering the fact and neglect the high order term , one obtains and thus




Substituting the above equations into (17), we have