Skyrmion soliton motion on periodic substrates by atomistic and particle based simulations
J. C. B. Souza N. P. Vizarim C. J. O. Reichhardt C. Reichhardt P. A. Venegas


VECA: Reliable and Confidential Resource Clustering for Volunteer Edge-Cloud ComputingThis material is based upon work supported by the National Science Foundation (NSF) under Award Number: OAC-2232889. Any opinions, findings, and conclusions or recommendations expressed in this publication are those of the authors and do not necessarily reflect the views of the NSF.
    

Hemanth Sai Yeddulapalli1,
Mauro Lemus Alarcon1,
Upasana Roy1,
Roshan Lal Neupane1, 

Durbek Gafurov1,
Motahare Mounesan2,
Saptarshi Debroy2,
Prasad Calyam1

1University of Missouri-Columbia, USA; 
2City University of New York, USA

Email: 
1{hygw7, lemusm, u.roy, neupaner, durbek.gafurov, calyamp}@missouri.edu;

2mmounesan@gradcenter.cuny.edu; 
2saptarshi.debroy@hunter.cuny.edu




    
=======================================================================================================================================================================================================================================================================================================================================================================================================





§ INTRODUCTION

Solitons are nonlinear wave perturbations <cit.>
that have been observed in a wide range of different science fields,
including
mathematics <cit.>,
chemistry <cit.>,
magnetism <cit.>,
and biology <cit.>.
Due to their nonlinear spin dynamics,
magnetic systems are particularly well able to stabilize magnetic solitons
or nonlinear magnetic textures,
which can take the form of
magnetic vortices <cit.>,
magnon drops <cit.>,
magnetic skyrmions <cit.>,
hopfions <cit.>,
and bimerons <cit.>.

Solitons or kinks can also be stabilized in assemblies of
particles coupled to a
periodic substrate <cit.>.
When the number of particles equals the number of
potential minima, the system should be free of kinks;
however, if the number of particles
is slightly higher or lower so that the system is off commensuration,
localized kinks appear that depin under applied drive levels which can
be much lower than the drive at which the bulk of the particles depin.
Kink motion at incommensurate fillings on periodic substrates has been
studied for colloidal particles
<cit.>,
superconducting vortices <cit.> and various
friction models <cit.>.
Since skyrmions are also particle-like textures, when they are placed
on a periodic substrate, kinks or solitons could also be stabilized
in the skyrmion lattice.
Recent work by Vizarim et al.<cit.> has shown the
possibility of creating and moving a soliton along
quasi one-dimensional chains
of skyrmions using a particle based model.

After this study, Souza et al.<cit.>
demonstrated with an atomistic model
that soliton motion along skyrmion chains
is stable and that
the soliton exhibited higher velocities than
free skyrmions.
The work on the quasi-one-dimensional systems opened
the possibility of using soliton or kink motion in
magnetic skyrmions as an information transfer method
for new types of soliton-based devices
employing skyrmions. An open question
is whether soliton motion through skyrmion lattices
remains stable
in more realistic fully two-dimensional systems,
whether the soliton behavior can be captured using both particle-based
and atomistic models,
and where the two models agree or disagree.

Magnetic skyrmions are particle-like topologically
protected magnetic textures <cit.>
that exhibit many similarities to overdamped particles:
they minimize their repulsive interactions by forming a
triangular array, can be set in motion
by the application of external drives, and can interact
with material defects in a variety of ways
<cit.>.
The key difference between skyrmions and
other overdamped particles is the presence of a
non-dissipative Magnus force that causes
skyrmions to move in the absence of disorder at
an angle known as the intrinsic skyrmion Hall angle,
θ^int_sk, with respect to the external
driving force
<cit.>.
In order to simulate all of the degrees of freedom of a skyrmion,
it is necessary to
use computationally expensive models, such as 
the atomistic model <cit.>, that can capture
behavior such as skyrmion annihilation,
creation and deformation.
To mitigate the computational expense of skyrmion simulations,
Lin et al.<cit.> proposed
a particle-based model for skyrmions that
assumes the skyrmions remain rigid, an approximation that
is valid for low skyrmion densities and low external currents.

