Directional locking and the influence of obstacle density on skyrmion dynamics in triangular and honeycomb arrays

Overdamped particles driven over a periodic substrate can exhibit a series of directional locking steps of preferred motion. Depending on the substrate symmetry, the particles can flow in a direction which differs from that of the applied external drive. For a square obstacle array, locking motion appears along the angles ϕ = arctan(p/q), where p and q are integers. The most notable locking steps appear at ϕ = 0°, 45° and 90°. The locking steps differ for other obstacle array symmetries. For example, in triangular arrays, the most prominent locking angles are ϕ = 30° and 60°, and locking motion appears at angles of $\phi =\mathrm{arctan}\left(\sqrt{3}p/\left(2q+1\right)\right)$. Locking phenomena has been widely studied using colloidal particles [1–6], type II superconducting vortices [4, 7, 8], classical electrons [9], and more recently using magnetic skyrmions [10, 11] and active matter [12]. In the case of colloidal particles, the colloids are placed within a periodic array of posts or optical traps, and an external drive is applied. Locking is observed as either the driving force or the substrate is rotated. The width of the locking steps can be modeled as a function of the details of the colloid–substrate interaction such as the shape and size of the obstacles or traps. Directional locking also appears for quasiperiodic substrates, which exhibit five or seven locking directions [13, 14].

Most directional locking studies have focused on overdamped systems, and produce locking steps by changing the angle of an applied external drive or by rotating the substrate. For magnetic skyrmions, due to the Magnus force it is possible to observe directional locking with the driving direction and substrate orientation both fixed. Skyrmions in chiral magnets can be found in numerous materials and their sizes can vary from a micron to 10 nm [15]. Recently, stable skyrmions were found in materials at room temperature [16–19]. When skyrmions are subjected to an external drive they can be set into motion and undergo a transition from a pinned to a sliding state [15, 20–23]. They can be measured in experiments by direct observation of skyrmion motion [16, 24, 25], changes in the topological Hall effect [20, 26], x-ray diffraction [27] or neutron scattering [28].

Skyrmions are promising for many technological applications. For example, skyrmions could be used as bits for information storage in logical and memory devices [29, 30] where the bit state is associated with the presence or absence of a magnetic skyrmion. The main difference between skyrmions and overdamped particles is the important role played by the Magnus force, which may dominate the dynamics of a skyrmion system [31]. It is believed that the Magnus force is also the origin of the low value of the observed threshold depinning current for skyrmions [32, 33]. In the absence of defects in the sample, the skyrmion flows at an angle with respect to the external drive known as the intrinsic skyrmion Hall angle ${\theta }_{\text{sk}}^{\text{int}}$ [15, 22, 34]. This angle changes as a function of the ratio between the Magnus and the damping term. Experimentally observed skyrmion Hall angles range from a few degrees up to values close to 90°, but the higher angles may be achieved only in certain systems [34–38]. The intrinsic skyrmion Hall angle can be a problem for certain technological applications, such as race track devices, since it limits the maximum distance a skyrmion may flow before contacting the edge of the device. Thus, controlling the skyrmion motion is of significant interest in the scientific community [15]. One way to achieve such control is by introducing artificial pinning centers or obstacles to the sample. The skyrmion Hall angle in the presence of pinning is close to 0° at low drives, but it can increase in quantized steps as a function of the applied drive [10, 39]. For sufficiently high drives, the skyrmion Hall angle saturates to a value very close to the intrinsic Hall angle, ${\theta }_{\text{sk}}^{\text{int}}$ [10, 34, 39–42]. Most studies of skyrmions interacting with pinning involve randomly placed defects that act like obstacles or pinning centers; however, nanostructuring techniques make it possible to place the defects at specific positions, enabling the creation of ordered pinning arrangements [43, 44].

