Suppression of the skyrmion Hall effect in synthetic ferrimagnets with gradient magnetization

Magnetic skyrmions are promising building blocks for future spintronic devices. However, the skyrmion Hall effect (SkHE) remains an obstacle for practical applications based on the in-line transport of skyrmions. Here, we numerically study the static properties and current-driven dynamics of synthetic ferrimagnetic skyrmions. Inspired by graded-index magnonics, we introduce a linear gradient of saturation magnetization (Ms) in the skyrmion-hosting sample, which effectively modulates the skyrmion Hall angle and suppresses the SkHE. Micromagnetic simulations reveal that ferrimagnetic skyrmions could exhibit greater susceptibility to the variation of Ms as compared to their ferromagnetic counterparts. The Thiele analysis is also applied to support the simulation results, which elucidates that the Ms gradient dynamically modifies the intrinsic normalized size of skyrmions, consequently impacting the SkHE. Our results pave the way to the graded-index skyrmionics, which offers novel insights for designing ferrimagnet-based skyrmionic devices.

Following the principles inspired by graded-index optics [63], researchers have proposed a continuous modulation of magnetic parameters in spintronics [64][65][66].Although previous studies have delved into scenarios involving magnetic anisotropy gradients [67][68][69] and Dzyaloshinskii-Moriya interaction (DMI) gradients [70][71][72]to control skyrmion dynamics, models with saturation magnetization (M s ) gradients have only been extensively investigated in the context of magnonics [73][74][75][76][77], but have not yet been explored in the study of skyrmion dynamics.The M s of FiM materials exhibit higher sensitivity to temperature changes compared to FM materials, which provides a natural advantage in constructing M s gradients through thermal landscapes [78,79].
In light of the current fervor in both FiM skyrmions and graded-index spintronics, in this work, we numerically investigate the skyrmion dynamics in synthetic FiM with a M s gradient.It is demonstrated that introduction of the M s gradient can effectively regulate the skyrmion Hall angle.A theoretical analysis further reveals that the regulation is attributed to the variation of normalized skyrmion radius.The paper is organized as follows.In Sec.II, we describe the micromagnetic and analytical model.In Sec.III, we show and discuss the results of the numerical studies, which are summarized in Sec.IV.Specifically, in Sec.III A, we verify how static skyrmion configuration varies with changes in material parameters, preparing for subsequent dynamic studies.In Sec.III B, two types of current injection geometries are compared, and an overall trajectory of skyrmion motion induced by spinorbit torque (SOT) is provided.The details of skyrmion dynamics are presented, demonstrating an effective suppression of SkHE in the presence of the M s gradient.In Sec.III C, we analyze the physical mechanisms behind the simulation results based on the Thiele approach.

II. MODEL AND METHODOLOGY
The system considered in our simulations is a synthetic FiM Gd/Co/Pt multilayer [58][59][60][61].As depicted in Figs. 1 (a) and (b), the magnetization in Gd and Co layers are oppositely aligned due to the AFM coupling, leading to the opposite polarities of skyrmions in the two layers.This kind of magnetization configuration is known as the exchangecoupled bilayer-skyrmion [37,41,58].As shown in Fig. 1 (c), a nanoplate geometry with dimensions of 300 × 200 nm is considered to investigate skyrmion dynamics.This nanoplate includes regions at left and right sides with constant M s values (i.e.M s and M * s ), each spanning a width of 50 nm.In the central square region of 200 × 200 nm, spatially varying M s values with a linear gradient ∆ M s are introduced along the x-axis.Periodic boundary conditions are imposed in the y-direction to remove the potential boundary effect.With a single isolated skyrmion initially placed at the left side, a spin-polarized current is applied to drive this skyrmion into motion towards the +x direction.
The total Hamiltonian H of the system includes intralayer terms H L intra for respective layers (L = Co, Gd) and an interlayer term H inter , which is given by The intralayer Hamiltonian for each layer reads where m L i = 1 represents the normalized local magnetic moment at site i, and ⟨i, j⟩ extends to all the nearest neighbor sites in each layer.The first term represents the FM Heisenberg exchange interaction with A being the intralayer exchange constant.The second term represents the perpendicular magnetic anisotropy (PMA) with K being the anisotropy constant.The third term represents the DMI, where D i j is the DMI constant and ν i j is the unit vector between sites i and j.The fourth term H d represents the dipole-dipole interaction.On the other hand, the interlayer Hamiltonian which describes the interaction between Co and Gd layers reads [31] where c z is the thickness of the cell size, and σ is the bilinear surface exchange coefficient.
To explore the current-induced dynamics of FiM skyrmion, we numerically solve the Landau-Lifshitz-Gilbert (LLG) equation augmented with a spin-torque term τ, where γ is the gyromagnetic ratio, α is the Gilbert damping coefficient, and is the effective local field with µ 0 being the the vacuum permeability constant.We consider two strategies for the injection of spin-polarized currents [80].For the current-in-plane (CIP) geometry, the current flows through all the layers where the spin-transfer torques (STTs) are formulated in the Zhang-Li form [81]: Here, the first term is the adiabatic torque, and the second term is the nonadiabatic torque with β being the degree of nonadiabaticity.The STT coefficient is given by u = (gµ B PJ)/(2eM s ) [82], where g, µ B , P, e and J are the Landè factor, Bohr magneton, spin polarization factor, elementary charge and applied current density, respectively.On the contrary, for the currentperpendicular-to-plane (CPP) geometry, the charge current flows through the Pt layer and leads to an out-of-plane spin current due to the spin Hall effect [80].In this case, dampinglike spin-orbit torque (SOT) works, which is formulated as [83]: where θ H is the spin Hall angle, and m p is the normalized spin polarization vector.Here, the field-like SOT is excluded as its impact on the skyrmion dynamics is negligible [37].
The simulations are performed by MuMax3 finitedifference GPU accelerated code [84].The simulation volumes for static and dynamic cases are 100 × 100 × 1 nm 3 and 300 × 200 × 1 nm 3 , respectively, with a cell size of 0.5 nm, matching the lattice constant of Co to maximize the calculation accuracy.The default simulation parameters and their ranges of variation are given in Table I.

