Unveiling domain wall dynamics of ferrimagnets in thermal magnon currents: Competition of angular momentum transfer and entropic torque

Control of magnetic domain wall motion holds promise for efficient manipulation and transfer of magnetically stored information. Thermal magnon currents, generated by temperature gradients, can be used to move magnetic textures, from domain walls to magnetic vortices and skyrmions. In the past several years, theoretical studies have focused on ferroand antiferromagnetic spin structures, where domain walls always move toward the hotter end of the thermal gradient. Here we perform numerical studies using atomistic spin dynamics simulations and complementary analytical calculations to derive an equation of motion for the domain wall velocity in ferrimagnets. We demonstrate that in ferrimagnets, domain wall motion under thermal magnon currents shows a much richer dynamics. Below the Walker breakdown, we find that the temperature gradient always pulls the domain wall toward the hot end by minimizing its free energy, in agreement with the observations for ferroand antiferromagnets in the same regime. Above Walker breakdown, the ferrimagnetic domain wall can show the opposite, counterintuitive behavior of moving toward the cold end. We show that in this case, the motion to the hotter or the colder ends is driven by angular momentum transfer and therefore strongly related to the angular momentum compensation temperature, a unique property of ferrimagnets where the intrinsic angular momentum of the ferrimagnet is zero while the sublattice angular momentum remains finite. In particular, we find that below the compensation temperature the wall moves toward the cold end, whereas above it toward the hot end. Moreover, we find that for ferrimagnets, there is a torque compensation temperature at which the domain wall dynamics shows similar characteristics to antiferromagnets, that is, quasi-inertia-free motion and the absence of Walker breakdown. This finding opens the door for fast control of magnetic domains as given by the antiferromagnetic character while conserving the advantage of ferromagnets in terms of measuring and control by conventional means such as magnetic fields.


I. INTRODUCTION
Fundamental interest in the understanding of the interaction of thermal stimuli and magnetic domains has been propelled by its potential to impact recording and processing technologies for magnetically stored information [1,2]. Control of magnetic states by thermally generated stimuli is hence a growing field of research. Prominent examples include the fields of spin caloritronics [3,4], e.g., domain wall (DW) motion by temperature gradients [5][6][7][8][9][10][11][12], and the field of ultrafast spin dynamics, e.g., thermally induced magnetic toggle switching by ultrafast heat load in ferrimagnets (FIs) [13][14][15][16][17][18][19][20][21][22]. Other coherent means of manipulating magnetic textures, such as heat-assisted magnetic recording (HAMR) [23,24] and helicity-dependent all-optical switching (HD-AOS) [25][26][27][28][29][30]-albeit not primarily induced thermally-are facilitated tremendously by additional application of an ultrashort thermal excitation and subsequent demagnetization [31][32][33][34][35][36]. A key ingredient for the theoretical description of the aforementioned magnetothermal effects, especially thermally induced DW motion, lies in the understanding of transport processes for energy and angular momentum. While in metals thermal spin currents are also transported by electrons, in insulators magnons, low-energy magnetic excitations are responsible for the transport of angular momentum via the spin Seebeck effect (SSE) [37]. Notably, thermal magnons can be used to move magnetic textures, such as DWs, vortices, and skyrmions [38][39][40]. In previous works the DW motion of ferromagnets (FMs) and antiferromagnets (AFMs) induced by temperature gradients has been investigated thoroughly [7][8][9][10][11]. For instance, both experimental [5,41] and theoretical [7,8,10,12] studies on FMs have shown that a DW in a temperature gradient moves toward the hotter end of the sample. On a microscopic level, the hot sample region acts as a magnon source. Since ferromagnetic magnons carry spin, angular momentum conservation dictates that a magnon which is transmitted through a DW exerts an adiabatic spin transfer torque (STT) onto the wall. As a consequence, the DW moves in opposite direction to the magnon propagation direction, i.e., toward the source [7,42,43]. Differently to the mechanism based on angular momentum conservation, an alternative explanation based on thermodynamic arguments has been suggested. Since the DW-free energy decreases as the temperature increases [44], the so-called nonadiabatic entropic torque acts on the magnetization, pulling the magnetic texture toward the hotter region of the sample, thereby maximizing the entropy and minimizing the free energy [8,10]. The generality of the latter picture makes it also applicable to DWs in AFMs, in which thermal magnons do on average not carry angular momentum [9,11] but also to more complex systems such as spin-spirals and skyrmions [39].
Domain wall motion by thermal gradients in AFMs offers complementary properties to the motion in FMs. On the one hand, AFM DW motion can be faster due to the almost complete lack of inertia and the missing Walker breakdown, which limits the maximum velocity. On the other hand, a disadvantage of AFM DWs is the difficulty to manipulate, control, and measure by conventional means, such as external magnetic fields. This kind of conventional magnetization control is only possible in a subclass of AFMs, so-called weak ferromagnets (WFMs) such as the rare earth (RE) orthoferrites, for instance, in which the Dzyaloshinskii-Moriya interaction (DMI) induces a small net-magnetic moment, perpendicular to the Néel order parameter [45,46]. So-called pure AFMs, such as NiO, for instance, in which there is no net-magnetization in bulk, require more sophisticated means of excitation [47,48]. FIs can be seen as a generalization of both systems, FMs and AFMs, since one may selectively tune the relevant magnetic properties by modifying for instance the sample temperature or composition [49,50]. This allows for an enhanced control of the ferromagnetic-or antiferromagnetic-like character of the spin dynamics and enables us to potentially exploit the characteristically fast spin-dynamics of an AFM [49][50][51], while at the same time one can easily manipulate them by using magnetic fields and measure it by conventional detection methods such as the magneto-optical Kerr effect (MOKE) or x-ray magnetic circular dichroism (XMCD).
Naturally, the larger parameter space of the FI, which emerges from the (at least) two nonequivalent magnetic sublattices, also implies that its magnetization dynamics becomes more complex to understand, i.e., the properties of thermal magnon currents strongly depend on the underlying microscopic spin structure [52,53]. Thus, DW motion in FI driven by temperature gradients has been scarcely investigated so far [6] and previous works on DW motion in FIs (and synthetic AFMs) were focused on more controllable stimuli, such as electric currents [51,54] and magnetic fields [49].
In this work we study DW dynamics in FIs driven by thermal magnon currents in constant temperature gradients [53,55]. We use an atomistic spin model based on the stochastic Landau-Lifshitz-Gilbert (LLG) equation to simulate ferrimagnetic DWs in a temperature gradient. Our simulation results will be compared to the previously developed theory for DW motion in FMs [56][57][58] based on the collective coordinates approach. Depending on the strength of the thermal gradient and the base temperature, we find similarities in the DW dynamics to both the FM and AFM. For instance, we can find a Walker breakdown as observed for FMs [7,8], but we also find the quasi-inertia-free motion observed in AFMs [9]. However, in addition we find, a feature that is unique to the FI and, as far as we know, has so far been reported neither for the FM nor the AFM: a motion toward the cold sample region in the case of a FI below angular momentum compensation and above Walker breakdown. Using a theoretical model based on linear spin-wave theory (LSWT) we show that this peculiar motion is due to angular momentum transfer and not linear momentum transfer.