In this work we simulate for the first time the behavior of a single skyrmion interacting with triangular and honeycomb obstacle arrays, at zero temperature, while being driven with an external force that is applied in the x direction. We investigate the influence of the obstacle density on the system, which is a parameter that has been explored in few works [45, 46]. We also study the behavior for varied values of the obstacle density and the ratio of the Magnus to the damping term. We find that the skyrmion behavior is strongly dependent on the obstacle density, and that some dynamical regimes appear or grow in magnitude while others vanish. For both arrays, the skyrmion exhibits strong directional locking at θ_sk = −30° and −60° due to the substrate symmetry. We plot dynamical phase diagrams that summarize this strong locking effect and its stability. Our results indicate that it is possible to use the obstacle density as a way to control the skyrmion motion.

We consider the dynamics of a single skyrmion in a two-dimensional system of size L × L with periodic boundary conditions in the x and y directions. The skyrmion interacts with either a triangular or a honeycomb obstacle array, illustrated in figure 1. We obtain the skyrmion dynamics using a particle-based model for skyrmions [23], shown in equation (1), using an in-house simulation code based on standard molecular dynamics techniques [47].

$$\begin{equation}{\alpha }_{\mathrm{d}}{\mathbf{v}}_{i}+{\alpha }_{\text{m}}\hat{z}{\times}{\mathbf{v}}_{i}={\mathbf{F}}_{i}^{o}+{\mathbf{F}}^{D}.\end{equation} \tag{ 1 }$$

Zoom In Zoom Out Reset image size

In this equation, the first term on the left is the damping arising from the spin precession and dissipation of electrons localized in the skyrmion core, where α_d is the damping constant. The second term on the left represents the Magnus force, where α_m is the Magnus constant. The Magnus term produces a force that is perpendicular to the skyrmion velocity. The first term on the right of equation (1) is the interaction between the skyrmion and the obstacles. We model this potential with the Gaussian form [10, 11, 48] ${U}_{o}={C}_{o}{e}^{-{\left({r}_{io}/{a}_{o}\right)}^{2}}$, where C_o is the strength of the obstacle potential, r_io is the distance between skyrmion i and obstacle o, and a_o is the obstacle radius. Thus, the force between the obstacle and the skyrmion takes the form ${\mathbf{F}}_{i}^{o}=-\nabla {U}_{o}=-{F}_{o}{r}_{io}{e}^{-{\left({r}_{io}/{a}_{o}\right)}^{2}}{\hat{\mathbf{r}}}_{io}\enspace$, where ${F}_{o}=2{C}_{o}/{a}_{o}^{2}$. Throughout this work we set a₀ = 0.65 and F₀ = 1. For computational efficiency we set a cut-off at r_io = 2.0; beyond this length the interaction is negligible. All distances are normalized by the screening length ξ, and obstacle densities are given in terms of 1/ξ². The last term in equation (1) is the external dc drive, ${\mathbf{F}}^{D}={F}^{D}\hat{\mathbf{d}}$, where $\hat{\mathbf{d}}$ is the direction of the applied external driving force. In this work we fix $\hat{\mathbf{d}}=\enspace \hat{\mathbf{x}}$. We measure the skyrmion velocity parallel, $\left\langle {V}_{{\Vert}}\right\rangle$, and perpendicular, $\left\langle {V}_{\perp }\right\rangle$ to the drive. For a skyrmion moving in the absence of obstacles, in the overdamped limit of α_m/α_m = 0, the skyrmion flows in the same direction as the driving force. For finite values of α_m/α_m, however, the skyrmion motion follows the skyrmion Hall angle, ${\theta }_{\text{sk}}=\mathrm{arctan}\left(\left\langle {V}_{\perp }\right\rangle /\left\langle {V}_{{\Vert}}\right\rangle \right)\enspace =\mathrm{arctan}\left({\alpha }_{\text{m}}/{\alpha }_{\mathrm{d}}\right)\enspace$. The external driving force is increased in small steps of δF = 0.001 and we wait 10⁵ simulation time steps between drive increments to ensure that the system has reached a steady state. We normalize the damping and Magnus coefficients according to ${\alpha }_{\mathrm{d}}^{2}+{\alpha }_{\text{m}}^{2}=1$.