A. Static FiM skyrmion without a M s gradient
We begin by verifying the configuration of a static FiM skyrmion in a system without a M s gradient.For this purpose, we utilize a square geometry system of 100 × 100 nm 2 with infinite boundaries.A bilayer skyrmion is initialized in the system and the equilibrium state is then determined using the conjugate gradient method.Firstly, while maintaining the default M s , we investigate the impact of variation in material parameters A, D, and K on the stability and size of the skyrmion.The results are depicted in Figs. 2 (a)-(c), where in each figure two parameters are varied with the remaining parameter being fixed at its default value.Overall, the skyrmion radius R varies from a few nanometers to approximately 30 nm within the selected range, which is positively correlated with D, while negatively correlated with A and K.When skyrmion size becomes excessively large or small, the system tends to change into multi-domain or uniform states.It is noteworthy that the stability of skyrmion does not necessarily require such a large PMA because similar configurations can be obtained by applying a vertical external magnetic field [88].Based on the above results, we set A = 15, D = 2.5 and K = 2.0 as default values, resulting in an initial skyrmion radius of 5.5 nm.
Secondly, we fix A, D and K at their default values and vary M s within a small range of ±5% (i.e., M Co s ranging from 1.386 to 1.512 while M Gd s maintains half of the M Co s value) [60] to observe variation in the skyrmion configuration.The m z profiles along the diameter direction under different M s are depicted in Fig. 2 (d).The dots are extracted from simulation data, while the solid curves are results of fitting with [86] where w is the width of the 360°domain wall of the skyrmion, considered as a fitting parameter, and R/w is the normalized skyrmion radius.Given that our simulation results agree well with the theoretical formula, we thus extract R and R/w from Fig. 2 (d) and plot their variations with respect to M s in Figs. 2 (e) and (f).For comparison, calculated values of R in a singlelayer FM are also plotted, where a notably higher sensitivity of FiM skyrmion size to M s is observed.Note that R/w also increases with increasing M s , and approximately exhibits a linear relationship, which has an important meaning in the subsequent theoretical analysis (see Sec. III C).