A. Atomistic spin model
We model the most simple kind of FI, that is, a twosublattice FI with a rock salt structure (G-type magnetic ordering) as depicted in Fig. 1. Our atomistic spin model is based on an extended Heisenberg Hamiltonian, for normalized magnetic moments S i = μ i /μ i . J i j denotes the isotropic Heisenberg exchange coupling and T is the biaxial on-site anisotropy with easy z axis and hard x axis (0 d y i < d z i = 0.5 meV). We use the following exchange parameters for the interaction between spins located on the same sublattice A/B, J AA = 16 meV, J BB = 0.5 meV, and between spins on different sublattices, J AB = −6 meV. These values are characteristic for ferrimagnetic RE-transition metal (TM) alloys [59], which are testbed materials in the field of ultrafast spin dynamics [13][14][15][16][17] and are receiving increasing attention in the field of spintronics [54].
The time evolution of the spins is computed with the stochastic LLG equation [60] 1 with the effective field H i = −∂H/∂S i + ζ i , containing both the deterministic field from the spin Hamiltonian H, Eq. (1), and the stochastic field ζ i in the form of Gaußian white noise, For the atomic magnetic moments we use μ A = 4μ Bohr and μ B = 5μ Bohr . Here, for simplicity, we assume that the Gilbert damping and gyromagnetic ratios are the same for both sublattices. The gyromagnetic ratios are γ i = 2μ Bohr /h = 1.76 × 10 11 s −1 T −1 , and the Gilbert damping is set to α G = 0.01. By numerical integration of Eq. (2), for a range of temperatures we calculate the thermal average of the sublattice-specific as well as the net angular momentum (Fig. 1). The angular momentum compensation temperature, at which the net angular momentum is zero, is found at T A = 107 K (Fig. 1). Moreover, our numerical calculations allow us to determine the Curie temperature of the system, T C = 616 K, in the range of ferrimagnetic RE-TM alloys [61,62]. We assume a lattice constant of a = 250 pm. Similar models have been already used in the literature to model the spin dynamics of ferrimagnetic systems, the most prominent example being GdFeCo alloys used for ultrafast toggle switching [13][14][15][16][18][19][20] and HD-AOS [25,29,30]. Despite their potential key role on such a switching process, the study of DW motion under thermal gradients of such kinds of materials [6] has gained far less attention.

B. Computation of domain wall dynamics
In our simulations a DW is placed in a constant temperature gradient and the magnetization is relaxed to a base temperature of T 0 . The base temperature determines the remanent angular momentum and thus enables us to tune the magnetic properties of the FI. During this relaxation phase (t < 0) we set α G = 1 which efficiently suppresses any DW dynamics. At t = 0 we set α G to 0.01, which releases the DW instantaneously. The wall coordinates, i.e., angles ν and positions Z ν (ν = A,B), are tracked by fitting the wall profiles and to the simulation data. ν is hereby the wall width and m ν (z) is the saturation magnetization for which we assume a linear correction in z to improve the fitting accuracy, compared to a spatially constant saturation magnetization. Absorbing boundaries in the form of enhanced Gilbert damping are applied in the longitudinal direction, whereas the transverse boundaries are periodic in order to have bulklike properties. Due to the sizable inter-sublattice coupling J AB , the deviation of the DW coordinates of the two sublattices A and B from each other are relatively small, such that for the tracking of the DW coordinates it is not necessary to distinguish the DW variables Z, , for the two sublattices A and B. A simulation setup of a typical DW profile in a temperature gradient is shown in Fig. 2.
We use a comparably large grid cross section of 96 × 192 spins, with a length of 480 spins in the direction of the temperature gradient to reduce thermal fluctuations of the data [63]. To handle the large computational effort of almost 10 7 spins in total with simulation times of several hundreds of picoseconds, we use a previously developed, highly efficient GPU accelerated atomistic spin dynamics simulation routine based on the NVIDIA CUDA C-API [64].