In figures 2(a) and (b) we show ⟨V_||⟩, ⟨V_⊥⟩, and θ_sk as a function of F_D for a sample with a triangular obstacle array where ρ_t = 0.128 and α_m/α_d = 1.0. In the absence of obstacles, the skyrmion flow would be along the intrinsic Hall angle, ${\theta }_{\text{sk}}^{\text{int}}=-45{\circ}$. As shown in figure 2(b), when obstacles are present the skyrmion direction of motion exhibits a series of quantized locking steps up to F^D = 1.7, followed at higher drives by a continuous saturation of θ_sk to the intrinsic Hall angle. For F^D ⩽ 0.307 the skyrmion motion is locked to θ_sk = 0°, as shown in figure 3(a), and figure 2(a) indicates that over this range of drives, ⟨V_||⟩ is increasing with F^D while ⟨V_⊥⟩ = 0. For 0.307 < F^D ⩽ 1.041 the skyrmion locks to θ_sk = −30°, which is a preferred direction of motion due to the symmetry of the triangular obstacle array. Here the skyrmion moves by one lattice constant in the x direction for every lattice constant in the y direction, giving ${\theta }_{\text{sk}}=\mathrm{arctan}\left(\sqrt{3}p/\left(2q+1\right)\right)=\mathrm{arctan}\left(\sqrt{3}/3\right)=30{\circ}$. The corresponding skyrmion trajectory in figure 3(b) flows around the upper contour of one obstacle in each row before proceeding to the next row. In figure 3(c) we illustrate the motion for locking to θ_sk = −36.56° at F^D = 1.1. In this case the skyrmion trajectories are straighter due to the higher velocities. At F^D = 1.5 in figure 3(d), the skyrmion trajectories are along θ_sk = −40.8°. For higher values of the applied drive, the skyrmion direction of motion continuously approaches the intrinsic skyrmion Hall angle.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

In addition to producing steps, the locking is also visible in the time series of the particle velocities. In figure 4(a) we show the parallel and perpendicular velocities at F^D = 0.5 when the system is in a directionally locked phase with θ_sk = −30°, while figure 4(b) shows the same quantities at F^D = 2.5 in a non-step region. In the locked phase, the velocities are exactly periodic with the same signal appearing during each cycle. In the non-locked region, there are more oscillations for the same time because the drive is higher; however, the patterns are no longer exactly repeating, although there is still a periodicity in the signal due to the periodic array. The differences can be seen more clearly by examining the Fourier transform S(ω) of the velocity signals for the x direction velocity component V_||, as shown in figure 5. The peaks in S(ω) are much stronger in the locking regime than in the non-locking regime. This suggests that if an ac driving is added to the dc driving, the locking step regions should show strong Shapiro step phenomena whereas in the non-locking regions, the Shapiro step effect will be absent or reduced.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Turning next to the honeycomb obstacle lattice, in figure 6 we plot ⟨V_||⟩, ⟨V_⊥⟩, and θ_sk versus F^D for a system with ρ_h = 0.123 and α_m/α_d = 1.0. There is a pinned phase for F^D ⩽ 0.076 when the skyrmion is trapped between the obstacles. For higher F^D , the skyrmion depins and moves along θ_sk = −30°, as illustrated in figure 7(a). Note that unlike the triangular array, there is no regime with θ_sk = 0° since the vacancies in the honeycomb lattice trap the skyrmion at low drives. For F^D > 1.294, there is a series of transitions in the skyrmion motion, which ultimately saturates at a value close to the intrinsic Hall angle when F^D = 3.0. As an example of the motion found on these closely spaced short steps, in figure 7(b) we show the skyrmion trajectory at F^D = 1.7 where θ_sk = −40.8°.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