Using atomistic simulations and particle based simulations,
we study the dynamical behavior of soliton
motion in magnetic skyrmion lattices on square
and triangular substrates


just away from commensuration that are subjected to an
external driving force.


For the square substrate, both models produce
soliton motion along a 45^∘ angle;
however, the atomistic model exhibits an
additional 30^∘
soliton motion that is absent in the particle based model. At higher
drives, the entire skyrmion lattice depins, and the transitions
between the the different soliton and skyrmion flow phases are visible
as changes in the transport curves and average Hall angle of
the kink or skyrmion motion.
For the triangular substrate, we also find regimes of stable soliton
motion, but the models
show substantial differences.

The atomistic model exhibits soliton motion
along a 30^∘ angle for a wide range of external
driving forces, whereas the particle model produces
soliton motion along a 45^∘ angle for a small
range of external driving forces.
In the particle model, the trajectory of the soliton is
much more meandering, resulting in flow around an average
angle of 
45^∘, while
in the atomistic model,
the finite size of the skyrmions reduces the amount of
meandering flow that occurs, causing motion along
30^∘ to be more stable.



§ METHODS





We model Néel skyrmions in thin films
with a magnetic field applied perpendicular
to the film at zero temperature, T=0 K, with
periodic boundary conditions along the x and y directions.
We use two distinct models, the atomistic model
and the particle based model. Common
substrate defect arrangements of N_m defects are used in both
models. The square array of defects is modeled as
ϕ_S(x, y)=A/4[cos(2π x/a_0)+cos(2π y/a_0) + 2],
where A is the defect strength and a_0 is the substrate
lattice constant. The values of A and a_0 are different
between each model and are listed below in the subsections describing
the individual models. 
The triangular array of defects is modeled as
ϕ_T(x, y)=∑_i=1^3A/2[cos(2π b_i/a_0) + 1], with b_i=xcos(θ_i)-ysin(θ_i)+a_0/2 and
θ_1=π/6, θ_2=π/2, θ_3=5π/6.
As in the square array, A is
the defect strength and a_0 the substrate lattice constant, with
different
values of A and a_0 used for each model.
For both models we choose values of a_0 such that
there are N_m=36 minima in the defect arrangement.
We select the number of skyrmions N_sk to be just above
commensuration with the substrate,
N_sk=N_m+1=37.
In fig. <ref>(a,b), we show a three-dimensional rendering
of the square and triangular defect arrays. Each potential minimum captures a
single skyrmion, and we add one additional skyrmion to the sample in
order to create
a kink or soliton
that can depin at a much lower drive
than the commensurate skyrmions. 



 §.§ Atomistic Simulations

The atomistic model tracks the state of individual
atomic magnetic moments <cit.>.
The Hamiltonian describing the interactions of
thin films at T=0K with an underlying square
atomic arrangement with lattice parameter a=0.5 nm is
given by<cit.>

    ℋ=   -∑_i, j∈ NJ_ij𝐦_i·𝐦_j
                    -∑_i, j∈ N𝐃_ij·(𝐦_i×𝐦_j) 
       -∑_iμ𝐇·𝐦_i
                    -∑_i K(x_i, y_i)(𝐦_i·𝐳̂)^2  .


The first term on the right side is the exchange interaction
between the nearest neighbors contained in the set N,
with an exchange constant of J_ij=J between magnetic moments
i and j.
The second term is the interfacial Dzyaloshinskii–Moriya
interaction, where 𝐃_ij=Dẑ×𝐫̂_ij is the Dzyaloshinskii–Moriya
vector between magnetic moments i and j and 𝐫̂_ij
is the
unit distance vector between sites i and j.
The third term is the Zeeman interaction with an applied external magnetic
field 𝐇.
Here μ=ħγ is the magnitude of the magnetic moment
and γ=1.76×10^11T^-1s^-1 is the electron
gyromagnetic ratio. The last term represents the 
perpendicular magnetic anisotropy (PMA) of the sample, where
x_i and y_i are the spatial coordinates of the i
magnetic moment. We use
K(x_i, y_i)=ϕ_S(x_i, y_i) for a square array of defects
and K(x_i, y_i)=ϕ_T(x_i, y_i) for a triangular array of defects.
In ultrathin films,
long-range dipolar interactions can be neglected <cit.>.