A. Adiabatic and nonadiabatic spin transfer torques
The theory of DW motion driven by spin-polarized electric currents in FMs is well established. Initial works have suggested that thermal magnon currents can be viewed as spin currents whose amplitude is proportional to the temperature gradient [8]. For FIs, the question is to what extent a similar picture holds and how theory has to be modified in order to account for the particular properties of FIs.
In a temperature gradient, the spatial variation of the stochastic noise ζ i (T i ), Eq. (3), in the effective field of the LLG equation (2) can be interpreted as sources of thermal magnons. This thermal magnon current acts on a DW in a similar way as the STT used to describe DW motion under spin-polarized electric currents in micromagnetic models [56,57]. In those models, the LLG equation is augmented by two additional torque terms to take into account the interaction of a spin-polarized electron current on the magnetization. The so-called adiabatic torque T ad = −(u · ∇)m ∝ −(∇T · ∇)m (6) and the nonadiabatic torque The parameter β eff is the dimensionless nonadiabaticity and by definition specifies the ratio between the two STTs. The adiabatic torque can be related to angular momentum conservation [7,42,43]: When a magnon current (or spinpolarized electron current) passes through the wall, their polarization is continuously rotated by 180 • , thereby changing the magnons' angular momentum. To obey angular momentum conservation in the combined domain + magnon system, one domain has to grow in size, i.e., the DW has to movethe direction depending on the relative polarization of the particle current and magnetization with respect to the DW. The adiabatic torque amplitude in this case is simply given by (cf. Ref. [43]) where J ∝ ∂T /∂z (see also Ref. [65] and Appendixes A and B) is the spin current density and l fu = (l A − l B )/2 the angular momentum per unit volume (l ν = μ ν m ν /γ ν ), see Fig. 1. The difficulty here is to find an expression for the spin current J for a FI.
On the other hand, Schlickeiser et al. [8] introduced the concept of nonadiabatic entropic torque due to the spatially varying exchange stiffness ∇A eff = (∂A eff /∂T )∇T in a FM.
Here we adapt their model for the FM to the FI, and we find the following expression for the nonadiabatic entropic torque strength (see Appendix C for the temperature dependence of the exchange stiffness A eff ): One of the main results of the present work is the demonstration of the validity of these relations, Eqs. (8) and (9), for FIs by comprehensive comparison to atomistic spin dynamics simulations of the DW motion under a temperature gradient.

B. Dynamics of the domain wall
The conceptual idea behind Eqs. (8) and (9) is that the dynamics of a FI can be viewed as an effective FM with angular momentum given by l fu -a model which has recently been employed in a similar fashion by Kim et al. [49] for field-driven FI-DW motion. Such a model should be valid for the low wall velocities in thermal gradients; although some deviations might occur close to T A since we do not take into account inertial effects proportional toZ, and¨ [48,66].
The dynamics of a rigid, ferromagnetic DW can then be described by the two collective coordinates that are the wall position Z DW and tilting angle DW [57,58]. We can adapt the corresponding equations of motion and rewrite them for the ferrimagnetic DW, leading tȯ where K ⊥ = K yy − K xx is the in-plane anisotropy; α ⊥ eff is the transverse Gilbert damping parameter, known from magnetic resonance for instance [67,68]; and DW = − sign(l fu ) DW is the signed wall width [69]. The coupled equations (10) and (11) have two kinds of steady-state solution.
Below the Walker breakdown "current" the driving stimulus is insufficient for DW to overcome the potential barrier, scaling with K ⊥ , and hence the wall angle becomes stationary and, consequently, the motion linear. In this case, the DW velocity reads where in the second equality we used the relation [70] commonly derived under magnetic resonance conditions and related to its linewidth in FIs. We note that Eq. (13) is a generalization of the already-known relations for the velocity of the DWs for FMs [8] and AFMs [9]. Importantly, the difference to those cases is that the effective damping in FIs has a nonmonotonous behavior. The relatively simple expression in Eq. (14) holds if the magnetization is not too close to the compensation point, at which an apparent divergence occurs and where more sophisticated damping models would be necessary [67,68]. Thus, very close to this point we expect some deviations in the wall velocity and precession which we derive from the ferromagnetic model here in this section. Although, qualitatively, this singular behavior accelerates the spin dynamics around T A and is the reason FIs with angular momentum compensation points are becoming so relevant for applications and functionalities related to the speed of the spin dynamics [49,51,54]. Additionally, while the micromagnetic exchange stiffness and its temperature dependence for FMs and AFMs is somehow known [71], the specifics of A eff in FIs remains an open problem-especially at elevated temperatures. Accordingly, above the Walker breakdown, the wall angle DW precesses continuously as the repelling force is limited by K ⊥ /|l fu |. The corresponding wall velocity then reads [58] where the positive sign is for the case (β eff − α ⊥ eff )u < 0. In the limit of a vanishing Walker threshold u W (∼K ⊥ ) [Eq. (12)] or high driving currents u u W , this solution converges to Interestingly, the velocity in this limit and in the small damping regime α ⊥ eff 1 can be approximated by the simple relation, V DW ≈ u = Ja 3 /l fu . We note that the small damping condition, α G |l A − l B |/|l A + l B |, holds in a wide range of temperatures when the system temperature is not too close to the compensation temperature.
To evaluate the equations presented in this section, we refer the reader to Appendixes A-C. First, in Appendix A we derive the dispersion relation of the FI from the LLG equation. In Appendix B this dispersion relation is used to calculate the thermal spin current density J, from which we can then derive the adiabatic STT parameter u. Finally, in Appendix C we compute the nonadiabatic STT for this system using the effective exchange stiffness of the FI. Altogether this allows us to compute the effective DW velocity and precession in a self-consistent way from the LLG equation [Eq. (2)] and the spin Hamiltonian [Eq. (1)] without the need of additional parameters. Figure 3 shows the wall velocity V DW we obtained from our atomistic spin dynamics simulations as a function of the temperature gradient ∂T /∂z, i.e., the amplitude u of the driving STT (the calculation of the steady-state wall velocity is described in detail in Appendix D). We study a range of temperature gradients such that we can investigate the dynamics below and above the Walker breakdown. Additionally, we consider four different base temperatures T 0 , ranging from below to above the compensation temperature T A = 107 K, that will allow us to assess the role of the net angular momentum in the dynamics of the DW.