To study the influence of the obstacle density on the skyrmion dynamics, we fix α_m/α_d = 1.0 and consider five samples with different obstacle densities for both the triangular and honeycomb arrays. For each sample we then separately vary α_m/α_d in order to determine the influence of the Magnus force on the behavior. In figure 8, we plot θ_sk versus F^D for varied obstacle densities for the triangular and honeycomb arrays, showing a clear influence of obstacle density on the skyrmion dynamics. In figure 8(a), at very low values of obstacle density such as ρ_t = 0.032, the skyrmion motion is very close to the intrinsic skyrmion Hall angle, which for α_m/α_d = 1.0 is ${\theta }_{\text{sk}}^{\text{int}}=-45{\circ}$. That is, the skyrmion flow is nearly unaffected by the presence of the widely spaced obstacles. As more obstacles are added to the sample, the skyrmion dynamics become richer and multiple locking steps of different sizes emerge. For the triangular array there are three main locking directions of 0°, −30°, and −40.8°. The θ_sk = 0° locking occurs only for higher obstacle densities of ρ_t > 0.128, while the θ_sk = −30° locking is present for all density values simulated in this work, although its extent is limited to 0 ⩽ F^D ⩽ 0.05 when ρ_t = 0.032. The θ_sk = −40.8° locking step is also present for all density values, but it is more robust for low densities. In the honeycomb array, we observe only two main locking directions, θ_sk = −30° and −40.8°. The θ_sk = −30° locking is very robust for all density values, but the phase with θ_sk = −40.8° is not very robust and becomes weak at high obstacle densities.

Zoom In Zoom Out Reset image size

We use the results obtained in figure 8 to construct dynamic phase diagrams of the main locking phases as a function of obstacle density ρ_t or ρ_h versus F^D , as shown in figure 9. In both cases, the extent of the pinned phase increases as the obstacle density is increased. The pinned phase is wider in the honeycomb lattice due to the trapping of the skyrmion inside the lattice vacancy sites. For the triangular lattice in figure 9(a), the θ_sk = 0° phase is absent for very low obstacle densities and increases in width up to a maximum at ρ_t = 0.288. In contrast, the θ_sk = −40.8° phase is largest at low obstacle densities. For the honeycomb lattice in figure 9(b), the θ_sk = −30° step retains a nearly uniform width over the entire range of obstacle densities considered, although the location of the step shifts to higher F_D with increasing ρ_h. The θ_sk = −40.8° phase is considerably diminished in extent compared to the triangular array.

Zoom In Zoom Out Reset image size

In figure 10 we plot θ_sk versus F^D for different values of α_m/α_d at varied obstacle densities for both the triangular and honeycomb arrays. Note that the α_m/α_d = 1.0 case is already shown in figure 8.

Zoom In Zoom Out Reset image size

Figure 10 indicates that the locking angles which appear change as α_m/α_d is modified. When α_m/α_d = 0.5, the intrinsic skyrmion Hall angle is ${\theta }_{\text{sk}}^{\text{int}}=-26.57{\circ}$. As a result, the skyrmion can lock strongly to the substrate symmetry direction of θ_sk = −30° even though $\vert {\theta }_{\text{sk}}\vert { >}\vert {\theta }_{\text{sk}}^{\text{int}}\vert$ for this particular case. This obstacle-induced motion at an angle with a magnitude larger than that of the intrinsic skyrmion Hall angle was observed in previous works for a square obstacle array [10, 39], but here we show that for a triangular obstacle array it depends upon the obstacle density. As illustrated in figure 10(a), the θ_sk = −30° locking step only exists for very low obstacle densities, such as ρ_t < 0.072, and disappears at higher obstacle densities. In figure 11(a) we show the skyrmion trajectory for the θ_sk = −30° step at F^D = 1.5, ρ_t = 0.072, and α_m/α_d = 0.5. In figure 10(a), the number of dynamic phases which appear is related to the obstacle density. At ρ_t = 0.128, the system has a very rich set of dynamic phases, while for other densities the number of dynamic phases is reduced. For very low ρ_t the skyrmion motion is very close to the intrinsic Hall angle value, so there are few dynamic phases. On the other hand, when ρ_t becomes too high, the reduced spacing between adjacent obstacles pinches off many possible skyrmion trajectories, eliminating the corresponding dynamic phases. We find that for a given value of F^D , varying the obstacle density can cause the skyrmion motion to follow a variety of different directions. For example, at F^D = 0.75, the skyrmion can move along θ_sk = 0° or θ_sk = −30° depending on the obstacle density. For ρ_t ⩽ 0.072, the skyrmion locks to θ_sk = −30°, while for ρ_t > 0.072, it locks to θ_sk = 0°. This opens the possibility of designing devices in which regions with distinct obstacle densities coexist in order to force the skyrmion to follow a designated trajectory. The pinning geometry is important, however, as figure 10(b) indicates that in a honeycomb obstacle array, there is only a single dynamic state of motion along θ_sk = −30° regardless of the obstacle density.

