Driven Superconducting Vortex Dynamics in Systems with Two-Fold Anisotropy in the Presence of Pinning

We examine the dynamics of superconducting vortices with two-fold anisotropic interaction potentials driven over random pinning and compare the behavior under drives applied parallel and perpendicular to the anisotropy direction. The number of topological defects reaches a maximum near depinning and then drops with increasing driving force as the vortices form one-dimensional chains. This coincides with a transition from a pinned nematic to a moving smectic aligned with the soft direction of the anisotropy. The system is generally more ordered when the drive is applied along the soft direction of the anisotropy, while for driving along the hard direction, there is a critical value of the anisotropy above which the system remains aligned with the soft direction. We also observe hysteresis in the dynamics, with one-dimensional aligned chains persisting during a decreasing drive sweep to drives below the threshold for chain formation during the increasing drive sweep. More anisotropic systems have a greater amount of structural disorder in the moving state. For lower anisotropy, the system forms a moving smectic-A state, while at higher anisotropy, a moving nematic state appears instead.


Introduction
A wide range of systems can be described effectively as an assembly of particles interacting with each other and with quenched disorder, leading to the appearance of depinning and multiple sliding phases [1,2]. Such systems include vortices in type-II superconductors [3,2], colloidal particles [4,5,6], active matter [7,8], magnetic skyrmions [9,10], pattern forming systems [11,12], and sliding Wigner crystals [13,14]. In each case there is a threshold for motion or a depinning transition above which the particles can depin into a disordered or fluctuating state containing numerous topological defects [15,2]. One of the most studied systems of this type is vortices in type-II superconductors, which can exhibit dynamical transitions into more ordered moving phases such as a moving crystal [16], anisotropic crystal [17], or moving smectic [18,19,20,21]. These transitions are associated with changes in the structure factor [17,18,19,20,21], the number and orientation of topological defects [18,20,21], and the noise characteristics [22,21,23,24,25]. In a two-dimensional (2D) system driven over random disorder, the fluctuations experienced by the particle due to its motion over the quenched disorder are anisotropic, leading to the formation of a moving smectic state [19]. Beyond superconducting vortices, moving smectics have also been studied in other 2D systems driven over quenched disorder including Wigner crystals [14] and frictional systems [26].
In most of these systems, the particle-particle interactions are isotropic, so in the absence of quenched disorder an isotropic crystal appears. For example, in the case of type-II superconductors with isotropic repulsion, the vortices form a triangular lattice [3]. There are, however, many examples of particle-like systems that have two-fold anisotropic interactions, including colloidal particles in tilted magnetic fields [27,28], dusty plasmas [29], electron liquid crystal states [30,31,32], skyrmions [33,34,35], and superconducting vortices [36,37,38,39,40,41,42,43]. Anisotropic vortexvortex interactions can arise from anisotropy in the material or nematicity in the substrate, or it can be induced by a tilted field. Theoretical work on vortex liquid crystal systems with two-fold anisotropy showed that these systems can form smectic-A states and exhibit two step melting transitions [39,41]. Magnetic skyrmions have many similarities to superconducting vortices and typically form a triangular lattice under isotropic conditions [44,45]. Two-dimensional anisotropic skyrmions can produce what are called skyrmion liquid crystals with smectic [35] or more specifically smectic-A ordering [46]. Far less is known about the behavior of driven states of anisotropic crystals in quenched disorder, such as what dynamical ordering transitions appear and what differences arise when the driving is applied parallel or perpendicular to the anisotropy direction. In this work we study the dynamics of superconducting vortices with a twofold anisotropy potential that are driven over quenched disorder both parallel and perpendicular to the anisotropy for varied values of the anisotropy and the quenched disorder strength. In addition to superconducting vortices, our results should also be relevant to the wider class of assemblies of particles with twofold anisotropic interactions moving over quenched disorder.