A. Overview of the domain wall dynamics
Below the Walker breakdown (filled symbols) we find that the wall velocity V DW scales linearly with the thermal torque amplitude u as is expected according to Eq. (13). The positive sign of the velocity hereby indicates a motion toward the hotter end as previously predicted for FMs and AFMs [7][8][9][10]. Above the Walker breakdown the situation is vastly different: Here we find that below the compensation point the wall moves toward the cold region (V DW < 0) in the limit u u W and above the compensation point to the hot one (V DW > 0) for any value of u. Nevertheless, the theoretical wall velocity predicted by Eqs. (13) and (15) nicely fits our simulations results for all four temperatures displayed in Fig. 3. Thus, different propagation directions of the DW motion are found for temperatures below and above the angular momentum compensation temperature T A . This implies that around T A (i) the adiabatic torque parameter u(T ) changes sign due to the sign change of J/l fu in Eq. (8) and (ii) the nonadiabaticity β eff (T ) changes sign likewise since the product β eff u, Eq. (9), is strictly positive (for ∂T /∂z > 0).
Another intriguing observation is the apparent increase in the Walker threshold u W for Figs. 3(a)-3(c), where in Fig. 3(c) the threshold was actually too high to be determined by our simulations. This is highly counterintuitive since the critical current in a biaxial magnet is determined by the in-plane anisotropy K ⊥ [see Eq. (11)], which decreases quickly with temperature K ν (T ) ∼ K ν (0)m 3 ν (T ) [72]. Thus, the Walker threshold u W (T ) is expected to decrease monotonically with T , which we observe only at even higher temperatures shown  in Fig. 3(d). However, in the FI, the expression u W ∝ K ⊥ /|l fu | in Eq. (12) increases, since the net-angular momentum l fu (T ) decreases faster than the individual sublattice order parameters m ν (T ) and hence faster than K ⊥ (T ). These insights about the Walker threshold could have impact in the design of ferrimagnetic devices with improved functionalities, as the temperature dependence of the individual order parameters can be readily tuned by material engineering techniques, e.g., by modification of sample composition.

B. Diversity of temperature dependence of the domain wall dynamics
The DW dynamics below the Walker breakdown (u < u W ) behave as one would expect from previous works in FM but with effective parameters accounting for the fact that the FI is composed of two antiferromagnetically coupled sublattices. Thus, it is worthwhile to further investigate the range of validity of this idea. As shown above, the Walker threshold can be controlled via the perpendicular anisotropy parameter, d y A . Therefore, in order to investigate the regime below Walker breakdown (u < u W ) we have to consider systems with biaxial anisotropy. On the other hand, one of the main results of this work is the demonstration that a DW motion toward the cold end of the sample is possible above the Walker breakdown (u u W ). Since for a uniaxial FI (d y ν = 0) the wall motion is always above the Walker breakdown (u W = 0), we can investigate the validity of Eq. (16) for the wall velocity without mixing effects coming from the presence of perpendicular anisotropy. To study the DW dynamics below and above the Walker breakdown more thoroughly, we computed the temperature dependence of the steady-state DW velocity for the biaxial (d y A = 25 μeV) and uniaxial (d y A = 0) FI, shown in Figs. 4(a) and 4(b), for a fixed temperature gradient of ∂T /∂z = 260 K μm −1 .

Domain wall velocity in biaxial systems
As we could already expect from the data in Fig. 3, the wall velocity below the Walker breakdown is only weakly sensitive to the base temperature, as can be seen in Fig. 4(a). The wall velocity decreases only slightly at higher temperatures. However, this effect is already known from previous studies and can be related to the changing equilibrium magnetic properties [11,71]. We can estimate the expected DW velocity by evaluating Eq. (13) (see also Appendix C). For the simulation parameters described in Sec. II A, this yields an STT of β eff u = 8.6 ms −1 and an expected DW velocity of about 95 ms −1 (black cross) and is in good agreement with the simulation data (colored squares), especially when we consider the rough approximations we applied. Note that although the entropic torque (9) alone is proportional to 1/|l fu | and hence expected to diverge at the angular momentum compensation point (white marker color) where l fu → 0, the DW velocity Eq. (13) remains finite since the damping coefficient diverges in the same fashion: α ⊥ eff ∝ 1/|l fu |.