Zoom In Zoom Out Reset image size

At α_m/α_d = 2.0, shown in figures 10(c) and (d), ${\theta }_{\text{sk}}^{\text{int}}=-63.44{\circ}$. For the triangular array in figure 10(c), when ρ_t ⩽ 0.072 the skyrmion depins directly into the θ_sk = −60° step, illustrated in figure 11(b). This is a very stable angle of motion for the triangular array due to its symmetry. When ρ_t ⩾ 0.128, additional steps emerge at θ_sk = 0° and −30°. For the honeycomb lattice in figure 10(d), when ρ_h ⩽ 0.266 the skyrmion depins into the θ_sk = −30° state and follows the symmetry direction of the honeycomb lattice. In contrast, at ρ_h = 0.342 the skyrmion depins into the θ_sk = −40.8° state shown in figure 11(c). The number of dynamic phases increases as the obstacle density increases. For example, the θ_sk = −49° step illustrated in figure 11(d) is absent for ρ_h ⩽ 0.123. In figures 10(e) and (f) we show samples with α_m/α_d = 3.0, where the intrinsic Hall angle is ${\theta }_{\text{sk}}^{\text{int}}=-71.57{\circ}$. For each array, the locking steps at θ_sk = −30° and −60° are significantly reduced in size. As the magnitude of the Magnus term increases, the skyrmion moves at a larger angle with respect to the driving direction, reducing the robustness of locking steps along the preferred symmetry directions of the obstacle arrays.

Both the triangular and honeycomb obstacle arrays should have preferred directions of motion along θ_sk = −30° and −60° due to the array symmetry. In this section we quantify the robustness of these dynamic phases as a function of α_m/α_d for arbitrary obstacle densities. We also explore the possibility of topological selection, since different species of skyrmions with distinct Magnus forces may coexist in a given sample [17, 19, 49–53]. Our results are for the dynamics of a single skyrmion, but we expect them to remain valid for a system with a low density of skyrmions [10]. We fix the obstacle density for the triangular (ρ_t = 0.072) and honeycomb (ρ_h = 0.123) arrays and vary α_m/α_d to investigate how the θ_sk = −30° and −60° steps evolve. In figure 12 we plot dynamic phase diagrams as a function of α_m/α_d versus F^D highlighting these selected locking steps. For the triangular array in figure 12(a), both steps have a range of α_m/α_d for which the step extends over nearly the entire window of F_D . This range is wider for the θ_sk = −60° step than for the θ_sk = −30° step. This could be of interest for technological applications since the direction of the skyrmion motion remains unchanged on each step. We find that when α_m/α_d ⩽ 1.2, the skyrmion only locks at θ_sk = −30°, while for α_m/α_d ⩾ 1.3, the skyrmion only locks with θ_sk = −60°. That is, the θ_sk = −30° and −60° steps do not coexist for a given value of α_m/α_d. This feature could be employed to perform topological selection, since skyrmions with a stronger Magnus component could lock to a different direction than skyrmions with weaker Magnus components. In the honeycomb array, figure 12(b) shows that the θ_sk = −30° and −60° steps coexist over the range 1.5 ⩽ α_m/α_d ⩽ 5.0. Here, the skyrmion can be switched from θ_sk = −30° motion to θ_sk = −60° motion by a fine adjustment of the external force. For α_m/α_d < 1.5 the skyrmion can only lock to θ_sk = −30°, while for α_m/α_d > 5.0, the skyrmion locks only to θ_sk = −60°. As mentioned in the previous section, the obstacle density plays an important role in determining the skyrmion dynamics, so we expect that changing the obstacle density may result in different phase diagrams.