B. Suppression of the skyrmion Hall effect
We first consider the case with constant M s .As described in Sec.II, both CPP and CIP injection geometries are initially considered for current-driven skyrmion motion.For CIP injection, we further examine two cases with different nonadiabatic degrees, β = 0 and β = α.The dependence of skyrmion velocity v =∥v∥ on the current density J for different scenarios is shown in Fig. 3 (a).The overlapping symbols indicate that skyrmions in the two layers always move as bounded entities.This means that the interlayer surface exchange σ is sufficiently strong, thereby preventing their decoupling under high current.Moreover, the CPP injection turns out to be of much greater efficiency in driving skyrmions compared to the CIP injection, which is consistent with reported findings on FM skyrmions [80] and AFM skyrmions [32].Therefore, in the following sections, we only consider the SOT-induced skyrmion motion for the sake of simplicity, with a moderate density of 1 GA/m 2 .
Next, we introduce the M s gradient in the central region of the nanoplate.The linear ∆ M s is defined as where ∆ M s varies between ±0.05 and is controlled by changing M * s at the right-side region.For various values of ∆ M s , snapshots of the driven skyrmions near the left and right markers (yellow chain lines) of the gradient region are presented in Fig. 3 (b).At 8.5 ns, the skyrmion reaches the left marker and begins to move along different trajectories under the influence of different ∆ M s .We find significant variations in the skyrmion size and velocity during the moving process.Positive ∆ M s leads to an expansion of the skyrmion and a decrease in velocity, while negative ∆ M s results in its contraction and an increase in velocity.The trajectories of the skyrmion present a scattered distribution.When ∆ M s = 0, there is a small but nonzero skyrmion Hall angle Φ, consistent with a previous report [58].As ∆ M s varies, the trajectory exhibits a monotonic and regular change, providing visual evidence that the magnetization gradient can regulates the skyrmion Hall angle Φ.
To learn more details of the skyrmion motion under the M s gradient, we further extract the data of skyrmion dynamics, which display the variations in skyrmion position Y (X), radius R, and velocity components v x and v y as functions of the position X in Fig. 4. Here, to provide more precise outcomes, we only display the cases where ∆ M s varies from −0.02 to +0.03, within the central gradient range.The skyrmion center is determined by averaging the position coordinates of the 360°domain wall with m z = 0.The results in Fig. 4 (a) are similar to those in Fig. 3 (b), but we can observe that the trajectories of the skyrmion center are remarkably smooth.It is noteworthy that there is almost no displacement in the ydirection when ∆ M s = +0.01. Figure 4 (b) reveals linear dependence between R and X with a slope determined by ∆ M s , where positive (negative) ∆ M s tends to induce the expansion (shrinkage) of the skyrmion.Importantly, the dynamically changing values of R agree well with those in the equilibrium states in Fig. 2 (e).This agreement is expected when the skyrmion velocity is relatively low because the driven skyrmion can be relaxed to approach the minimal-energy configuration at equilibrium at each moment.Therefore, we believe that the findings about the static configuration in Fig. 2 (e) are also applicable for the analysis of the driven skyrmions.The velocity components, v x and v y , are crucial parameters in determining Φ.Hence, we plot them separately in Fig. 4 (c).When ∆ M s is positive (negative), both v x and v y increase (decrease) as compared with those at ∆ M s = 0.The differ- ence appears in their spatial-position dependence, that is, v y is nearly independent of the position, whereas v x exhibits nearly linear dependence on X.
After obtaining the skyrmion instantaneous velocity v(v x , v y ), we can then calculate the skyrmion Hall angle Φ at each moment by As the motion of skyrmion is a continuous process, we are more interested in its total displacement over a certain period of time.Therefore, we further calculate the average value of Φ throughout the entire motion, presenting it alongside the instantaneous values in Fig. 5 (a). Figure 5 (a) is a violin plot with the instantaneous and average values of Φ being represented by hollow and solid symbols, respectively.The curves on both sides of the violin shape display the probability den-

C. Analysis based on Thiele equation
In this section, we theoretically analyze of the simulation results.While the size of a Néel skyrmion may change, its central symmetry in shape remains unchanged, which allows us to model it as a rigid point-like particle.Therefore, its motion can be qualitatively elucidated by a modified Thiele equation [89] The three terms arise from the precessional, damping, and SOT terms in the LLG equation, respectively.Here, G = nẑ is the gyromagnetic coupling vector, D is the dissipative tensor, and the tensor B quantifies the efficiency of the spin Hall torque over the two-dimensional configuration.The parameters n and d correspond respectively to the topological invariant as a sum of the solid angles of magnetizations, and the component of the disspation tensor as an indicator of the magnetization rotation length-scale, given by For a FM skyrmion with a simple profile, d is approximately given by d ≈ π 3 R/w [23], while for a FiM skyrmion with R roughly as large as w, it is given by d ≈ 2πR/w [57].This means that d ∝ R/w is expected to scale likewise with the normalized skyrmion radius R/w presented in Fig. 2 (e).In our model, due to the presence of spatially dependent ∆ M s , the skyrmion internal configuration R/w changes with X, resulting in the variation of d = kR/w during the skyrmion motion, where k is the proportionality coefficient.The solutions of Eq. ( 10) are where B 0 is a constant derived from the tensor B. Equations ( 13) and ( 14) can be used to explain the results of instantaneous velocity in Figs. 4 (c) and (d).When α is small enough, the term α 2 d 2 can be neglected, so that v y becomes independent of d, manifested in Fig. 4 (d) as v y almost unrelated to X.As for v x , the term αd cannot be neglected, so that v x is related to d, i.e., correlated with R/w, as seen in Fig. 4 (c), where v x exhibits an approximate linear correlation with X.
Regarding the variation of the skyrmion Hall angle Φ, it follows the relation The above Eq.( 15) can be used to explain the simulation results shown in Fig. 5 (a).The original variable ∆ M s affects the internal skyrmion configuration R/w, causing an incremental change ∆ d during the motion process.This further leads to an increment of ∆ Φ.When ∆ M s takes an appropriate value, around +0.01 in this work, ∆ Φ precisely compensates for the original Φ, thus suppressing the skyrmion Hall effect.To visually demonstrate the relationship in Eq. ( 15), we plot the function of ∆ tan Φ with respect to ∆ w/R, as shown in Fig. 5 (b).This function reveals the correlation between the changes in Φ and the normalized skyrmion radius.We find that the simulation results for different ∆ M s are uniformly distributed around the theoretical line, and the linear relationship conforms to the description provided by Eq. ( 15).