The time evolution of the individual
atomic magnetic moments is given by the LLG
equation <cit.>


    ∂𝐦_i/∂ t=-γ𝐦_i×𝐇^eff_i
                                 +α𝐦_i×∂𝐦_i/∂ t
                                 +pa^3/2e(𝐣·∇)𝐦_i  .

Here
𝐇^eff_i=-1/ħγ∂ℋ/∂𝐦_i
is the effective magnetic field including all interactions from
the Hamiltonian, α is the Gilbert damping, and the last
term is the spin-transfer-torque (STT), where 
p is the spin polarization,
e is the electron charge,
and 𝐣=j𝐱̂ is the applied current density.
The STT current assumes that the conduction electron
spins are always parallel to
the magnetic moments 𝐦 <cit.>, and
the driving force from the STT current <cit.>
acts perpendicular to 𝐣.
We fix
α=0.3, J=1 meV, D=0.5J, and μ𝐇=0.6(D^2/J)(-𝐳̂).
The resulting skyrmions 
move at an intrinsic skyrmion Hall angle of
θ_sk^int=64^∘ with respect to the driving
force exerted by external currents.
For both
the square and triangular defect arrays
we use A=0.1J and a_0=14 nm.
The sample dimensions are 84 nm × 84 nm for the square
array of defects, and (2/√(3))84 nm × 84 nm for the
triangular array of defects. The difference in
sample size is required to properly apply boundary conditions.



 §.§ Particle Based Simulations

The particle based simulations are governed by the
equation of motion<cit.>

    α_d𝐯_i+α_m𝐳̂×𝐯_i=∑_i≠ j𝐅_sk(𝐫_ij)+𝐅_P(𝐫_i)+𝐳̂×𝐅_D  ,

where 𝐯_i is the velocity of skyrmion i.
The first term on the left side is the damping term,
α_d, which can be written<cit.>
as α_d=-α𝒟,
where α is the Gilbert damping and 𝒟 is the
dissipative tensor. The second term is the Magnus force
where the Magnus strength α_m can be
written<cit.>
as α_m=-4π Q, where Q is the skyrmion charge. The
ratio α_m/α_d determines the intrinsic skyrmion
Hall angle θ_sk^int=arctan(α_m/α_d).
In order to match
the particle based θ_sk to the atomistic
θ_sk,
we use values of α_m and α_d such
that θ_sk^int=arctan(α_m/α_d)=64^∘.
The first term on the right side of eq. <ref> is the skyrmion-skyrmion
interaction given by
𝐅_sk(𝐫_ij)=-U_skK_1(r_ij)𝐫̂_ij, where U_sk=1 is the
interaction strength and K_1 is the modified first order Bessel function.
The second term is the interaction with the underlying
substrate potential, given by
𝐅_P(𝐫_i)=-∇ϕ_S(𝐫_i)
for the square array of defects and
𝐅_P(𝐫_i)=-∇ϕ_T(𝐫_i)
for the triangular array of defects.
In both defect arrays, the potential strength is A=4 and
the substrate lattice constant is a_0=6.
The last term is the interaction with an external drive,
𝐅_D=F_D𝐱̂, which is in accordance
with the action of an STT current on
magnetic skyrmions <cit.>.
Our simulation box is of size
36 × 36 for the square defect array
and (2/√(3))36×36 for the triangular defect array.



§ SQUARE ARRAY





We first compare the dynamics of skyrmions
interacting with a square array of defects in atomistic and
particle-based simulations.
In fig. <ref>(a), we plot ⟨ v_x⟩ and ⟨ v_y⟩
versus applied current j, and in fig. <ref>(b), we show the
corresponding effective θ_sk=arctan(⟨ v_y⟩/⟨ v_x⟩)
versus j from the atomistic simulation.
We observe four dynamical phases: a pinned phase I with no
skyrmion motion,
phase II_a where a soliton moves along 45^∘,
phase II_b in which a soliton
moves at 30^∘, and phase V in which
all of the skyrmions depin and move without locking to any direction.
The skyrmion motion in phases
II_a and II_b is illustrated in fig. <ref>(a, b).




