Anisotropic critical behavior of current-driven skyrmion dynamics in chiral magnets with disorder

The dynamic pinning effects are significant in manipulating skymions in chiral magnetic materials with quenched disorder. Through numerical simulations of the non-stationary current-driven dynamics of skyrmions with the Landau–Lifshitz–Gilbert equation, the critical current, static and dynamic critical exponents of the depenning phase transition are accurately determined for both adiabatic and non-adiabatic spin-transfer torques and with different strengths of disorder, based on the dynamic scaling behavior far from stationary. We find that the threshold current is insensitive to a small non-adiabatic coefficient of the spin-transfer torque, but dramatically reduced for a large one. The critical exponents indicate that the critical dynamic behavior is robust for different spin-transfer torques in the perpendicular component of the Hall motion, while exhibits a weak universality class in the direction of the driving current. The anisotropic behavior around the depinning phase transition provides a quantitative analysis of the drive-dependent skyrmion Hall effect in experiments. Further, the theoretical analysis using the Thiele’s approach is presented, and the critical current and the static exponents support the simulation results.


Introduction
Skyrmions are topological field configurations, originally proposed to describe the resonance states of baryons in nuclear physics [1]. Such topological spin configurations have been observed in the chiral magnetic materials [2][3][4][5][6]. Skyrmions can be stabilized by the Dzyaloshinskii-Moriya (DM) interaction in ultrathin ferromagnets, and have been recently observed in the triangular crystal form over a wide range of temperatures in various materials [5][6][7][8]. In addition, it has been theoretically and experimentally demonstrated that skyrmions can be driven by the current with an ultra-low current density, which is at least five orders of magnitude smaller than that typically applied in experiments on the current-driven domain-wall dynamics [9][10][11][12][13]. The topological protection, nanoscale structure, low energy cost and low Joule heating of skyrmions are very promising characteristics for future applications in advanced electronic and spintronic devices [14,15].
In the magnetic materials, skyrmions can be driven by current via the spin transfer torque (STT) effect which describes the interaction between the electrical current and the local magnetization. Theoretically the STT can be introduced in the Landau-Lifshitz-Gilbert equation (LLG) by adding two additional terms, which represent the adiabatic and non-adiabatic STTs, respectively [16,17]. Skyrmions generally generate a transverse motion with respect to the current, i.e. the so-call skyrmion Hall effect [11,[18][19][20][21][22]. When a current with an adiabatic STT is applied in the absence of defects, there exists no intrinsic pinning of skyrmions, in contrast to the domain wall and helical magnetic structure. However, the quench disorder in the experimental material will induce a pinning force for skyrmions [10,23]. Thus, the collective motion of skyrmions will exhibit a threshold current which separates the pinning state and the sliding state. Recently there have been various theoretical studies on the depinning phase transition of skyrmions [13,18,[24][25][26][27], and the experiments have also detected the phase transition in thin films with quench disorder, including that for the individual skyrmion [10,15,20,23,[28][29][30][31][32][33].
The skyrmion dynamics could be tackled in some respects with the particle model based on a modified Thiele equation. For example, there have been various efforts to understand the depinning phase transition of skyrmions [24,26,34,35]. The static and dynamic phases of skyrmions driven by the magnetic force have been numerically simulated [24]. The collective transport dynamics of skyrmions around the depinning phase transition exhibits different characteristics for the directions perpendicular and parallel to the driving force [35][36][37]. However, the particle model is a phenomenological model, in which a skyrmion is considered as a particle, and the detailed microscopic structure and distortion of the skyrmion are not concerned. Moreover, the particle model can not well describe the critical dynamic behaviors induced by the non-adiabatic STT of the electric current.
The LLG equation is fundamental in the theoretical study of micromagnetics, and may provide a more comprehensive description of the dynamical properties of skyrmions. The pinning phenomenon of an individual skyrmion driven by current has been investigated [12,25,38]. The pinning and depinning states, and the multiplication and annihilation of the current-driven skyrmions induced by quenched disorder have also been simulated based on the LLG equation [39,40]. Recently, experiments reveal that the skyrmions in disordered materials exhibit a drive-dependent skyrmion Hall effect just above the depinning threshold [21,23,28,41], and the simulation study suggests that the drive dependence is accounted for through the anisotropic behavior in the parallel and perpendicular directions of the skyrmion motion [42]. However, there is not any quantitative study of the anisotropic critical behavior on the depeinning phase transition. Moreover, an effective method to reduce the threshold current in experiments is very significant in the application of skyrmions. Disorder generally exits in magnetic materials prepared for experiments, and the insensitivity of the threshold current to the non-adiabatic coefficient and the disorder was suggested [13,26]. But such a critical current is usually dependent of the non-adiabatic coefficient in other magnetic structures, such as domain walls and helical states [13,17,43]. Our conjecture is that the insensitive behavior of skyrmions to the non-adiabatic coefficient possibly exists only in a restricted regime, and a larger non-adiabatic coefficient may effectively reduce the threshold current. Meanwhile, due to the difficulties in the numerical simulation of the LLG equation, the characteristics of the depinning phase transition of skyrmions, such as the critical exponents, the universality class and its dependence on the non-adiabatic coefficient, have not been obtained. The comprehensive understanding of the critical dynamics of skyrmions also provides a guidance to find the materials with an appropriate non-adiabatic coefficient in the efficient manipulation of skyrmions.
Due to critical slowing down, it is very difficult to simulate the stationary state around the depinning phase transition. The nonstationary dynamic approach looks novel and efficient in tackling the dynamic phase transitions, since the measurements are carried out in the short-time regime of the dynamic evolution [35,[44][45][46][47]. In addition, it may avoid the errors induced by the finite time step Δt in the molecular dynamics simulations of the stationary state. In this paper, we present numerical simulations based on the LLG equation for the skyrmions driven by the electric current around the depinning phase transition in a thin film with quench disorder. The critical current and both the static and dynamic critical exponents are accurately determined with the nonstationary dynamic approach, and the dynamic effect of the non-adiabatic STT is investigated. Further, the theoretical analysis of skyrmions around the depinning phase transition is also carried out, and the critical current and critical exponents from the numerical simulations and the theoretical approach are compared. In section 2, the model and the equation of motion are described, and in section 3, the dynamic scaling analysis and the theoretical solution of the current-driven dynamics of skyrmions are formulated. In section 4, numerical simulations and theoretical results around the depinning phase transition are presented.

Equation of motion
We consider the Heisenberg model with the DM interaction for the chiral magnetic materials. The Hamiltonian on the two-dimensional square lattice is given as where m i is a three-dimensional normalized spin at site i, and á ñ ij represents the nearest-neighbor sites. The first term is the ferromagnetic exchange interaction, and a typically value J=1 meV is adopted in this paper. In our notation, J is used as the unit of other parameters. The second term is the DM interaction, and the third one is the Zeeman effect induced by the external magnetic field h which is set in the −z direction, perpendicular to the x-y plane. The competition between the ferromagnetic and the DM interaction results in the helical state for a small h and the ferromagnetic state for a large h, while the skyrmion crystal state (SkX) emerges for an intermediate h.
In order to obtain a stabilized SkX phase, the DM interaction and the external magnetic field are set to D/J=0.4 and h/J=0.04, respectively [13,40]. The easy-axis anisotropy is introduced in the fourth term with impurities, and Λ is a set of random sites, which represents the quenched disorder. The strength of disorder is described by K/J, and the density of the random sites is 5% of the total lattice sites in our simulations [21].
The current-driven spin dynamics at zero temperature is described by the LLG equation with the STT where γ=gμ B /ÿis the gyromagnetic ratio, α is the Gilbert damping constant, and for convenience, m just denotes m i in equation (1). In this paper, the time integration of the discrete LLG equation is done using the fourth-order Runge-Kutta method, and the central difference discrete method is adopted in the numerical simulation. The time step is typically set to Δt=0.05, and a smaller Δt has been also used to confirm the results. The effective magnetic field H eff acting on the spin m is calculated with the Hamiltonian The last two terms in equation (2) describe the STT effect from the coupling between the spin-polarized electric current and the local magnetization. The third term is the adiabatic spin-transfer torque, which arises from the electron moments staying aligned with the magnetization, when the electrons flow through a region with a nonuniform magnetization. u is a current vector directing along the direction of the conduction-electron motion with the amplitude u=pa 3 j/2eM, where p, a, M are the lattice constant, the spin-polarization and the magnitude of the local magnetic moment, respectively, while j is the magnitude of the spin-polarized electric current. In this paper, u is put along the x direction. The fourth term describes the spin-transfer torque of the non-adiabatic effect, and b is called the non-adiabatic coefficient. We use b instead of β for the non-adiabatic coefficient just to avoid the confusion with the critical exponent β introduced below. The non-adiabatic STT is perpendicular to the adiabatic STT and the magnetization. When b is equal to zero, there only exists the adiabatic STT. In general, the value of the non-adiabatic coefficient b is dependent of the materials, and b can be large for some materials, for example, b=0.96(4) in the ultrathin film Co 70 Fe 30 /Pd [48]. The time t and the current density j are in units of ÿ/J;0.66 ps and 2eMJ/pa 2 ÿ;1.2×10 12 A m −2 , respectively, with the lattice constant a=1 nm, the spin-polarization p=0.4 and the magnitude of the local magnetic moment M=1. In our simulations, the simulation box is typically set to be L x ×L y =288×288 with periodic boundary conditions in both directions. Simulations with other smaller and larger sizes of the simulation box have been also performed, to confirm that the finite-size effect has been negligibly small. In order to obtain the initial spin configurations of a perfect skyrmion crystal, we use the replica Monte Carlo simulation to generate a spin configuration at a low temperature to reduce the computational complexity of the LLG equation, and further update it with the LLG equation without the driving current, i.e. u=0, until the sufficient convergence of the spin configuration [13,19]. With such initial configurations at hand, we then apply a steady electric current in the x direction, and start the numerical simulation of the relaxation dynamics with the LLG equation. The maximum time in our simulations is t max =10000, which is sufficient to observe the dynamic scaling behavior of the order parameter at the depinning phase transition, and measure the critical current and critical exponents with reasonable errors. Simulations up to longer times will lead to a higher accuracy, but it is difficult with our computer resources due to the high computational complexity induced by the finite-Δt and finite-size effects. The number of total samples of the quenched disorder for average is about 10000 to ensure the statistical errors of the measured quantities are reasonably small.

Dynamic scaling analysis
Skyrmions are characterized by the topological number defined by x y Under this definition, a single skyrmion results in Q=1, and the total number of skyrmions is given by Q. On a discrete lattice, Q may fluctuate slightly around, but the fluctuation is very small. It has been shown that the swirling structure of a skyrmion induces an emergent magnetic flux [10]. In two-dimension materials, this emergent magnetic flux will generate a magnetic field perpendicular to the film plane. In our simulation, the emergent magnetic field B e =(0, 0, Q) is directed along the z direction. Thus skyrmions drifting with a velocity v d must induce an emergent electric field in complete analogy to the Faradays law, [10,40]. To evaluate the velocity of the moving skyrmion, we measure the emergent electric field [10,40], The average velocity of the skyrmions is , 6 x y e y x e and  á ñ represents the statistical average of all samples.
In the stationary state, the collective motion of skyrmions driven by a current exhibits a dynamic depinning phase transition. The threshold current u c separates the pinned liquid state and the moving amorphous glass state of the skyrmions, with the zero and nonzero average velocities respectively. Assuming the depinning phase transition to be of second order, and based on the renormalization group arguments, therefore, there should exist a non-stationary dynamic scaling form already in the macroscopic short-time regime, after a microscopic time scale t mic [49,50]. Here t mic is the time that the dynamic system needs to sweep away the microscopic details such as those of the initial states and interactions, and it is sufficiently large in the microscopic sense but still small in the macroscopic sense. The macroscopic dynamic scaling behavior emerges only for tt mic . In the critical regime and for a sufficiently large simulation box, the dynamic scaling form of the order parameter, for example, the average velocity of the skyrmions in the perpendicular direction, is described by where τ=(u − u c )/u c is the reduced current, β ⊥ and ν ⊥ are static critical exponents, and z ⊥ is the dynamic critical exponent. At the critical current, τ=0, it leads to a power-law behavior, Therefore, searching for the best power-law behavior of v y (t, τ), one may locate the transition current u c , and then measure the exponent β ⊥ /ν ⊥ z ⊥ . The exponent 1/ν ⊥ z ⊥ can be extracted from the logarithmic derivative of equation (7), Further, the Binder cumulant of the velocity is defined as At the critical current, based on the finite-size scaling analysis and assuming the non-stationary spatial correlation length ξ ⊥ (t) is small, the Binder cumulant should scale as where d=2 is the spatial dimension. In the scaling regime, ξ ⊥ (t) usually grows in a power law,

Theoretical solution
With the Thiele's approach, the LLG equation in equation (2) can be reduced to a relatively simple equation for the collective motion of skyrmions [51]. Assuming that the spin textures of skyrmions are rigid during the current-driven motion, and the skyrmions collectively move with a steady velocity v d , one may write the local magnetization at ( ) Under these assumptions, the LLG equation is reduced to the so-called Thiele equation, which describes the collective motion of skyrmions in the stationary state [13,52] v Due to the special rotation symmetry of the skyrmion, ( ) =  G 0, 0, , and p =  Q 4 . In our simulation, we have computed the total skyrmion number in the presence of disorder, and it just changes about one percent compared to that without disorder, although the skyrmions may be distorted due to the disorder.  is a dissipative force tensor, and its components are given by In the absence of the distortion of the skyrmion,  is typically diagonal, i.e. d = mn mn  , and  is a scalar. In fact, the dissipation tensor will not be perfectly diagonal, if skyrmions are distorted due to the disorder. In most previous works, however, the non-diagonal components of the tensor are neglected in order to simplify the theoretical analysis [13,40,52]. We assume that  remains constant in the stationary state for convenience, and p =  Q 6.976 can be evaluated at h/J=0.04 without quenched disorder. F pin is the pinning force induced by impurities, and its phenomenological form was suggested in [10,52] v where v v | | = d d ,  represents the strength of the pinning force, v pin is a characteristic velocity. In other words, F pin is in the opposite direction of v d . The function f (v d /v pin ) parameterizes the dependence of the pinning force on the velocity [52]. For a small current u, one expects f (v d /v pin )≈1.
Equation (13) Taking into account v d 2 =v x 2 +v y 2 , we can explicitly derive the expression for v x , v y and v d , In [13], similar results have been also obtained except for the magnitude of the velocity, v d . In our work, we will determine the critical current from the solution of v d . Suppose the depinning phase transition is of second order, the magnitude of the velocity of skyrmions tends to zero, i.e. v  0 d , at the critical current u c . Thus we obtain the critical current u c from equation (20), The constant  can be estimated from the critical current u c with the adiabatic STT (b=0) in numerical simulations.

Simulations
In numerical simulations, we first set the non-adiabatic coefficient b=0.00, 0.005, 0.05, 0.5, and 1.0, respectively, with a fixed strength of disorder, K/J=1.0, for a comprehensive understanding of the dynamic effect of the STTs. Then we extend the simulations to other strengths of disorder with the non-adiabatic coefficient b=0.5. We measure the average velocities and their fluctuations for the collective motion of skyrmions in both directions perpendicular and parallel to the driving current, and determine the critical current u c and all the static and dynamic critical exponents. Below, the detailed analysis of the simulations is presented for the non-adiabatic coefficient b=0.5. The results for all values of b and for different strengths of disorder are summarized in tables 1 and 2.
In figure 1, the time evolution of the structure of skyrmions in a magnetic thin film with quenched disorder is plotted at the critical current u c =0.0240. At t=0, the initial configuration is a perfect skyrmion crystal, which is the spontaneous ground state of skyrmions without quenched disorder. We have also considered other initial states of skyrmions, for example, the disordered state, or the glass-like state. The resulting dynamic behaviors are essentially the same, but the correction to scaling is relatively small for the perfect crystal initial state. The driving current is suddenly switched on when the simulation starts. As the time evolves, skyrmions may be distorted due to the quenched disorder, meanwhile the crystal structure gradually deforms and tends to a glass-like state [24]. If the driving current is sufficiently large, in principle, the configuration of a skyrmion could be broken, and skyrmions may be annihilated or generated [40]. In this work, however, the driving current is weak or medium, and the number of skyrmions is conserved. A skyrmion may change its shape, but still keep the topological number very close to Q=1 due to the topological protection.
In figure 2(a), the perpendicular velocity v y (t) versus the time t is plotted for different values of the current u. In the inset, the coefficient R 2 of determination, which measures the quality of the power-law fit, is shown as a function of the current, and the critical current u c =0.0240(1) is determined from the maximum in the time regime tt mic ∼300. From the slope of the curve, one measures the critical exponent β ⊥ /ν ⊥ z ⊥ =0.68(1), according to equation (8). The errors include both the statistical ones and the fluctuations along the time direction. In figure 2(b), the logarithmic derivative ∂ τ ln v(t, τ) at the critical current is displayed for both the perpendicular and parallel directions. With equation (9), the slopes of the curves yield the critical exponents 1/ν ⊥ z ⊥ =0.69(1) and 1/ν P z P =0.59(1) respectively. Compared with the symbols sizes, error bars of the data points are small and invisible in figure 2. Errors of the critical current and critical exponents include the errors of statistical samples and the fluctuations in the time direction.
The spatial correlation length ξ(t) can be computed from the Binder cumulant of the velocities, based on equations (10) and (11), which is shown in figure 3(a). The correlation length in the perpendicular direction exhibits an almost perfect power-law behavior, and the slope of the curve gives the dynamic critical exponent 1/z ⊥ =0.56 (1). However, there exits a strong correction to scaling at early times for the correlation length in the parallel direction. Fitting the curve with a power-law correction ξ x (t)∼  t z 1 (1+c/t), it yields the dynamic exponent 1/z P =0.59 (2).
Theoretically one could also consider to determine the critical current from the velocity parallel to the driving current. Due to the correction to scaling of the correlation length ξ x (t) in the parallel direction, however, a visible deviation from the power-law behavior at the critical current u c =0.0240 is observed. Such a phenomenon is similar to that in the particle model of skyrmions [35]. To circumvent this difficulty, the  (2) 1.82 (4) 2.05 (15) 0.91 (7) 3.45 (17)   correlation length ξ x (t) rather than the time t should be taken as the dynamic scaling variable. Thus the dynamic scaling form of the parallel velocity is rewritten as At the critical current, τ=0, the parallel velocity obeys a power law, v In the inset of figure 3(a), this power-law behavior is displayed at the critical current u c =0.0240. One thus measures the critical exponent β P /ν P =1.27(2) from the slope of the curve.
In figure 3(b), the critical current u c measured from the LLG simulation is plotted versus the non-adiabatic coefficient b. Through the theoretical analysis based on the Thiele's approach to the depinning phase transition in the stationary state, we obtain also the relation between the critical current u c and the non-adiabatic coefficient b in equation (21). For comparison, we estimate the strength of the phenomenological pinning force, i.e. the constant , from the critical current u c for the adiabatic STT (b=0). In figure 3(b), the red solid curve represents the theoretical result of u c based on equation (21), which coincides well with the numerical simulation of the LLG equation. In [13], it is reported that the critical current is insensitive to the non-adiabatic coefficient b. Our result does support this statement in the small b region. But the critical current u c is dramatically reduced for a large b.
In table 1, the critical current and critical exponents of the perpendicular and parallel directions for different values of the non-adiabatic coefficient b are summarized. The critical exponents β ⊥ , ν ⊥ and z ⊥ are robust for different b including b=0, indicating a strong universality class for the collective motion of skyrmions in the perpendicular direction. In other words, the characteristic of the Hall motion is independent of the nonadiabatic coefficient b. For the parallel direction, however, it exhibits a weak universality class. At b=0, the  critical exponents β P and z P are significantly larger than β ⊥ and z ⊥ . As the non-adiabatic coefficient b increases, β P and z P decrease, while ν P is opposite. These results are very different from those of the particle model [35]. In the particle model, the Magnus force is phenomenologically introduced to describe the Hall motion of the particle-like skyrmions in the perpendicular direction, and it is very much simplified compared to the adiabatic and non-adiabatic STT interactions in the LLG equation, which particularly depend on the microstructure of each skyrmion. The particle model leads to two weak universality classes, that is, the critical exponents for both perpendicular and parallel directions vary with the strength of the Magnus force. And the universality class in the parallel direction tends to the one of vortices when the Magnus force is small. In the present paper, we generally investigate the current-driven skyrmion dynamics based on the LLG equation with the STTs, and focus on the influence of the non-adiabatic coefficient b on the critical dynamics of skyrmions. The simulation results reveal the anisotropic critical behaviors in the perpendicular and parallel directions as shown in table 1, and the change of the critical exponents in the parallel direction is induced by the non-adiabatic STT. The values of the critical exponents are also very different from those of the particle model.
In the stationary state, the velocity of skyrmions exhibits a power-law behavior just above the critical current v∼τ β , where τ=(u − u c )/u c is the reduced current. With the theoretical solution in equations (18) and (19) for the velocities in two directions, the critical exponents β P and β ⊥ can be extracted, respectively. In figures 4(a) and (b), the steady velocities versus the reduce current τare shown for different values of the non-adiabatic coefficient b. The slopes of the curves give the static exponents β ⊥ and β P . β ⊥ is equal to 1.00 for both the adiabatic and non-adiabatic STTs, i.e. independent of b=0 and a non-zero b, and this result is consistent with the numerical simulation in table 1 based on the LLG equation in equation (2). For the parallel direction, β P =1.00 for the non-adiabatic STT, independent of the non-zero b. However, β P =2.00 for the adiabatic STT (b=0) is very different.
In contrast to the theoretical analysis, β P from the numerical simulation of the LLG equation decreases with the increasing non-adiabatic coefficient b, and there is a crossover from the universality class of the adiabatic STT to that of the non-adiabatic STT. This may be induced by the distortion of skyrmions. In the theoretical solution in equations (18) and (19), v x is proportional only to v d 2 for b=0 around the depinning phase transition, while v y relies on v | | d independent of b. This possibly results in two distinct universality classes for the parallel direction with and without the non-adiabatic STT, respectively. Moreover, the non-adiabatic STT and the distortion of skyrmions influence mainly the critical behavior of skyrmions in the driving direction, and only weakly in the perpendicular component of the Hall motion.
In table 2, the critical current and critical exponents for different strengths of disorder with a non-adiabatic coefficients b=0.5 are summarized. The critical current drops rapidly when the strength of disorder reduces. The static critical exponents β and ν decrease with the decreasing strength of disorder in both perpendicular and parallel directions, while the dynamic critical exponents z is opposite. Although the anisotropic critical behavior in two directions still exists, it becomes obviously weaker as the strength of disorder reduces. One may expect that such an anisotropy will eventually vanish in the absence of disorder.
The critical exponent β ⊥ for a weaker disorder deviates from the result of the theoretical analysis, in contrast to that for a stronger disorder. The critical exponent β P decreases with a weaker disorder, but it may become smaller than the theoretical result β P =1.00. The constant  in the theoretical solution can be estimated from u c of the adiabatic STT, which is equal to 0.126πQ, 0.0809πQ, 0.0411πQ for the strength of disorder K/J=1.0, 0.75, 0.5, respectively. With the disorder-dependent , one can also theoretically calculate the critical current u c for different strengths of disorder. Compared to figure 3(b), the deviation of the theoretical critical current from the numerical one becomes somewhat larger for a weaker disorder, e.g. about one percent. All these results imply that the theoretical solution presented in section 3.2 is only approximately valid for a strong disorder, where the pinning phenomenon could be described by a mean-field-like ansatz.
It is worth noting that, the skyrmion Hall angle Θ in the this film without defects is independent of the current density, while that with quenched defects increases with the increasing current density, and approaches a constant value calculated in the absence of the pinning effect. This drive-dependent skyrmion Hall effect around the depinning phase transition has been experimentally observed [21,23,28,41]. In our results, different critical behaviors in the directions perpendicular and parallel to the driving current indicate a nonlinear relation between the skyrmion Hall angle Θ and the driving current u, because of the ratio ( ) . In general, from the critical exponents given in table 1, the dependence of the skyrmion Hall angle on the driving current increases with the decreasing non-adiabatic coefficient b. When only with the adiabatic STT (b=0), | | ( ) 1.10 15 , and the dependence is the strongest. On the other hand, in consistence with the theoretical solution, | | in the LLG simulation may vanish when the non-adiabatic coefficient b is large enough, that is, the skyrmion Hall effect may be independent to the driving current around the depinning phase transition. Furthermore, the drive-dependent skyrmion Hall effect becomes weaker as the strength of disorder reduces, and it will vanish in the absence of disorder. The quantitative result of the anisotropic dynamic behavior around the threshold current explains well the drive-dependent skyrmion Hall effect in experiments. Our results indicate that a large non-adiabatic coefficient may not only efficiently reduce the threshold current of skyrmions, but also significantly weaken the drive-dependent Hall effect. And a weak disorder could also achieve a similar effect as the large non-adiabatic coefficient does. The lower drive dependence makes the manipulation of skyrmions more controllable in experiment.

Conclusion
With the LLG equation, we numerically simulate the nonstationary dynamic behaviors of skyrmions driven by currents in a magnetic thin film with quenched disorder. Based on the dynamic scaling form far from stationary, the critical current is convincingly located from the power-law behavior of the skyrmion velocity, and the static and dynamic exponents are then accurately determined with and without the non-adiabatic STT. Such a threshold current obtained in our simulation is insensitive to a small non-adiabatic coefficient b of the STT, but dramatically reduced for the large b. Importantly, the current-driven skyrmions exhibit anisotropic critical dynamic behaviors in the directions perpendicular and parallel to the driving current. Such an anisotropic characteristic provides a quantitative analysis of the drive-dependent skyrmion Hall effect in experiment around the depinning phase transition. The results for different strengths of disorder indicate that the anisotropy becomes weaker as the strength of disorder reduces. In practical applications, it may be a guidance to find the material with an appropriate non-adiabatic coefficient to manipulate more efficiently the skyrmions in the presence of disorder.
The theoretical analysis using the Thiele's approach with a phenomenological pinning force is also presented. The critical current u c and the static exponent β ⊥ in the perpendicular direction agree well with the simulation results, at least for a relatively strong disorder. The static exponent β P in the parallel direction is equal to 2.00 and 1.00 for the adiabatic and non-adiabatic STT, respectively. In the numerical simulation with the LLG equation, however, there exists a crossover inbetween, which may be caused by the distortion of skyrmions during the time evolution.