IV. CONCLUSION
In conclusion, we have investigated the static properties and current-induced dynamics of the synthetic FiM skyrmions.Compared to the FM skyrmions, equilibrium FiM skyrmions exhibit a higher sensitivity of their static profiles to the variation in M s .Inspired by the graded-M s in the field of gradedindex magnonics, we have explore the influence of a M s gradient on the skyrmion dynamics by particularly focus on the SOT-induced motion.It is shown that the presence of a M s gradient significantly affects the trajectory, velocity, and size of the driven skyrmion.Notedly, a 5% variation in M s gradient can control the skyrmion Hall angle within a range of ±10°.Specifically, a 1% M s gradient completely suppresses the SkHE.We have analyzed the simulation results using Thiele theory, demonstrating that the M s gradient regulates the skyrmion Hall angle by influencing the normalized skyrmion radius.In this study, we have incorporated the M s gradient from the magnonics to the skyrmionics.Our findings are expected to provide a new perspective to overcome the SkHE, which has been an obstacle to the technical applications of skyrmions so far.

FIG. 1 .
FIG. 1.(a) Schematic configuration of a ferrimagnetically exchangecoupled bilayer skyrmion.Arrows represent the magnetization, and the color represents their z-component m z .(b) Cross-sectional view of the Gd/Co/Pt heterojunction.Only the magnetization in free layers are explicitly modeled.(c) Schematic diagram of the geometry for skyrmion dynamics.The central region has a linear magnetization gradient ∆ M s , and the left and right regions have constant M s .

FIG. 2 .
FIG. 2. (a)-(c) Variations of equilibrium skyrmion radius R versus material parameters A, D, and K.The color represents the skyrmion size, and the white region indicates the absence of a stable skyrmion.The red circles mark the situation with default parameters, i.e., K = 2.0, D = 2.5, and A = 15.(d) Magnetization component m z plotted as a function of position X along the skyrmion diameter under different M s .The symbols are simulation results, and the solid lines are fitted curves.(e) Variations of skyrmion radius R for both FM and FiM cases versus M s .(f) Variations of the ratio of skyrmion radius to domain wall width R/w versus M s .The symbols are simulation or analytical results, and the dashed lines are used to guide the eyes.

FIG. 3 .
FIG. 3. (a) Skyrmion velocity v as functions of the driving current density J for both CPP and CIP injections.The Co and Gd layers are represented by hollow and solid symbols, respectively.(b) Snapshots of the SOT-induced skyrmion dynamics for various ∆ M s .The dashed lines separate regions with constant M s or finite M s gradient, and the arrows indicate the direction of skyrmion motion.The color wheel indicates the in-plane magnetization.∆ M s and time information are indicated in white text.

FIG. 4 .
FIG. 4. (a) Variations of the skyrmion center coordinates (X,Y ) in the central region during the motion process.(b) Variations of the skyrmion radius R versus position X for various ∆ M s .(c) Variations of the skyrmion velocity components, v x and v y , versus position X.The dashed lines represent v x , and the solid lines represent v y .The color of the line indicates different ∆ M s .

FIG. 5 .
FIG. 5. (a) Violin plot with hollow and solid symbols depicting the instantaneous and average skyrmion Hall angles Φ for various ∆ M s .The side-profile lines represent the probability density distribution of instantaneous values.The dashed line illustrates the variation of Φ versus ∆ M s .(b) The increment of skyrmion Hall angle ∆ tan Φ plotted as a function of the increment of the reciprocal of the normalized skyrmion radius ∆ w/R.The symbols represent simulation results, with colors indicating various ∆ M s .The solid line is a linear fitted line based on Eq. (15).

TABLE I .
Simulation parameters.
[60]r simplify, M Gd s is always set to an estimated value of 50% of M Co s[60].b Values enclosed in square brackets denote the ranges of variation.c For simplify, an average value of g = 2.1 is employed for both layers.