Figure <ref>(c, d) shows ⟨ v_x⟩, ⟨ v_y⟩,
and θ_sk versus F_D from particle based simulations
of skyrmions
interacting with a square array of defects.
We again observe four dynamical phases, which are a
pinned phase I,
phase II_a with soliton motion along 45^∘,
phase IV_a in which all of the skyrmions
depin and move along
θ_sk=45^∘,
and phase V where all
of the skyrmions have depinned but move without locking to
any direction.
The skyrmion motion in phases II_a and IV_a is illustrated
in fig. <ref>(c, d).

Figure <ref>(a) shows the skyrmion trajectories
during the phase II_a soliton motion along
45^∘ from
the atomistic model.
Here, the extra skyrmion shares an anisotropy minimum with another
skyrmion, and the skyrmion-skyrmion interaction force between the
two lowers the depinning force at which one of the skyrmions can
escape from the minimum and become an interstitial skyrmion.
The interstitial skyrmion moves across the anisotropy landscape
until it reaches another anisotropy
minimum containing a pinned skyrmion.
Through skyrmion-skyrmion interactions,
the interstitial skyrmion pushes the pinned skyrmion
out of the anisotropy minimum and takes up residence in the minimum.
The displaced skyrmion becomes the new
interstitial skyrmion, and this pattern of motion repeats indefinitely.
In fig. <ref>(b) we illustrate the
skyrmion trajectories during the phase II_b soliton
motion along 30^∘ using the atomistic model.
The mechanism
of motion is identical to that observed in fig. <ref>(a),
but now the interstitial skyrmion encounters a pinned
skyrmion along the 30^∘ angle line instead of
the 45^∘ angle line.

Figure <ref>(c) shows the skyrmion trajectories during the
phase II_a soliton motion at 45^∘ obtained from
the particle based model. The
behavior is identical to that found in the atomistic model,
where the skyrmion-skyrmion interactions cause the extra skyrmion
trapped inside a potential minimum to depin at low
F_D.
Since the skyrmion is treated as
a point particle, the trajectory differs in detail from the
atomistic trajectory shown
in fig. <ref>(a); however,
the mechanism of exchange between interstitial and pinned skyrmions
remains the same.
In fig. <ref>(d) we show
the skyrmion trajectories in phase IV_a from the
particle based model
where all of the skyrmions have depinned and
move along 45^∘.
Here the extra skyrmion does not play
a major role, since all
of the skyrmions are 
interacting with the potential in an ordered manner and following
a 45^∘ trajectory. The skyrmion lattice is similar
to a moving crystal but contains a localized defect produced by
the extra skyrmion.

We note that previous work on skyrmions
moving over a two-dimensional square periodic  
substrate under an increasing drive
showed a directional locking effect in which
the skyrmion motion locked to particular symmetry angles of the 
substrate <cit.>.
Continuum models for individual
skyrmions on antidot lattices
also produce similar directional 
locking
<cit.>. 
The results we describe here are different
in that the motion is not of individual
continuously moving skyrmions
but of kink traveling through a skyrmion lattice,
so the locking is a collective effect instead of a single
particle effect.
Additionally, the applied drive needed
to produce kink motion is substantially lower
that the drive at which an isolated skyrmion would depin from
the substrate potential minimum. The reduction in the
driving
threshold would be particularly useful for applications. 



§ TRIANGULAR ARRAY