Methods
In previous computational modeling of vortices as point particles interacting with pinning the vortices had a pairwise isotropic repulsive potential that is proportional to the zeroth order Bessel function, U (r) = K 0 (r) [47], causing the vortices to form Equipotential lines for the vortex-vortex interaction in Eq. (1) for increasing anisotropy parameter K 2 = 0.0 (a), 0.5 (b), 1.0 (c), and 1.5 (d). The deep pinch points along the y-axis cause chain states to form that affect the dynamics of the system. a triangular ground state in the absence of quenched disorder. Point particle models of vortices with two-fold, four-fold, and six-fold anisotropic interactions have also been considered [48,49] in the context of triangular to square and other rotational transitions. In anisotropic systems with a total of n a anisotropy axes, the vortexvortex interaction potential has the form where r = |r i − r j | is the distance between two vortices at positions r i and r j and the angle between the vortices with respect to the x-axis is θ = tan −1 (r y /r x ) with r = r i − r j , r x = r ·x and r y = r ·ŷ. Here A v is an isotropic vortex interaction strength that we use as a normalization parameter. The magnitude of the anisotropic contribution to the vortex interaction is given by K a . In this work we consider twofold anisotropy (n a = 2) or K 2 . Specifically, we examine the dynamics of systems with different anisotropy strengths driven over quenched disorder. In Fig. 1 we illustrate the equipotential lines for the two-fold anisotropic vortex potential with K a = 0.0, 0.5, 1.0, and 1.5. For K 2 = 0.0 in Fig. 1(a), the interaction is isotropic and the vortices form a triangular lattice. As the anisotropy increases, the potentials become more elongated, implying that the vortex-vortex interaction forces are strongest along the x-direction and weaker along the y-direction. In previous work, anisotropy was introduced by multiplying the force components of an isotropic interaction potential by different factors in the x-and y-directions [41]. Although this approach produced anisotropic diffusion, smectic ordering appeared only for very large differences in the multiplication prefactors, making this an unrealistic representation of anisotropic systems. A much better realization is the 2D anisotropic potential of the type shown in Fig. 1, which produces a much more complicated force configuration than a simple multiplication factor would.
We consider a two-dimensional system of size L × L with L = 72λ, where λ is the London penetration depth. The dynamics of vortex i are governed by the following overdamped equation of motion: Here η is the damping constant, which is set to unity. The vortex-vortex interaction force is F vv = −∇(U ) = (−∂U/∂x, −∂U/∂y). With a twofold anisotropy and φ a = 0, the force is There are a total of N v vortices in the sample. Each vortex also experiences forces F p i from the substrate, which is modeled as N p parabolic pinning traps placed in random but non-overlapping positions. Each pinning site is of radius r p and can exert a maximum force of F p . The vortex-pin interaction is directed toward the center of the pinning site and is given by ij . For the parameters we consider, an individual pinning site can capture at most one vortex. In this work we fix r p = 0.5λ and F p = 0.5.
The initial vortex positions are obtained using simulated annealing. Starting from a high temperature where the vortices are in a liquid state, we lower the temperature to zero in a series of steps. Thermal forces arise from Langevin kicks where k B is the Boltzmann constant. We begin the annealing process at F T = 4.0 where the vortices are rapidly diffusing, and gradually cool the system to F T = 0.0. The temperature is reduced by ∆F T = −0.05 every 10 4 simulation time steps.
After annealing we apply a driving force in either the x-or y-direction. Here x is the hard direction of the anisotropy along which the vortices are more repulsive, and y is the soft direction. We start at F D = 0.0 and increase the force in increments of ∆F D = 0.05 every 10 4 simulation time steps up to a maximum drive of F max = 1.5. We then decrease the drive by ∆F D every 10 4 time steps until F D = 0.0 again. We use a pinning density ρ p = N p /L 2 ranging from ρ = 0 to ρ = 0.48225 and fix the vortex density ρ v = N v /L 2 to ρ v = 0.44. For all results discussed here, values are obtained by averaging the results of 10 simulations for each set of parameters.
Increasing or decreasing the value of K 2 changes the magnitude of the energy potential experienced by each vortex, which is equivalent to a change in the effective vortex density. To eliminate effects arising from a density difference, we define an effective magnetic field that is proportional to the two-dimensional integral of the interaction potential: Using K 2 = 0 and A v = 2.0 as a reference, we set A v = 4/(2 + K 2 ) for each individual molecular dynamics simulation, such that all simulations have the same value of B eff . The vortex lattice configurations are analyzed after annealing and during the drive sweep processes by calculating the structure factor, S(k), and by using a Voronoi polygon construction. This yields the local coordination number z i of each vortex, which is used to compute the fractions P n = 1 Nv Nv i=1 δ(z i − n) for n = 5, 6, and 7. The most useful parameter is the fraction of defects P d = 1 − P 6 , which provides a measure of the disorder of the system.