Zoom In Zoom Out Reset image size

Using triangular or honeycomb obstacle arrays, the skyrmion motion can be locked into a range of different directions, with the strongest locking occurring for θ_sk = −30° and −60°. Both of these locking directions are robust against changing the obstacle density or the value of α_m/α_d. For a given value of the external applied drive, different obstacle densities can produce different directions of skyrmion motion. This property could be harnessed to build devices containing regions with different obstacle densities which can steer the skyrmion along a desired trajectory. We expect that similar results would appear for other periodic array geometries but that the angles of motion would differ depending on the substrate symmetry. The impact of thermal fluctuations would be an interesting issue to address in a future study. It is known that temperature can modify phase transition points, or even cause them to vanish [48]. Thermal creep could also become relevant [41]. If attractive pinning sites are used instead of repulsive obstacle sites, we expect very similar results, but have found that attractive pinning sites produce less pronounced locking steps [10]. Our results are based on a point-like model for skyrmions [23]; however, actual skyrmions have internal modes that can be excited which can modify the skyrmion dynamics. Additionally, the skyrmions may be distorted in shape by the driving force or through interactions with the obstacles [33]. Such effects could be explored using continuum-based simulations. Although we consider only a single skyrmion, we expect our results to be general for the case of multiple skyrmions at sufficiently low density, so that skyrmion–skyrmion interactions remain unimportant. If the skyrmion density is higher, it is possible that new dynamic phases will arise due to the collective behavior [54]. Lattice commensuration effects [55] could also become important. When the skyrmion lattice is commensurate with the pinning or obstacle array, the skyrmions can arrange themselves in such a way that the skyrmion–skyrmion interactions cancel out, producing dynamics that are similar to those found in the single skyrmion case.

We have investigated the dynamics of a single skyrmion driven over triangular and honeycomb obstacle arrays under zero temperature in order to determine the effects of changing the obstacle density and develop possible new ways to control the skyrmion motion. We show that the skyrmion exhibits a series of directional locking effects that can be quantized or continuous as a function of the applied drive. For low obstacle densities, the depinning forces are very weak and the skyrmion tends to move very close to the intrinsic Hall angle, reducing the number of dynamic phases. At higher obstacle densities, the depinning force is larger and a richer variety of directional locking phases appear. The main difference between the triangular and honeycomb obstacle arrays is the absence of a θ_sk = 0° locking step in the honeycomb array. The vacancies in the honeycomb lattice trap the skyrmions more effectively and prevent them from moving until the drive is too large to permit motion along the θ_sk = 0° direction. For the triangular array, θ_sk = 0° steps are more prominent at higher obstacle densities and lower values of α_m/α_d. Both arrays show pronounced locking at θ_sk = −30° and −60° due to the array symmetry. These locking steps appear over a wide range of obstacle densities and α_m/α_d values. We demonstrated the robustness of these phases for a variety of α_m/α_d values, and discussed possibilities for switching the skyrmion motion between the locking steps or using the locking effects to perform topological selection of different skyrmion species.

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 US Department of Energy (Contract No. 892333218NCA000001). NPV acknowledges funding from Fundação de Amparo à Pesquisa do Estado de São Paulo - FAPESP (Grant 2017/20976-3).

All data that support the findings of this study are included within the article (and any supplementary files).