We next compare the dynamical behavior of skyrmions on
a triangular defect array
in atomistic simulations and particle based simulations.
For the atomistic simulations,
fig. <ref>(a) shows ⟨ v_x⟩ and ⟨ v_y⟩
versus j, while in fig. <ref>(b) we plot the corresponding
θ_sk versus j.
We observe three dynamical phases.
At low drives, we find a pinned phase I.
There is a transitional phase III_a of disordered soliton
motion, where soliton transport occurs but does not follow a
well defined direction.
At higher drives, phase II_b appears in which the soliton
moves along
30^∘. Over the range of j values that we consider, we
never observe a completely depinned state, but we expect that for
larger values of j,
the pinned skyrmions would eventually escape from
the anisotropy minima and move throughout the sample.
The skyrmion trajectories
in phases III_a and II_b are illustrated
in fig. <ref>(a, b).




In fig. <ref>(c), we plot ⟨ V_x⟩ and
⟨ V_y⟩ versus F_D for a particle based simulation
of skyrmions driven over a triangular defect array, and
in fig. <ref>(d) we show the corresponding
θ_sk versus F_D.
Five dynamical phases appear.
At low drives, we find a pinned phase I.
Phase II_a is soliton motion along
45^∘, and it is followed by a transitional phase III_b in which
a small number of skyrmions are present simultaneously and move
chaotically across the system.
In phase IV_b,
all of the skyrmions are moving along 30^∘,
and phase V consists of
all of the skyrmions moving without locking to any direction.
Illustrations of the skyrmion trajectories
for phases II_a and IV_b appear in fig. <ref>(c, d).

Unlike what we found for the square defect array,
here the atomistic model and the particle based model
do not exhibit good agreement.
A pinned phase is present in both cases, but as we increase
the drive, in the atomistic simulations
we observe a wide phase II_b 30^∘ soliton
motion over the
range
4.8×10^10A m^-2≤ j≤ 12×10^10A m^-2.
In comparison, for the particle based simulations
there is a small window of phase II_a
45^∘ soliton motion over the range
0.66≤ F_D≤ 0.88.
The phase II_b 30^∘ angle soliton motion is
absent in the particle based simulations.
This large difference in behavior is likely due to the fact that
in the atomistic simulations, the skyrmion has a
finite size and is able to deform, as is visible by
comparing the phase II_a flow in fig. <ref>(a)
to that in fig. <ref>(c), or
comparing the phase II_b flow in fig. <ref>(b)
to the phase II_a flow in fig. <ref>(c).
For higher values of F_D, the behavior of the atomistic and
particle based simulations diverge significantly.
The particle based simulation
produces phase IV_b flow in which
all of the skyrmions move along 30^∘, followed by
phase V flow in which all of the skyrmions move without locking
to a particular direction,
while the atomistic simulations remain trapped in phase II_b
with soliton motion along 30^∘.

Figure <ref>(a) shows the trajectories of skyrmions
interacting
with a triangular defect array from atomistic
simulations performed at j=4×10^10A m^-2, corresponding to
phase III_a in fig. <ref>(a, b). The soliton mechanism of motion
described previously still occurs here, but the soliton does not
lock to any angle and gradually works its way all around the sample.
When we increase the external current to j=8×10^10A m^-2,
corresponding to phase II_b in fig. <ref>(a, b),
the skyrmions move as illustrated in fig. <ref>(b).
Again, we observe a soliton motion, but unlike what was shown
in fig. <ref>(a),
the motion occurs along a well defined angle of 30^∘.
This motion
appears
over the range 4.8×10^10A m^-2≤ j≤ 12×10^10A m^-2.

In fig. <ref>(c) we illustrate the skyrmion trajectories for motion
over a triangular defect array from
particle based simulations
with F_D=0.8, corresponding to the phase II_a flow
in fig. <ref>(c, d).
For this value of
F_D, a soliton moves across the sample at 
45^∘.
We observe the 45^∘ soliton motion only
in the particle based simulations and do not find it in the
atomic simulations.
When we increase F_D to F_D=1.5,
corresponding to phase IV_b in fig. <ref>(c, d),
all of the skyrmions flow along 30^∘,
as shown in fig. <ref>(d).
Similar to what we found in
fig. <ref>(d), the skyrmion lattice travels as a moving
crystal that contains
a localized defect produced by
the extra skyrmion.

