Bilayer skyrmion dynamics on a magnetic anisotropy gradient

Magnetic skyrmion transport has been primarily based on the use of spin torques which require high current densities and face performance deterioration associated with Joule heating. In this work, we derive an analytical model for energy efficient skyrmion propagation in an antiferromagnetically-coupled bilayer structure using a magnetic anisotropy gradient. The interlayer skyrmion coupling provides a strong restoring force between the skyrmions, which not only prevents annihilation but also increases their forward velocity up to the order of km s–1. For materials with low Gilbert damping parameter, the interlayer skyrmion coupling force can be amplified up to ten times, with a corresponding increase in velocity. Furthermore, the analytical model also provides insights into the dynamics of skyrmion pinning and relaxation of asymmetric skyrmion pairs in bilayer-coupled skyrmion systems.


Introduction
Magnetic skyrmions are particle-like chiral magnetization states that are promising as an information carrier due to their attractive characteristics such as nanometre dimensions, and topological stability [1][2][3][4][5][6][7][8][9][10][11]. In the implementation of a magnetic skyrmion racetrack memory, the transport techniques explored thus far have been mainly limited to the use of spin-polarized electrical current injection which utilizes spin transfer torque (STT) and spin-orbit torque (SOT). Experimental transport measurements based on such mechanisms showed that high current densities on the order of 10 12 A m −2 are required to achieve practical speeds of about 100 m s −1 [12,13]. The elevated temperatures at such high current densities can drastically reduce the propagation speed of skyrmions while the risk of skyrmion annihilation increases significantly [14][15][16].
A wide range of alternative propagation techniques has been proposed not only to overcome the issue of Joule heating but also for applications in insulating skyrmions. These alternatives include magnetic field gradients [17][18][19], electric field gradients [20], temperature gradients [21][22][23], spin waves and magnons [24][25][26]. These systems face various challenges in their implementation and scalability, while spin waves face attenuation and scattering about each skyrmion along a long nanowire.
The use of a magnetic anisotropy gradient for skyrmion propagation is favoured for its ease of electrical integration and control [27][28][29][30][31]. A magnetic anisotropy gradient is known to exert a force to propagate skyrmions towards the region with lower Ku [32][33][34]. By applying an electric field across the ferromagnetic layer, its magnetic anisotropy can be modulated due to the voltage controlled magnetic anisotropy (VCMA) effect [35][36][37].
In our work, we derive an analytical model for the dynamics of VCMA gradient-driven skyrmion pair in synthetic antiferromagnetic bilayer structure. The synthetic antiferromagnetic bilayer structure aims to surpass the velocity limit of the monolayer structure due to the skyrmion Hall effect [8,[38][39][40]. A restoring force that is almost ten times greater than the original driving arising from the interlayer exchange coupling raises the maximum possible velocity of skyrmions to the order of km/s. The VCMA gradient transport operates by generating a magnetic anisotropy strength K u ( )gradient across the skyrmion using an external electric field and shifting the K u gradient along the wire as illustrated in figure 1. The magnetic anisotropy potential generates the VCMA gradient force (F VGF ) which drives the skyrmion from the region of high K u to the region of low K .
u The sawtooth configuration of the K u gradient confines the skyrmion to the region of lower K u as it propagates with the K u gradient. The K u gradient and its propagation along the wire is achieved by the application of the corresponding electric field across the region. As this system only utilizes voltage across the track to generate a gradient for propagation as opposed to a continuous flow of current for STT, less energy is required for its operation.
As shown in figure 1(b), the bilayer structure consists of two antiferromagnetically-coupled magnetic layers separated by a spacer layer. The driving force in this system F VGF acts along the wire and towards the region of lower magnetic anisotropy strength. The trajectory of a skyrmion is generally not along the direction of the driving force due to the skyrmion Hall effect generated by the Magnus force or gyroscopic force. The opposite sign of the skyrmion topological number in the top and bottom magnetic layer results in the gyroscopic force acting in opposite directions, hence the skyrmion in the top layer is driven towards the upper edge while the skyrmion in the lower layer is driven towards the lower edge. The bilayer system introduces a bilayer coupling force (F BCF ) between the two skyrmions on each layer as they are driven apart. The skyrmions achieve steady motion along the wire when the gyroscopic force that acts perpendicular to velocity is mitigated by F BCF and the dissipative force that acts against velocity is balanced by F . VGF This results in the trajectory as seen in figure 1(b).

Methods
Numerical solutions for the skyrmion propagation were evaluated by micromagnetic simulation using MuMax3. The time evolution of the magnetization state was resolved following the Landau-Lifshitz-Gilbert  field H eff is given by the summation of effective fields that includes external magnetic fields, demagnetization field, exchange field, thermal field, Dzyaloshinskii-Moriya interaction (DMI) and anisotropy field. Material parameters used for micromagnetic simulation are exchange stiffness A exIntra =15 pJ m −1 , a=0.1 , M s =580 kA m −1 , and DMI strength D=3 mJ m −1 . These material parameters are used as the control for investigations on the dependence skyrmion velocity on individual material parameters. A range of values for perpendicular magnetocrystalline anisotropy, K u =0.80-1.00 MJ m −3 values were used where skyrmions can be nucleated and sustained [3,45]. The interlayer exchange stiffness was set at the value of A ex Inter =0.45 pJ m −1 which is possible by strong antiferromagnetic coupling spacer such as Ru [46][47][48].
The simulation cell size was set at 1 nm by 1 nm by 0.4 nm. The bilayer was constructed by having two magnetic layers separated by a single non-magnetic layer of 0.4 nm thickness which represents a material such as Ru or Ir to antiferromagnetically couple the adjacent ferromagnetic layers. The gradient magnetic anisotropy is formed by varying the values of magnetic anisotropy cell by cell along the slope.

Analytical model for VCMA gradient driven skyrmions
The dynamics of bilayer skyrmions on a VCMA gradient observed above can be modelled by the Thiele equation using the corresponding forces acting on the skyrmions. Assuming that the magnetization profile is rigid, the skyrmion equation of motion can be modelled by the Thiele equation derived from Landau-Lifshitz Gilbert (LLG) equation [4,27,42,[49][50][51]. The Thiele equation is given by, where G is the gyromagnetic coupling vector, v is the drift velocity of the centre of mass of the skyrmion, D  is the dissipative tensor, d is the thickness of the magnetic layer, and F is the force acting on the skyrmion.
The first term in equation (2) is referred to as the gyroscopic force which acts perpendicular to the direction of propagation v. This force results in the skyrmion hall effect (SHE) which is the deviation in direction of motion from its driving force. The first term is also known as the Magnus force term due to the resemblance in the force generated when a spinning body travels through a viscous fluid that acts perpendicular to the velocity of the body [50,51].
) where Q is the topological number of the skyrmion. The value of Q can be either +1 or −1 depending on the chirality and core polarity of the skyrmion which can be either in the z +ˆor z -ˆdirection.
The second term in equation (2) is the dissipative force which opposes the motion of the skyrmion and can be inferred from the opposing direction of this force with v. The dissipative tensor is given by which depends on the skyrmion magnetization profile. The third term is a general term used to describe any external forces acting on the magnetic skyrmion where the force E F = - is given by the energy gradient. In the case of our work on bilayer skyrmions on a VCMA gradient, the two forces acting on the skyrmions are VCMA gradient force (F VGF ), and bilayer coupling force (F BCF ) as shown in figure 1(b). The edge repulsion force is not considered in our calculations as the skyrmion motion investigated are all distant from the edges of the wire.
Equation (2) can be decomposed into two equations given by In these two equations, the skyrmion Hall effect manifests as a non-vanishing v y when F x is non-zero and F y is zero. The skyrmion Hall angle q which is the angle at which the skyrmion propagates relative to the direction of driving force F, can be given by In a bilayer skyrmion system, F F The bilayer coupling force provides a force in the transverse direction that mitigates the skyrmion hall effect leading to v 0 m s . where the ratio of the forces is From equations (5) and (6), skyrmion motion parallel to the wire is dependent on both longitudinal and transverse forces. As a is typically very small, the ratio of F y to F x is large, and the skyrmion Hall angle is large. In our bilayer structure where 82.8 , | and the bilayer coupling force becomes the dominant contributor to v .
x Equation (6) shows the inverse relationship between skyrmion velocity and .
a By being strongly coupled to material parameters, a huge opportunity is presented for material optimization for high skyrmion speeds.
3.1. VCMA gradient force VGF arises from the inhomogeneous distribution of magnetic anisotropy energy generated using VCMA. The magnitude of VGF can be quantified by where K u is the magnetic anisotropy strength, m z is the magnitude of magnetization along the z-axis, and d is the thickness of the magnetic layer. From equations (8) and (9), F VGF acts towards the region of lower K u as observed in the simulation results due to the term which is dependent on skyrmion size that decreases with increasing K u [27]. Hence, the VCMA gradient favours and exerts a stronger force on larger skyrmions.

Bilayer coupling force
Bilayer coupling force F BCF arises from the interlayer exchange coupling of the bilayer skyrmions which forms an energy potential well which is a function of the distance between the centres of the overlapping skyrmions. The potential well provides the restoring force between the skyrmions when they are driven apart. Total exchange energy of the system is given by the integral over the volume of the magnetic layer where A ex is the exchange stiffness [52]. Assuming that the skyrmion profile is rigid, the x m 2 2 ¶ ¶ / and y m 2 2 ¶ ¶ / terms do not have any dependence on the lateral separation between the skyrmions. Therefore, the interlayer exchange energy is calculated using only the interaction between the cells directly above or below. For the antiferromagnetic coupling between the layers, A ex Inter is negative. The second order derivative is then simplified to where m 1 is the normalised magnetization on layer 1, m 2 is the normalised magnetization on layer 2, h is the separation distance between layer 1 and layer 2.
where r is the radial distance from the centre of the skyrmion, and R is the radius of the skyrmion. Our magnetization profile approximation shows good agreement within the range of K u investigated while at lower K , u the magnetization profile deviates from our approximation to resemble that of a chiral bubble domain instead.
Using the fitted magnetization profile given in equation (12) where s is the effective separation between the two coupled skyrmions. As the skyrmions are driven apart and separation increases, the skyrmions experience a minor distortion in their profile near its maximum magnetization away from the other skyrmion. The effective separation between the two skyrmions can be calculated by taking the weighted position of the skyrmion across its centre. Figure 2 shows the presence of a restoring force. Within the range of separation below half their radius, the force profile resembles a harmonic oscillator where the gradient of the force is negative and linear. For this range of separation, the restoring force is independent of skyrmion radius and given by figure 2, as the bilayer skyrmions are increasingly driven apart by a larger gyroscopic force, the restoring force reaches a maximum at s kR = where k 1.15 = is a dimensionless constant. This is intuitive based on the origin of the force which relies on the overlap of the skyrmions that are finite in size. Beyond separation of kR, F BCF gradually decreases to zero as the skyrmions become increasingly decoupled and eventually cease to interact with each other when they no longer overlap. Using equation (14) at the separation of kR, we find the relationship between the theoretical maximum restoring force and skyrmion radius to be

As shown in
Using the material parameters considered for this work, the maximum velocity that can be sustained is 10.3 km s −1 . This maximum velocity can be increased when F BCF is complemented by other forces which act towards the centre of the wire such as edge repulsion force. Edge repulsion force that acts only in the vicinity of the edge requires a wire with a sufficiently thin width to be effective in complementing F . BCF Enhancement of the edge repulsion force to prevent annihilation can be achieved by engineering the edges using additional magnetic material and VCMA [32,33,[53][54][55].

Skyrmion velocity dependence on magnetic anisotropy strength and gradient
We investigated the dependence of the skyrmion speed on both the steepness of the K u gradient and the value of K u at the centre of the skyrmion. The range of K u gradient investigated ranges from 500 GJ m −4 to 10 TJ m −4 which corresponds to K u change of 500 J m −3 and 0.01 MJ m −3 per nm respectively. This range is limited by the optimal range of K u required to stabilise the skyrmions. The size of skyrmions which is dependent on K u results in minor expansion of the skyrmions when they are subjected to an increasing K u gradient.  (14) is represented by a solid line. figure 3 verifies the linear relationship between skyrmion speed and K u gradient shown in equation (8). Figure 3 also shows a decrease in speed with increasing K u value at the centre of the skyrmion. This is due to the decrease in F VGF which decreases with skyrmion size as discussed previously.

Results presented in
The fitted solid and broken line represents the analytical solution derived using F VGF and F BCF from equation (8) and equation (14) respectively. The accurate agreement between the numerical and analytical solutions supports the validity of our model for both F VGF and F . BCF The accuracy of our model for F BCF is limited by the magnetization profile approximation used in equation (12) that is accurate for high K u skyrmions.
The maximum K u gradient investigated at 10 TJ m −4 is possible based on the large K u variation of up to 0.29 MJ m −3 observed experimentally [56]. While this observation has yet to be observed in the system of our interest, it nonetheless supports the possibility for the realization of such a K u gradient.

Skyrmion velocity dependence on the material parameter
VCMA-driven skyrmion velocity dependence on Gilbert damping parameter and saturation magnetization was investigated at K u =0.80 MJ m −3 and K u gradient=2000 GJ m −4 . Figure 4 shows the strong inverse relationship between skyrmion velocity and a which is consistent with equation (6). The decrease in a from 0.1 to 0.01 resulted in an increase in velocity from 120 m s −1 to 1200 m s −1 shown in figure 4, surpassing the velocities shown in figure 3. Even at such velocity, the skyrmion pair remains coupled and protected from edge annihilation. equation (17) predicts that F BCF can overcome the Magnus force up to velocity of 10.3 km s -1 . This relationship makes a a critical material parameter to optimize for higher skyrmion velocity. The range of material parameter investigated was restricted to skyrmions which did not undergo severe deformation during propagation and retained an eccentricity of less than 0.2.
However, the dependence of skyrmion velocity on saturation magnetization is not inversely related as equation (6) suggests. This disagreement can be explained by the increase in the skyrmion size with saturation magnetization that results in an increase in driving force by the K u gradient. The increase in skyrmion size can be attributed to the increase in demagnetization with M , s which in turn decreases the effective out-of-plane  magnetic anisotropy. As shown in the inset of figure 4, skyrmion diameter increases more than linearly with saturation magnetization and this results in the velocity of equation (6) to have a relationship with a saturation magnetization of the order of more than 1. Hence, skyrmion velocity was observed to increase with saturation magnetization.
The variation of skyrmion velocity that increases with DMI strength, as shown in figure 5, is the result of the increase in skyrmion diameter similar to the dependence observed with saturation magnetization. While saturation magnetization only results in a small increase in velocity, DMI holds a stronger dependence and results in more than double increase in velocity when DMI increases from 2.6 to 3.2 mJ m −2 . However, the magnetic skyrmion is only stable within the DMI range of 2.4 to 3.6 mJ m −2 , thus voltage-modulation of DMI is an important consideration when utilizing VCMA as it can potentially destabilize and annihilate the skyrmions [57,58].

Skyrmion velocity field in a bilayer structure
The F BCF quantified in equation (14) is the restoring force exerted by each skyrmion on the other to restore the completely overlapped state. Applying F BCF to the Thiele's equation given in equation (2), the velocity field of a skyrmion due to the F BCF exerted by the other skyrmion can be acquired. The trajectory of a skyrmion in the proximity of another pinned skyrmion follows a spiral configuration that ends when the skyrmions lie directly over one another. Figure 6 shows the velocity field of a skyrmion due to the bilayer coupling force from a pinned skyrmion in the other layer. The result is shown in figure 6 can be qualitatively derived from the direction of the force which acts radially towards the centre of the skyrmion and the skyrmion hall angle given in equation (5) which describes the relative angle between the driving force and velocity. The colour in figure 6 represents the strength of the bilayer coupling force with separation between the skyrmions as shown previously in figure 2.
In the bilayer structure, the trajectory due to bilayer coupling force between a freely moving skyrmion (skyrmion 1) on one layer and another skyrmion (skyrmion 2) on the other layer, which has different material parameter are shown in figure 7. From equation (6), skyrmion velocity is material parameter dependent and an asymmetry in those parameters such as saturation magnetization and Gilbert damping parameter can result in spiralling trajectories. The case where the difference in the material parameters results in the ratio of the velocity of skyrmion 1 to skyrmion 2 to be 0.2, 0.8, and 1.0 are presented in figure 7. When both skyrmions have identical material parameters (1.0), the two skyrmions converge towards each other in a straight trajectory. Increasingly asymmetric or massive skyrmion 2 for the cases of 0.8 and 0.2 multiple results in the trajectory to become increasingly spiral and centred around skyrmion 2. At the limit where skyrmion 2 is stationary, the trajectory is as shown in figure 6, where skyrmion 2 does not move while skyrmion 1 spirals towards the centre of skyrmion 2 to the point where they overlap.

Outlook
The experimental realisation of this work will first require the observation of magnetic skyrmions in a SAF configuration which has yet to be reported. Beyond the search for skyrmion-stable SAF systems, several other considerations must be considered. The formation of a continuous electric field gradient can be achieved by adopting the multiplexed 3-level racetrack design with high, average and low voltage step-wise gradient presented by Wang et al and Liu et al [27,31]. In this design, the dimensions of each repeating unit of K u gradient along the racetrack is dependent on the skyrmion size, maximum K u variation, and the size of the electrodes. The electrode width is optimized at approximately the diameter of the skyrmion to achieve best propagation efficiency [31]. The maximum K u gradient applicable is then dependent on the maximum K u variation of the system. In addition, DMI and RKKY interactions can both be tuned by external voltages. From the recent work by Suwardy et al which showed a large DMI variation of 1.3 mJ m −2 , the voltage-induced DMI can play a significant role in affecting the stability of the skyrmions [58]. Hence, the maximum voltage and the corresponding anisotropy gradient which can be applied may be limited.
The bilayer skyrmion pair can also be driven using a K u gradient in only one of the ferromagnetic layers instead of both ferromagnetic layers. This has been verified by our simulation and details are available in appendix .

Conclusion
The numerical results from our work on magnetic anisotropy gradient driven skyrmions in an antiferromagnetically coupled bilayer structure have shown a significant superiority over the monolayer structure due to the restoring force between the skyrmions which arises from interlayer antiferromagnetic exchange interaction. This restoring force confines the skyrmions to travel along the proximity of the wire centre and prevent edge annihilation. For the range of separation between the coupled skyrmions of less than half its radius, the restoring force resembles a harmonic oscillator which motivates the possibility of nano-oscillator and other applications. The driving force generated by the magnetic anisotropy gradient is dependent on skyrmion size which can be maximised with a low magnetic anisotropy strength or high saturation magnetization. Skyrmion velocity can also be maximised by utilizing its inverse dependence on Gilbert damping parameter. From our model of the restoring force, we derive the velocity field of a skyrmion which returns a spiralling configuration and explored the trajectory of two skyrmions with discrepancies in material parameters which returned straight or spiralling trajectories. Dynamics of asymmetric synthetic antiferromagnetically coupled skyrmions with material parameter discrepancy can be investigated further in the future.