Domain wall velocity in uniaxial systems
More fascinating dynamics can be found for the uniaxial FI, where predominantly the adiabatic STT drives the DW, see Fig. 4(b). Here we can make two clear observations: (i) the direction of motion of the wall changes sign very close to the compensation temperature (white marker color). Below the compensation point the wall moves to the colder sample regions, i.e., copropagates with the magnon current, whereas above the compensation point the regular motion toward the thermal source is obtained. (ii) The absolute value of the wall velocity drastically increases toward the compensation point, which is supported by Eq. (15) for a magnon current density J which depends only weakly on temperature. This assumption, in particular that J does not change sign at T A , is hereby motivated by its derivation from the dispersion relation (see Appendix A), which does not depend on temperature in a qualitative way, as long as T T C [73]. However, to fully understand the DW dynamics for the freely precessing wall, far above the Walker breakdown, we first need a better understanding of the origin of the spin current density J. Until now, it is not even clear what the sign of the spin current density J is in the FI [52]. For a single sublattice FM the magnetic moment of the magnon is given by the reduction of the S z spin component with respect to the saturation value (m z = ±1) and is hence antiparallel to the ground-state magnetization [7,42]. In the FI the netmomentum of the magnon will be determined by the ratio of the two components μ A S z A /γ A and μ B S z B /γ B relative to each other.
The classical spin-wave amplitudes for the uniaxial FI at T = 0 K can be calculated from LSWT, following Refs. [55,74] (see Appendix A). We find that the lowfrequency branch (σ = −1) in the dispersion relation carries momentum parallel to m z A = +1 and the high-frequency . The magnon current density J σ of the two branches (σ = ±1) of the dispersion relation carry opposite angular momentum. We find that for both base temperatures T , the high frequency branch σ = +1 (light blue) dominates over the low-frequency branch σ = −1 (orange), see Appendix A. This leads to a polarization of the net-magnon current density J = J + + J − parallel to the l B sublattice angular momentum (in the source region), i.e., |J + | > |J − |. Due to the change in the magnon polarization when passing through the wall, to satisfy angular-momentum conservation in the combined domain + magnon system, the domain with the net-angular momentum l fu in the "down" direction has to grow in size. Below T A this is the domain on the hot side (left), whereas above T A it is the domain on the cold magnon side (right). Hence, for the same spin current density J, we obtain different DW propagation directions above and below T A . branch (σ = +1) parallel to m z B = −1. The question now is as follows: Which of these two branches dominates the net spin current? Lower frequency implies higher thermal population-in the classical model that is a population according to a Rayleigh-Jeans distribution n 0 k,σ = k B T /hω k,σand at the same time longer lifetimes τ k,σ . However, the high-frequency branch has a much steeper dispersion relation and hence much higher group velocities v k,σ = ∂ω k,σ /∂k and propagation lengths ξ k,σ = |v k,σ |τ k,σ .
To answer this question, we calculate the thermal magnon current density J quantitatively by solving the following kspace integral (see Appendix A and B for the derivation): The dispersion relation ω k,σ and the magnon propagation lengths ξ k,σ can be written in closed form expression, and ϑ k is simply the angle between the k vector and the z direction.
Thus the numerical solution of Eq. (17) poses only minimal computational effort, and we find that the high-frequency branch of the dispersion clearly dominates the net-magnon current density J. Thus, we have the peculiar situation in which below the compensation point T A , the net-magnon current has a polarization parallel to the ground-state angular momentum l fu -a situation opposite to the case of a simple FM-leading to the opposite direction of DW motion, that is, toward the cold sample regions for T < T A (see Fig. 5, left). We can use the spin current Eq. (17) to compute the adiabatic STT quantitatively, yielding u = −15.7 ms −1 at T = 0. Moreover, we can now calculate the nonadiabaticity parameter β eff using our previously determined nonadiabatic torque parameter β eff u = +8.6 ms −1 yielding β eff = −0.55.
In combination with Eq. (16) we can now predict a wall velocity of about V DW (T = 0) = −14.8 ms −1 for the 260 K μm −1 thermal gradient at T = 0, which, as already expected, can be well approximated by V DW ≈ u. This value matches quite well with the results from our atomistic spin dynamics simulation, although we are unable to simulate T = 0 exactly. Now that we have gained some insight on the role of the adiabatic STT on the DW motion at low temperatures, we can address the domain wall velocity in the uniaxial FI at finite temperature. In particular, the DW velocity in Fig. 4(b) presents an apparent asymmetry close to the angular momentum compensation point T A . As we discussed above, the magnon current density J above and below T A should not change drastically as it is derived directly from the dispersion relation, cf. Ref. [73]. Thus, we argue that the temperature dependence of u = Ja 3 /l fu is mostly due to l fu (T ). Consequently, we expect the adiabatic STT to act approximately antisymmetrically on the DW in the vicinity of T A , implying min[V DW ] ≈ − max[V DW ], as indicated in Fig. 5. However, the wall velocity from our numerical simulations clearly does not follow this symmetry-we find min[V DW ] = −54 ms −1 and max[V DW ] = +142 ms −1 .
We can trace back the origin of such asymmetry by considering Eq. (16) for the domain velocity in the temperature regime close to T A . We know that both the adiabatic (8) and nonadiabatic STTs (9) scale via ∼1/|l fu |; thus, the temperature dependence of β eff should be mostly due to the softening of the exchange stiffness A eff (T ), i.e., weak for T T C (cf. Appendix C). This allows us to assume |β eff | ∼ const and only include a temperature dependence in the form of the necessary sign change at T A . The temperature dependence of V DW (T ) in Eq. (16) is thus due to α ⊥ eff (T ) [Eq. (14)] and u(T ) = Ja 3 /l fu (T ) [Eq. (8)]. Using these assumptions we can calculate V DW (T ) from our theoretical model, shown as the black line in Fig. 4(b). This model excellently describes the temperature dependence of the wall velocity over the full temperature range, including its asymmetry. We conclude that below T A , the entropic torque and the angular momentum transfer work against each other, whereas above T A they act in the same direction. Naturally, the angular momentum transfer becomes less important if angular momentum conservation is broken, that is, when α ⊥ eff becomes large in the vicinity of T A . At the same time, the contribution of the nonadiabatic term ∼β eff α ⊥ eff /(1 + α ⊥ 2 eff ) increases. These findings demonstrate that a ferrimagnetic DW can be pushed away from a thermal magnon source by angular momentum transfer-an effect which in FMs and AFMs can be achieved only by a less efficient linear momentum transfer, i.e., magnon reflection [12,[75][76][77].

C. Emergence of torque compensation temperature
Aside from the DW velocity, Figs. 4(c) and 4(d) also show the temperature dependence of the dynamics of both the wall tilting (biaxial) and precession (uniaxial). For the biaxial FI the steady-state tilting angle DW gradually decreases with temperature until it reaches zero at a temperature of about 125 K after which it increases again with opposite sense of rotation. A similar behavior is found for the wall precession DW =˙ DW in the uniaxial FI. The wall precession changes sign at the very same temperature of 125 K, suggesting that this phenomenon is independent of the in-plane anisotropy K ⊥ . In the following we define this point of completely suppressed DW tilting and precession the torque compensation temperature T T . However, unlike in the biaxial FI, in the uniaxial FI there is an additional rapid increase of the wall precession frequency DW at the angular momentum compensation point. Such an increased wall precession close to T A was recently predicted for field-driven DW motion in ferrimagnetic GdFeCo by Kim et al. [49]; however, in their case, the sign change of the wall precession coincides with the angular momentum compensation point. For the case of thermal magnon current-driven DW motion T T differs from T A , implying that field-driven DW motion is fundamentally different from thermally-induced motion.
We can understand the existence of the torque compensation by comparison to the AFM. In the AFM, the symmetry of the nonadiabatic STT can only lead to propagation V DW of the wall along the temperature gradient. A rotation of the wall angle DW on the other hand does not occur, since both sublattices try to rotate in opposite directions, canting the sublattice magnetizations instead of tilting the DW angle [9]. In the FI this is also the case, but unlike in the AFM the torques from the two sublattices will in general not have equal magnitude and are thus not fully compensated.
Another explanation of these results is provided by the coupled equations of motions for the collective coordinates Z DW and DW , Eqs. (10) and (11). Crucial is hereby the role of the nonadiabaticity parameter β eff (T ), shown in Fig. 4(e). In the previous section we already determined a crude estimate of |β eff | ≈ 0.55 by calculating the nonadiabaticity at low temperature from Eq. (9). We can use this number again to solve Eq. (11) for the steady-state precession frequency, shown in Fig. 4(d), where we find excellent agreement with the simulation results.
A more rigorous approach to compute the nonadiabaticity β eff (T ), including its temperature dependence, can be computed from the ratio R = V biaxial DW /V uniaxial DW of the wall velocity V biaxial DW below Walker breakdown, Fig. 4(a), and the one for the freely precessing wall V uniaxial DW Fig. 4(b). By dividing Eqs. (13) and (16) we can then solve for β eff (T ) using R(T ) from the numerical simulations: The magnonic torque u drives the wall precession via DW ∝ (β eff − α ⊥ eff )u. Thus, the DW precession is expected to cease for β eff (T ) = α ⊥ eff (T ), or, in other words, the critical gradient u W , Eq. (12), diverges for β eff (T ) → α ⊥ eff (T ). The intersection point β eff (T ) = α ⊥ eff (T ), shown in Fig. 4(e), is in very good agreement with the torque compensation point shown in Figs. 4(c) and 4(d) which have been determined directly from wall tilting and precession, respectively. Note that by definition this intersection point also marks the temperature at which the wall velocities above and below Walker breakdown, depicted in Figs. 4(a) and 4(b), coincide. For our parameters that is a steady-state velocity at T T = 125 K of about 90 ms −1 in both panels. Furthermore, we are now able to generalize one of our findings, namely that the torque compensation point T T is found above the angular momentum compensation temperature T A : If the adiabatic STT mediated by the thermal magnon current u acts repulsive on the DW, then the nonadiabaticity β eff has to be negative to ensure that the nonadiabatic STT (9), i.e., the product β eff u, remains positive (for ∂T /∂z > 0) [8][9][10]. Thus, the term (β eff − α ⊥ eff )u in Eq. (11) can only be zero for an attractive adiabatic STT, since α ⊥ eff is strictly positive.

D. Domain wall motion in time domain
For discussing the DW motion in the time domain it is helpful to compare our results on the FI's dynamics to the previous works on FMs and AFMs [7][8][9]11,57,58]. For the FM we can do this even in a quantitative manner by simply switching the sign of the intersublattice coupling J AB to get a ferromagnetic exchange between the sublattices A and B. The steady-state wall velocity V DW below the Walker breakdown appears to be more or less unaffected by the sign change of J AB , as can be seen in the top panel of Fig. 6. This is indeed expected from Eq.  This is especially problematic for very low gradients as for instance in Figs. 6(a), 6(d) and 6(g).
The steady state below the Walker breakdown is characterized by a constant tilting angle DW , where the torques of the nonadiabatic STT are balanced by the anisotropy torques, see Eqs. (10) and (11). During the initial rotation of the DW to this angle the velocity increases to its steady-state value, and hence one can interpret it as an inertial mass of the wall [78]. As mentioned, these torques are partially compensated in the FI, greatly reducing the tilting angle and therefore also the effective inertia of the DW. For the same reason the Walker breakdown u W , at which the wall starts to rotate continuously, is shifted to much higher critical gradients. At T = 58 K we find a threshold gradient in the FI of about k B |∂T /∂z| FI W /d A FI y ≈ 1.8 nm −1 compared to k B |∂T /∂z| FM W /d y A,FM ≈ 0.15 nm −1 in the FM. At the torque compensation point the FI resembles an AFM, for which there is no tilting and hence the wall can move quasi-inertia-free, i.e., without a relevant acceleration time [9]. Ultrafast DW acceleration in the FI is not only found at exactly the torque compensation point T T but also slightly below, due to the diverging wall precession DW close to the angular momentum compensation point T A [see Fig. 4(d)]. Thus, even though the wall has to tilt by a finite angle, the steady-state angle is reached on ultrashort timescales of only few picoseconds.
It should be noted, though, that there are other effects in an AFM that can be attributed to a mass of the DW [48,66]. However, these effects are much smaller and proportional to the velocity of the wall, which is here restricted by feasible temperature gradients.

V. CONCLUSION
To summarize our results, we calculated the DW dynamics of a FI in a thermal gradient using both large-scale atomistic spin dynamics simulations based on the stochastic LLG equation and analytical calculations based on LSWT. Our simulation results are in good agreement with our theoretical findings that we derived from LSWT. Whereas in the thoroughly studied ferromagnetic systems the adiabatic and nonadiabatic STT lead qualitatively to the same result [7,8,10]-a motion to the hotter sample region-a ferrimagnetic DW reacts differently to these two kinds of torques. The nonadiabatic torque leads to a consistent motion toward the hotter end, as is the case for the FM and AFM and can be explained by the free-energy minimization via an entropic torque [8,9]. On the other hand, the adiabatic STT can either push or pull the ferrimagnetic DW away from or toward the spin-wave source, depending on whether the temperature is below or above the angular momentum compensation point. In the FI the copropagation of the DW with the magnon current at low temperature is not due to linear momentum transfer resulting from magnon reflection [12,[75][76][77], but due to the angular momentum transfer from the transmitted magnons. Moreover, the FI shows another distinct characteristic point, besides the angular momentum and magnetic compensation points, that is a torque compensation point at which we find a reversal of the DW rotation. Consequently, at the torque compensation point the Walker breakdown is strongly suppressed, which suggests that high DW velocities and ultrafast DW acceleration should be achievable at this point.
Finally, we want to mention that first experimental evidence on copropagation of a DW with a thermal magnon current, induced by ultrashort laser pulses, has been reported recently by Shokr et al. [6]. In their work it is reported that DWs in a ferrimagnetic GdFeCo alloy will move away from the laser spot center, i.e., against the thermal gradient and toward the cold region, corroborating our findings for the DW motion above the Walker breakdown.
Note added in proof. Recently, similar equations of motion, as those derived in Sec. III B, have also been proposed by Okuno et al. [79] for DW motion in ferrimagnetic GdFeCo alloy, driven by spin-polarized electrical currents.

ACKNOWLEDGMENTS
The authors thank Philipp Graus and Johannes Boneberg for fruitful discussions. This work was financially supported by the Deutsche Forschungsgemeinschaft through the SFB 767 "Controlled Nanosystems." At the FU Berlin support by the Deutsche Forschungsgemeinschaft through SFB/TRR 227 "Ultrafast Spin Dynamics," Project A08 is gratefully acknowledged. Furthermore, U.R. acknowledges funding from the Deutsche Forschungsgemeinschaft via the Project No. RI 2891/2-1.

APPENDIX A: DISPERSION RELATION OF THE ROCKSALT-TYPE FERRIMAGNET
We start the derivation of the spin-wave dispersion by in- . This ansatz implies an x − y symmetry and thus vanishing in-plane anisotropy, K ⊥ = K yy − K xx = 0, in order to avoid dealing with squeezed magnon states [80]. We assume a groundstate magnetization of m z A = +1 and m z B = −1. Following Refs. [55,74], one can deduce the linearized LLG equation in k space in analogy to the FM and AFM cases: where the frequency matrix on the right-hand side is given by The matrix elements of Eq. (A2) are The structure factors C (n) k are related to the neighbor positions of shell n and can be expressed as 2a ν are hereby the lattice constants of the face-centered orthorhombic unit cell (see Fig. 1), and in the following we assume a ν = a for simplicity. The solution of Eq. (A1) is given by the eigenvalues of Eq. (A2) and can be computed in closed form for the 2 × 2 matrix: We can further simplify this expression by assuming α G 1, leading to where the frequencies ω k ± = { k ± } and damping rates λ k ± = { k ± } are the imaginary and real parts of the complex eigenvalues k ± , respectively. The frequency˜ k simply reads Note that for k = 0 Eqs. (A10) and (A11) coincide with the results of Kamra et al. [[68], Eq. (16) and (17)] for the magnetic resonance mode of a FI.
Next we want to derive the "amplitude" of the magnon. Although the absolute value of a magnon is not well defined in our semiclassical picture, it is sufficient to compute the relative amplitudes between the A and B sublattices, as the absolute values of the amplitudes will cancel for a thermal magnon distribution. These relative amplitudes are related to the (non-normalized) eigenvectors to the eigenvalues k ± of Eq. (A2). The classical equivalent to the magnon amplitude μ k ± μ A,B follows from simple geometrical considerations as or where S is a scaling parameter which quantifies the classical spin-wave amplitudes. The sublattice-resolved magnon amplitudes as defined inside the brackets of Eq. (A15) are shown in Fig. 7. One can clearly see that the sign of μ k σ does not depend on k but only on σ , since for a given branch σ , one sublattice is always excited much more strongly. In fact, apart from the modes close to the point, the magnon amplitude can be approximated by an excitation of only one of the two sublattices: For the low-frequency branch (orange) it is the B sublattice (top), whereas for the high-frequency branch (blue) it is the A sublattice.

Temperature step
First we suppose a system of an extended nanostrip with a temperature profile in the form of a step function T (z) = T 0 + T (−z). The system is then isotropic along the x and y directions and we only expect a net-spin current propagating along the z direction. Since the thermal magnon occupation in the classical limit follows a Rayleigh-Jeans distribution, n 0 k,σ = k B T /hω k,σ , one can drop the base temperature T 0 as long as it is low enough that it does not affect the effective magnetic parameters, i.e., as long as the dispersion relation (A10) and (A11), shown in Fig. 8 is still valid.
Note that unlike previous works that studied the action of spin waves on DWs [43,75], an effective one-dimensional model is not sufficient here, since thermally excited magnons with off-axis wave vector k = kẑ are relevant and due to the large grid cross section included in the numerical simulations. The macroscopic spin current density follows from integrating over all thermally excited modes in the BZ. The two sublattices of the checkerboard AFM (Fig. 1) are two fcc lattices with a magnetic lattice constant of 2a, respectively. Thus, the BZ is a truncated octahedron with q X = π/a [81]. For a DW at a position z > 0 away from the temperature step, we can restrict the k-space integral to the half space v z k > 0 to only include forward-propagating magnons: Each mode k, σ contributes with to the net current J. Here ∂ω k,σ /∂k = v k,σ is the group velocity of the mode, n k,σ is the magnon occupation number at the source, and the exponential factor accounts for the absorption of the current with propagation length ξ k,σ = |v k,σ |τ k,σ . The factor of two in the exponential accounts for the conversion of the spin-wave amplitude to the magnon number, proportional to the squared amplitude. ϑ k denotes the angle between k and the z direction and is needed to compute the actual propagated distance r k = z/ cos ϑ k .
Ignoring depletion effects at the interface, i.e., at the magnon source, we can assume that the magnon occupation at the source is given by the thermal population n k,σ ≈ n 0 k,σ and hence we get and we arrive at our preliminary result for the magnon current density due to a temperature step We should note that we assumed a fixed magnon amplitude ofh in the derivation, which for the quantum mechanical case is a reasonable assumption for the FI [80] but seems arbitrary for the classical case [see Eq. (A14)]. This is, however, no longer relevant, since the magnon amplitude eventually cancels when we put in the thermal occupation n 0 k,σ ∝ 1/hω k,σ . In Fig. 9 we compare the DW velocity calculated according to Eq. (B5) with our numerical simulations results. Since the base temperature is set to zero, the system is below the compensation temperature. The DW velocity is plotted as a function of distance z from a 1-meV temperature step. We find excellent quantitative agreement between numerical simulations and the LSWT. Furthermore, as for the thermal gradients, we also find the motion of the wall to be away from the magnon source.

Temperature gradient
The solution of the previous Appendix B 1 is easily applicable to temperature gradients by simply summing up over several temperature steps dT (z) = ∂T /∂z dz. Once more we will use the fact that the Rayleigh-Jeans distribution is linear in T (z). Therefore, in a constant temperature gradient, we have the same amount of magnons flowing from the right to the left (carrying spin −σh) as we have "magnon-holes" flowing from right to left (carrying spin +σh). This means we can again restrict the k-space integral over half of the BZ and multiply the result by 2, such that the final result for the spin current is the sum of equal contributions of magnon current and "magnon-hole" current, We obtain the final result Eq. (17) by performing the z integration, FIG. 10. Comparison of the semiclassical spin current emitted from a k B T = 1 meV temperature step, derived via the Rayleigh-Jeans distribution (RJ; red line) and the quantum statistical derivation based on the Bose-Einstein distribution (BE) for different base temperatures T 0 . At very low temperature the spin current is dominated by the low-frequency branch, since the high-frequency magnons are still frozen out. Thus the sign of the net-spin current is reversed in the vicinity of the source, indicated by the dashed line.

Quantum effects
One can further compute the spin current in the quantized form by replacing the Rayleigh-Jeans distribution n 0 k,σ with the Bose-Einstein distribution n BE k,σ . For simplicity we restrict this discussion to the case of a temperature step as the findings are expected to be qualitatively similar for the thermal gradient.
In quantum statistics, the magnon occupation number is no longer linear in the temperature, and hence the base temperature will be relevant. The resulting spin current (B2) emitted from our temperature step should thus be proportional to (B8) Figure 10 shows the DW velocity corresponding to a thermal spin current calculated with the correct quantum statistics for a set of base temperatures k B T 0 (the Curie temperature for the given exchange constants is T C = 616 K). At very low temperature only the lowest magnon energies will be occupied, i.e., the low-frequency branch of the dispersion (Fig. 8) will dominate the spin transport, despite the low propagation length. Moreover, it is implied that real systems can exhibit a sign change of the net-spin current J at very low temperature. Thus, for u W → 0, the DW velocity V DW = Ja 3 /l fu [Eqs. (8) and (15)] changes sign not only at the compensation point where l fu changes sign but also a second time when the temperature is sufficiently high to populate the long-range, high-frequency magnons of the upper branch-the ones that carry negative angular momentum (see Fig. 8). The overall magnon current and magnon accumulation at low temperature is greatly reduced with respect to the classical case. However, for higher temperatures and in particular at room temperature, we qualitatively retain the semiclassical magnon current derived with the Rayleigh-Jeans distribution. From this we can conclude that the semiclassical treatment is sufficiently accurate for describing most experiments, which are usually carried out near room temperature with magnetic materials of similar ordering temperature [61].

APPENDIX C: ENTROPIC TORQUE IN THE FERRIMAGNET
The entropic or magnetothermal torque in the FI can be defined in analogy to the FM case [8]. The effective exchange stiffness for our cubical FI is composed of the three exchange contributions A eff = A AA + A AB + A BB . In the molecular field approximation, for a magnetic texture along the (001) direction, these can be written as [71] The different numerical factors come from the symmetry of the shells which is fcc for the ferromagnetic exchanges A AA and A BB (eight neighbors with z = ±1), and simple cubic for the antiferromagnetic exchange A AB (two neighbors with z = ±1), see Fig. 1. Their temperature dependence is hereby assumed to be well approximated by the meanfield expressions A i j /A i j (0) = m i m j . It should be noted that although we chose exchange parameters with J AA J BB and thus J BB is not significantly affecting the magnetic ordering (T C for instance), it does add a non-negligible contribution to the entropic torque due to the faster demagnetization of m B compared to m A . The temperature dependence of the equilibrium magnetizations m i (T ) is taken from the data in Fig. 1. We find ∂m A /∂T ≈ −4.87 × 10 −4 K −1 for the strongly coupled sublattice A and ∂m B /∂T ≈ −2.15 × 10 −3 K −1 for the weakly coupled lattice B. The temperature derivatives which we obtained for the three exchange stiffness contributions are summarized in Table I.
The effective magnetocrystalline anisotropy density is K zz eff = (d z A + d z B )/2a 3 = 5.13 × 10 6 Jm −3 . This value is chosen rather high in order to (i) keep the DW width and hence the required computation grid small and (ii) to reduce the characteristic timescale of the DW acceleration which is proportional to the wall width DW [see Eq. (11)]. In our simulations we observe a DW width of about 1.6 to 2.2 nm (depending on temperature) which is in good agreement with the theoretical prediction of DW = A eff /2K zz eff = 1.50 nm.

APPENDIX D: COMPUTATION OF STEADY-STATE DOMAIN WALL DYNAMICS
Due to the different timescales involved in the ferrimagnetic DW dynamics, determining the steady-state velocity, precession, and tilting is challenging. On the one hand, simulation time should be as short as possible, in order to minimize the thermal drift, i.e., the error margins of the temperature, but at the same time one has to assure that the simulation time is sufficiently long for the wall to reach its steady-state motion.
For the data in Fig. 3 the steady-state dynamics were determined as follows: Below Walker breakdown, we simulated a fixed amount of time of 320, 128, 128, and 384 ps for Figs. 3(a)-3(d), respectively. These numbers reflect the acceleration timescales at the different base temperatures. The first 25% of this simulation time was hereby discarded in order to reach the steady-state velocity (and tilting), and the other 75% were used for computing the time average of V DW displayed in Fig. 3. The Walker breakdown was defined by the wall angle tilting by more than 45 • plus a 5 • error margin, in order to account for the diverging wall precession time at exactly the Walker threshold. Above the torque compensation point, Fig. 3(d), or very close to the Walker thresholds, the DW precession is slow and the precession period can be several hundreds of picoseconds. In this case, the steady-state velocity was time averaged over only one 180 • rotation to keep thermal drift as low as possible and ensure a well-defined temperature. For the faster precessing walls in Figs. 3(a) and 3(b), the precession period can be as low as a few tens of picoseconds, and hence we simulated several precession periods to improve the signal-to-noise ratio. In this case, the time average over 256 and 128 ps of simulation time was taken, respectively (rounded down to the next integer number of 180 • rotations).
In Fig. 4(a) and 4(c), the steady-state velocity was determined by fitting an exponential function ∼V DW (1 − e −t/τ ) to the velocity data using a 320-ps simulation time. This procedure was not applicable in Fig. 3, since the corresponding fits would not converge properly, especially for the lowest temperature gradients.
Finally, for the data in Figs. 4(b) and 4(d) we simply took the time average over a comparably short simulation time of 128 ps, due to the lack of inertia in the uniaxial FI.