The difference in the angle of soliton motion on a triangular
defect array between the atomistic and particle based models
can be explained
by the finite skyrmion size in
the atomistic simulations.
The barrier between anisotropy minima is enhanced when the skyrmion
has a finite size, whereas in the particle based simulations,
the pointlike nature of the skyrmions gives a reduced barrier
between anisotropy minima.
The larger barrier potential created by the finite
skyrmion size can be observed
by comparing the
trajectory of an interstitial skyrmion as it
passes between two anisotropy minima;
the trajectories in the atomistic simulation
exhibit fewer meanders compared to the trajectories in the particle
based simulations.
The increase in the barrier potential
forces the atomistic interstitial
skyrmion to move following the 30^∘ angle imposed by the
triangular lattice,
whereas the interstitial skyrmion in the
particle based model can travel along a wider range of paths
between substrate minima.
To compensate for the overly unconstrained mobility of
skyrmions the particle
based model, the barrier between potential minima can be increased
either by
increasing
the skyrmion-skyrmion interaction strength U_sk or
reducing the lattice constant of the substrate potential.
In principle, it may be possible to identify simulation parameters
for the particle model of the triangular substrate that would match
the behavior of the atomistic simulations and compensate for the
rigidity and vanishing size of the particle-based skyrmion model.



§ SUMMARY


We compared the results of atomistic simulations and particle based
simulations of soliton motion for skyrmion assemblies just past
commensuration on both square and triangular substrates.
For the square array, both models agree well at low and high
drives, but differ for intermediate drives.
The atomistic model produces
a pinned phase, soliton motion along 45^∘,
soliton motion along 30^∘, and unlocked flow of all skyrmions.
The soliton motion proceeds via the replacement
by an interstitial skyrmion of
a pinned skyrmion in an anisotropy minimum, with the depinned skyrmion
becoming the new interstitial skyrmion.
The particle based model produces
a pinned phase, soliton motion along 45^∘,
a phase where all of the skyrmions move
along 45^∘, and unlocked flow of all skyrmions.
The particle based model
does not exhibit the 30^∘ soliton motion
found in the atomistic
simulations, and the trajectories in the particle based model meander
more than the atomistic model trajectories due to the rigidity and
vanishing size of the particle based skyrmions.
For both models, the different dynamic phases are visible
as signatures
in the velocity-force curves, skyrmion Hall angle, and skyrmion trajectories. 

The atomistic and particle based models do not agree well on the
motion of skyrmions over
a triangular defect array.
The atomistic model produces
a pinned phase, a transitional phase in which a soliton moves with
no well defined angle,
and a regime of soliton motion at 
30^∘ that
spans a wide range of external drive values.
The particle based model exhibits 
a pinned phase, soliton motion along 45^∘,
a transitional phase in which all skyrmions participate in disordered
soliton motion,
a phase in which all of the skyrmions move
along 30^∘, and unlocked motion of all of the skyrmions.
Here only the pinned phase is common
between the two models.

Our results provide a better understanding of the
regimes in which the particle model is a good or a poor approximation
for the skyrmion motion.
We argue that it can be possible to mitigate the approximations
made in the particle based model by adjusting the strength of the
interactions between the skyrmions or modifying the lattice constant
of the substrate.
Our results will be beneficial for determining how to
control skyrmion soliton 
motion using a combination of anisotropy trapping and
external driving.

We gratefully acknowledge the support of the U.S. Department of
Energy through the LANL/LDRD program for this work.
This work was supported by the US Department of Energy through the Los Alamos National Laboratory. Los
Alamos National Laboratory is operated by Triad National Security, LLC, for the National Nuclear Security
Administration of the U. S. Department of Energy (Contract No. 892333218NCA000001). 

J.C.B.S acknowledges funding from Fundação de Amparo à Pesquisa do Estado de São Paulo - FAPESP (Grant 2023/17545-1).

We would like to thank Dr. Felipe F. Fanchini for providing the computational resources used in this work. 
These resources were funded by the Fundação de Amparo à Pesquisa do Estado de São Paulo - FAPESP (Grant: 2021/04655-8).

eplbib