Results
In Fig. 2 we plot the fraction of topological defects P d versus driving force F D for a system with K 2 = 0.55 and ρ p = 0.48225 under driving in the x-and y-directions. The system is in a disordered configuration after the annealing process with a large fraction, P d = 0.43, of vortices that are not sixfold coordinated. As F D increases, a depinning transition occurs near F D = 0.2 that coincides with a maximum in P d of P d = 0.45 for y-direction driving and P d = 0.47 for x-direction driving. The depinning threshold is slightly higher for driving in the x-direction. Above depinning, P d rapidly drops with increasing F D and approaches P d = 0.05 for F D > 1.0. The minimum value of P d is similar for driving in either direction at this value of K 2 . For driving in the y-direction, the decrease in P d as F D increases is correlated with dynamical ordering into a moving smectic phase containing a small number of dislocations that are aligned in the driving direction, as observed previously [18,20,21]. One distinction we find in the anisotropic system is that the smectic state for driving in the x-direction is not aligned with the driving direction but is instead aligned with the soft anisotropy or y-direction.
In Fig. 3(a) we show a Voronoi construction of a portion of the system from Fig. 2 with K 2 = 0.55 and ρ p = 0.48225 at F D = 0.0, where one-dimensional (1D) chains of vortices appear that are aligned in the y-direction. The corresponding structure factor S(k) in Fig 3(d) contains a set of diffusive peaks aligned in the k y direction along k x = 0, indicative of nematic ordering. Here the vortices are spaced more closely in the y-direction than in the x-direction due to the anisotropy of the repulsive vortex-vortex force, which is smaller along the y direction, permitting the vortices to approach each other more closely from this direction. In Fig. 3(c,f) when a drive of F D = 1.5 is applied along the y-direction, there are only a small number of dislocations present that are all aligned in the y-direction. The corresponding structure factor is still anisotropic but has sharp peaks in the y-direction indicative of a smectic phase. Figure 3(b,e) shows the same system with a drive of F D = 1.5 applied along the x-direction. A similar moving smectic appears that is perpendicular to the drive direction.
In Fig. 4 we plot P d versus F D for driving in the x− and y-directions for a sample with ρ p = 0.48225 at a larger anisotropy of K 2 = 1.45. There is a clear difference in the defect density, with y-direction driving producing much lower values of P d than x-direction driving over the range of drives we consider, indicating that the system is better ordered when the drive is aligned with the soft anisotropy direction or the natural smectic orientation of the system. For driving in the x-direction, P d reaches a minimum value of P d = 0.24, while we find a lower minimum of P d = 0.14 for driving in the y-direction. As the drive is decreased from its maximum value of F D = 1.5, the system starts to disorder again for driving in both directions, reaching nearly identical values of P d ≈ 0.32 at the pinning transition.
In Fig. 5(a,d) we show a Voronoi construction and structure factor for the system Figure 4. P d versus F D in a system with K 2 = 1.45 and ρp = 0.48225 for xdirection driving (blue) and y-direction driving (red). The system is more ordered for driving in the y-direction.
in Fig. 4 with K 2 = 1.45 at F D = 0, where a pinned nematic phase appears. Figure 5(b,e) shows the same system for driving in the y-direction with F D = 0.75. The vortices are more ordered, as indicated by the sharper peaks in S(k), and the system has formed a moving nematic. In Fig. 5(c,f), for driving in the y-direction at F D = 1.5, the number of defects has diminished and the peaks in S(k) are sharper. The vortices exhibit smectic ordering and form a series of non-overlapping 1D chains. In Fig. 6(a) we show a smaller region of the Voronoi construction from Fig. 5(a) at F D = 0.0 indicating that the vortices form chains that can break or intertwine. Figure 6(b) shows a schematic of the vortices with the anisotropic potentials forming a nematic structure. Here the system forms chains that can end or begin inside the sample. In Fig. 6(c) we show a small region of the Voronoi construction from Fig. 5(e) for y-direction driving with F D = 1.5 where the system forms a smectic state and the 1D chains do not overlap. Figure 6(d) shows a schematic of the vortex structure in this state, which is known as smectic-A in liquid crystal systems [46]. This is similar to the phase proposed for vortex liquid crystals with anisotropic potentials [39]. Here there are no breaks in the 1D chains. Individual chains can contain different numbers of vortices, producing dislocations that are aligned in the y-direction. In the K 2 = 1.45 system, driving in the x-direction produces a set of phases very similar to those found for driving in the y-direction, but the nematic phase persists up to higher drives.
We can also characterize the system by measuring the maximum and minimum number of defects generated during the drive cycle F D = 0 → 1.5 → 0 for varied K 2 . The maximum number of defects appear at the depinning threshold, while the minimum number are present at the highest drive of F D = 1.5. In Fig. 7(a) we plot the minimum and maximum values of P d versus K 2 for driving in the y-direction in samples with ρ p = 0.48225. For low K 2 < 0.8, at depinning there is a nematic state with P d ≈ 0.45. The vortices order into a moving smectic phase at higher drives with the defect density reaching minimum values of P d = 0.03 to P d = 0.07. For K 2 > 0.9, the system is less defected at the depinning transition but contains more defects in the driven reordered states. As the anisotropy increases, it becomes more difficult to destroy the chains of vortices at depinning, giving a lower density of defects at the depinning transition; however, it becomes easier for the chains to slide past one another at higher drives, creating a larger number of more persistent dislocations in the driven phase. In Fig. 7(b) we show the minimum and maximum values of P d versus K 2 for the same system under driving in the x-direction. Near K 2 = 0.2 there is a peak in the minimum number of defects corresponding to the critical anisotropy at which the smectic undergoes a transition from alignment in the x-direction to alignment in the y-direction. For K 2 > 0.2 the system forms a moving smectic aligned in the y-direction, while when K 2 > 0.9, a moving nematic appears at higher drives.
In Fig. 8 we show the Voronoi construction and structure factor for the system in Fig. 7(b) with K 2 = 0.2 for driving in the x-direction at F D = 1.5. Here, the system  does not form a nematic or smectic aligned in the y-direction, but instead adopts a polycrystalline ordering with a partial alignment in the x-direction. The ordering is more clearly visible in the structure factor, where two prominent peaks are aligned with k x . For even smaller values of K 2 , the system exhibits a strong smectic alignment along the x-direction for driving in the x-direction.
We can also characterize the system by measuring the total hysteresis in the form of the sum of the absolute differences in P d for a given drive F D , where we compare the value P + d for increasing current with P − d for decreasing current. The total hysteresis    is obtained by numerically integrating the absolute difference between the two curves, H = 1.5 as indicated by the shaded area in Fig. 9. We find that H is largest at lower values of the anisotropy, since in these systems the number of defects varies over a greater range and a nearly perfect lattice appears at high driving that remains more robust against pinning forces as the current is decreased. In Fig. 10(a) we plot the total hysteresis H versus K 2 for a system with ρ p = 0.48225 under driving in the x-and y-directions. For driving in the x-direction, there is a peak in H near K 2 = 0.2 corresponding to the smectic x-direction to y-direction realignment transition. When K 2 < 1.0, the hysteresis is largest for driving in the x-direction, but for larger K 2 , we find the largest hysteresis for driving in the y-direction. We observe similar behavior as we vary the pinning site density ρ p , with the overall magnitude of H gradually decreasing with decreasing ρ p . This is seen in Fig. 10(b), which shows H versus K 2 at a lower pinning density of ρ p = 0.09645, where the total hysteresis is lower but the same trends appear.
We have also considered the effect of changing the pinning density over the range ρ p = 0 to ρ p = 0.48225. In Fig. 11(a) we plot a heat map of the minimum value of P d (i.e., for F D = 1.5) as a function of K 2 versus ρ p for driving in the x-direction. The light blue line at K 2 ≈ 0.2 indicates the switching of the smectic from the x-direction to the y-direction alignment. As ρ p increases, the critical value of K 2 at which this transition occurs increases. For K 2 > 1.0 and ρ p > 0.1, the amount of disorder in the system increases dramatically and a nematic structure is present, while in other regions of the parameter space, we find a smectic state aligned in either the x-or ydirections. Figure 11(b) shows the heat map of the minimum value of P d for driving in the y-direction. In this case, the smectic structures formed by the vortices are always aligned in the y-direction. For K 2 > 1.0 and ρ p > 0.19, there is an increasing number of defects in the moving phase, but generally the system is in a smectic state.

Summary
We have examined the driven dynamics of vortices with two-fold anisotropic interaction potentials driven over quenched disorder. In general, the pinned states have nematic ordering and the driven states form moving smectic A phases. A more ordered smectic state containing fewer dislocations appears for driving along the soft anisotropy direction compared with driving along the hard anisotropy direction. When we cycle the drive, we observe hysteresis in the dynamics. Specifically, once the smectic state has formed, it can persist down to lower drives than those at which it appeared on the initial application. We also find that as the anisotropy increases, the system is generally less ordered in the high drive states since dislocations have a lower formation energy. We map out the dynamic phase diagram for this system as a function of varied anisotropy and disorder strength using as our characterization tools the orientation of dislocations and the features in the structure factor. Our results should be general to the broader class of driven systems with two-fold anisotropy driven over random disorder, which includes electronic liquid crystals, colloidal particles, and magnetic skyrmions.