Driven superconducting vortex dynamics in systems with twofold anisotropy in the presence of pinning

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 [2, 3], colloidal particles [4–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 enter a disordered or fluctuating state containing numerous topological defects [2, 15]. One of the most studied systems of this type is vortices in type-II superconductors, which can exhibit dynamical transitions at higher drives into ordered moving phases such as a moving crystal [16], anisotropic crystal [17], or moving smectic [18–21]. These transitions are associated with changes in the structure factor [17–21], the number and orientation of topological defects [18, 20, 21], and the noise characteristics [21–25]. In a two-dimensional (2D) system driven over random disorder, the fluctuations experienced by the particles due to their motion are anisotropic, leading to the formation of a moving smectic state [19]. Beyond superconducting vortices, moving smectics have been studied in other 2D systems driven over quenched disorder including Wigner crystals [14] and frictional systems [26].

In most of the above mentioned systems, the particle–particle interactions are isotropic, and 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 twofold anisotropic interactions, including colloidal particles in tilted magnetic fields [27, 28], dusty plasmas [29], electron liquid crystal states [30–32], skyrmions [33–35], and superconducting vortices [36–43]. Anisotropic vortex–vortex interactions can arise from anisotropy in the material or nematicity in the substrate, or it can be induced by a tilted field in uniaxial superconductors. Theoretical work on vortex liquid crystal systems with twofold 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]. 2D 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 along different anisotropy directions.

In this work we study the dynamics of superconducting vortices with a twofold anisotropy potential in the presence of pinning. Here driving forces are applied along the hard and soft directions of the vortex–vortex interaction as the magnitudes of the anisotropy and the pinning density are varied. The results may be applicable to members of the family of iron-based superconductors, where a strong electronic in-plane ansiotropy and elliptical vortices are observed [47]. However, while some studies of the anisotropic dynamics have been reported [48], achieving single domain samples and eliminating vortex pinning to twin boundaries remains a challenge. 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.

In previous computational modeling of vortices as point particles, the vortices had a pairwise isotropic repulsive potential that is proportional to the zeroth order Bessel function, U(r) = K₀(r) [49], causing the vortices to form a triangular ground state in the absence of quenched disorder. Some attempts were made to introduce effects of anisotropy by multiplying the force components of an isotropic interaction potential by different factors along orthogonal 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.

Recently, point particle models incorporating twofold, four-fold, and six-fold anisotropic interactions have been considered in the context of triangular to square and other vortex lattice rotational transitions [50, 51]. This work used a vortex–vortex interaction potential with n_a anisotropy axes,

$$\begin{equation}U(r,\theta )={A}_{\mathit{v}}{K}_{0}(r)\left[1+{K}_{\mathrm{a}}\,{\mathrm{cos}}^{2}\left(\frac{{n}_{\mathrm{a}}(\theta -{\phi }_{\mathrm{a}})}{2}\right)\right],\end{equation} \tag{ 1 }$$

where r = |r_i − r_j | is the distance between vortices at positions r_i and r_j and the angle between the vortices with respect to the x-axis is θ = tan⁻¹(r_y /r_x ) with r = r_i − r_j , ${r}_{x}=\mathbf{r}\cdot \hat{\mathbf{x}}$ and ${r}_{y}=\mathbf{r}\cdot \hat{\mathbf{y}}$. Here A_v is the isotropic vortex interaction strength that is used as a normalization parameter, as discussed below. The magnitude of the anisotropic contribution to the vortex interaction is given by K_a. The interaction potential in equation (1) produces a more complex interaction, including non-radial forces, between the particles than simple multiplicative factors, and hence provides a better realization of anisotropic systems.

In this work a twofold anisotropy (n_a = 2) is studied. Specifically, we examine the dynamics of systems with different anisotropy strengths driven over quenched disorder. Figure 1 illustrates the equipotential lines for the twofold anisotropic vortex potential with K₂ = 0, 0.5, 1.0 and 1.5. For K₂ = 0, the interaction is isotropic and the vortices form a triangular lattice with no preferred orientation. 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. Increasing or decreasing the value of K₂ changes the magnitude of the energy potential experienced by each vortex, which is equivalent to a change in the effective vortex density. This may be quantified by an effective magnetic field that is proportional to the 2D integral of the interaction potential:

$$\begin{equation}{B}_{\text{eff}}\propto {\int }_{0}^{\infty }r\,\mathrm{d}r{\int }_{0}^{2\pi }\mathrm{d}\theta \,U(r,\theta )\propto {A}_{\mathrm{v}}\left(2+{K}_{2}\right).\end{equation} \tag{ 2 }$$

Using K₂ = 0 and A_v = 2.0 as a reference, potential effects arising from a density difference are eliminated by setting A_v = 4/(2 + K₂) for each individual simulation, such that all have the same value of B_eff.

Zoom In Zoom Out Reset image size

Molecular dynamics simulations were carried out for a 2D system of size L × L with L = 72λ, where λ is the London penetration depth, with periodic boundary conditions imposed in both the x- and y-directions. For all simulations the number of vortices was N_v = 2280 corresponding to a vortex density ρ_v = N_v/L² = 0.4398/λ². The dynamics of vortex i are governed by the following overdamped equation of motion:

$$\begin{equation}\eta \frac{\mathrm{d}{\mathbf{r}}_{i}}{\mathrm{d}t}={\mathbf{F}}_{i}^{\mathit{v}\mathit{v}}+{\mathbf{F}}_{i}^{\mathrm{p}}+{\mathbf{F}}_{i}^{\mathrm{T}}.\end{equation} \tag{ 3 }$$

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

$$\begin{equation}{F}_{x}={A}_{\mathit{v}}\left[\mathrm{cos}(\theta ){K}_{1}(r)\left(1+{K}_{2}\,{\mathrm{cos}}^{2}(\theta )\right)-\frac{{K}_{2}}{r}{K}_{0}(r)\,\mathrm{sin}(\theta )\mathrm{sin}(2\theta )\right]\end{equation} \tag{ 4 }$$

$$\begin{equation}{F}_{y}={A}_{\mathit{v}}\left[\mathrm{s}\mathrm{i}\mathrm{n}(\theta ){K}_{1}(r)(1+{K}_{2}\;{\mathrm{c}\mathrm{o}\mathrm{s}}^{2}(\theta ))+\frac{{K}_{2}}{r}{K}_{0}(r)\;\mathrm{c}\mathrm{o}\mathrm{s}(\theta )\mathrm{s}\mathrm{i}\mathrm{n}(2\theta )\right].\end{equation} \tag{ 5 }$$

Each vortex also experiences forces ${\mathbf{F}}_{i}^{\mathrm{p}}$ from the substrate, which is modeled as N_p parabolic pinning traps placed in random but non-overlapping positions. Each pinning site has a 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 ${\mathbf{F}}_{i}^{\mathrm{p}}={F}_{\mathrm{p}}{\sum }_{j}^{{N}_{\mathrm{p}}}({\mathbf{r}}_{j}^{\mathrm{p}}-{\mathbf{r}}_{i})\,{\Theta}({r}_{\mathrm{p}}-\vert {\mathbf{r}}_{i}-{\mathbf{r}}_{j}^{\mathrm{p}}\vert )$. The individual pinning parameters are kept fixed at r_p = 0.5λ and F_p = 0.5, which allow a single pinning site to capture at most one vortex. Pinning densities ρ_p = N_p/L² ranging from zero to 0.4825/λ² were studied. Thermal forces are modeled by Langevin kicks ${\mathbf{F}}_{i}^{\mathrm{T}}$ with the properties ⟨F^T⟩ = 0.0 and $\langle {\mathbf{F}}_{i}^{\mathrm{T}}(t)\,{\mathbf{F}}_{j}^{\mathrm{T}}({t}^{\prime })\rangle =2\eta {k}_{\mathrm{B}}T{\delta }_{ij}\,\delta (t-{t}^{\prime })$ where k_B is the Boltzmann constant.

The initial vortex positions are obtained using simulated annealing. Starting from a high temperature F^T = 4.0 where the vortices are rapidly diffusing in a liquid state, the system is gradually cooled to zero temperature in decrements of ΔF^T = −0.05 every 10⁴ simulation time steps. After the initial annealing, a driving force F_D is applied 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. The drive starts at F_D = 0.0, increasing by ΔF_D = 0.05 every 10⁴ simulation time steps up to a maximum value of 1.5. The drive is then reduced back to zero at the same rate.

The structural properties of the vortex system are analyzed in real space 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}=\frac{1}{{N}_{\mathrm{v}}}{\sum }_{i=1}^{{N}_{\mathrm{v}}}\,\delta ({z}_{i}-n)$ for n = 5, 6, and 7. The most relevant parameter is the fraction of defects P_d = 1 − P₆, which provides a measure of the disorder of the system. All reported values of P_d were obtained by averaging the results of 10 simulation runs. The structure factor, $S(\mathbf{k})=\frac{1}{{N}_{\mathrm{v}}}\left\vert {\sum }_{i=1}^{{N}_{\mathrm{v}}}{\text{e}}^{\text{i}\mathbf{k}\cdot {\mathbf{r}}_{i}}\right\vert$, provides a complementary measure of the order in reciprocal space.

Figure 2(a) shows the fraction of topological defects P_d versus driving force F_D for a system with K₂ = 0.55 and ρ_p = 0.4825/λ². Initially 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.25 that coincides with a maximum in the defect density of P_d = 0.45 for y-direction driving and P_d = 0.48 for x-direction driving. Generally, the depinning threshold is slightly higher for driving in the x-direction due to the relatively larger average spacing between the vortices in that direction induced by the anisotropy. Above depinning, P_d decreases rapidly 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₂. As the drive is decreased from its maximum value of F_D = 1.5, the system returns to a disordered state for both directions. Increasing the twofold anisotropy to K₂ = 1.45 while keeping the pinning density constant leads to a distinction between drives along the x- and y-directions, as shown in figure 2(b). Here, there is a clear difference in the defect density, with y-direction driving producing much lower values of P_d, indicating that the system is better ordered when the drive is aligned with the soft anisotropy direction corresponding to the natural smectic orientation of the system. For driving in the x and y-directions, the defect fraction reaches a minimum value of P_d = 0.25 and 0.14 respectively Again, the system disorders as the drive is reduced, reaching nearly identical values of P_d ≈ 0.34 at the pinning transition.

Zoom In Zoom Out Reset image size

To complement the characterization in terms of the overall defect fraction in figure 2, the corresponding structural properties in both real and reciprocal space were examined. Figure 3(a) shows a Voronoi construction of a portion of the system from figure 2(a) for K₂ = 0.55, ρ_p = 0.4825/λ², and zero drive, where there are a large number of dislocations indicated by non-sixfold coordinated vortices. The corresponding structure factor S(k) in figure 3(d) contains a set of broad peaks along the vertical direction for k_x = 0, indicative of weak nematic ordering with one-dimensional (1D) chains of vortices aligned in the y-direction. The anisotropy of the vortex–vortex interaction is evident from the elliptical diffuse intensity, since an isotropic interaction would give rise to a circle. Figures 3(b) and (e) show the same system with a drive of F_D = 1.5 applied along the x-direction. Here, only a small number of vertically aligned dislocations are present. The corresponding structure factor is still anisotropic but now has sharp peaks in the y-direction indicative of a smectic phase. A similar moving smectic phase is seen in figures 3(c) and (f) when a drive of F_D = 1.5 is applied along the y-direction. In contrast to earlier work for isotropic systems where the smectic phase is always oriented along the direction of the drive [18, 20, 21], we observe an alignment along the soft anisotropy direction for K₂ > 0.2 regardless of the direction of the drive.

Zoom In Zoom Out Reset image size

Figure 4 summarizes the main structural properties observed for the higher anisotropy of K₂ = 1.45. A pinned nematic state still appears at zero drive as shown in figures 4(a) and (d), with numerous defects indicated by five and sevenfold coordinated vortices and broad peaks on the vertical axis in the structure factor for k_x = 0 together with diffuse lines on the left and right. Figures 4(b) and (e) show the same system for driving in the y-direction with F_D = 0.75. The vortices are more ordered, with fewer defects in the Voronoi construction and sharper peaks in S(k), and the system has formed a moving nematic. In figures 4(c) and (f), for an increased drive of F_D = 1.5 in the y-direction, the number of defects has diminished further and the peaks in S(k) are still sharper. The vortices are now organized in a series of non-overlapping 1D chains indicating smectic ordering.

Zoom In Zoom Out Reset image size

The difference between nematic and smectic phases are illustrated in more detail in figure 5. Figure 5(a) shows a smaller region of the Voronoi construction from figure 4(a), and figure 5(b) contains a corresponding schematic of the vortex configuration. Here the system forms chains that begin, end and intertwine inside the sample, forming a nematic structure. In figure 5(c) we show a small region of the Voronoi construction from figure 4(c) for y-direction driving with F_D = 1.5 where the system forms a smectic state and the 1D chains do not overlap. Figure 5(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. The lack of order perpendicular to the chains is evident from the structure factor in figure 4(f), which has diffuse vertical lines, rather than peaks, for k_x ≠ 0. Changing the drive to be along the x-direction for the K₂ = 1.45 system produces a very similar result, although the nematic phase persists up to higher values of F_D. Finally, we note that the transition between nematic and smectic states is gradual, and near the cross-over the classification can be ambiguous as seen in figures 4(b) and (e).

Zoom In Zoom Out Reset image size

To characterize the system across the entire range of anisotropies and drives, one can consider the variation in the defect density during the cycle F_D = 0 → 1.5 → 0. Figure 6(a) and (b) show the minimum and maximum values of P_d as a function of K₂ in samples with ρ_p = 0.4825/λ² for drives along the x-direction and the y-direction, respectively. The maximum falls at the depinning threshold and the minimum is at the highest drive of F_D = 1.5. Consider first the results in figure 6(a) for driving in the y-direction. For low K₂ < 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 in the range P_d = 0.03–0.07. For K₂ > 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. Turning now to the results in figure 6(b) for driving in the x-direction, one again observes a smaller number of defects at depinning with increasing K₂ due to the greater durability of the chains, but a greater number of defects appear at the maximum drive due to the preservation within the chains of the defects that were present at low drives. Furthermore, there is a peak in the minimum number of defects near K₂ = 0.2, where the smectic undergoes a transition from alignment in the x-direction to alignment in the y-direction. This reorientation transition arises due to the competition between the twofold anisotropy, which favors an alignment along the soft direction of the vortex–vortex interaction, and the motion of the vortices, which favors an alignment along the direction of the drive [18, 20, 21]. Thus for K₂ < 0.2, the system forms a moving smectic aligned in the x-direction, while for 0.2 < K₂ < 0.8, the system forms a moving smectic aligned in the y-direction. When K₂ > 0.9, a moving nematic aligned in the y-direction appears at higher drives.

Zoom In Zoom Out Reset image size

Figure 7 shows the Voronoi construction and structure factor for the system at the critical value K₂ = 0.2 for a drive of F_D = 1.5 in the x-direction. Here, the system does not form a nematic or smectic aligned in the y-direction, but instead adopts a polycrystalline ordering with domain boundaries decorated by strings of defects. One of the domains is oriented along the x-direction, as evident from the structure factor that has two prominent peaks on the horizontal axis for k_y = 0. For even smaller values of K₂, the system exhibits a strong smectic alignment along the x-direction for driving in the x-direction.

Zoom In Zoom Out Reset image size

The fraction of defects shows a hysteretic behavior for the cycle F_D = 0 → 1.5 → 0. This is quantified by numerically integrating the magnitude of the difference in the defect fraction for increasing $({P}_{\mathrm{d}}^{+})$ and decreasing $({P}_{\mathrm{d}}^{-})$ drive $H={\int }_{0}^{1.5}\vert {P}_{\mathrm{d}}^{+}({F}_{\mathrm{D}})-{P}_{\mathrm{d}}^{-}({F}_{\mathrm{D}})\vert \,\mathrm{d}{F}_{\mathrm{D}}$, indicated by the shaded areas in figure 8. Figure 9 shows the hysteresis for two different values of the pinning density which display qualitatively similar features. First, the magnitude of H increases with pinning density, with the largest values occurring for lower values of K₂ where P_d varies over a greater range and a nearly perfect lattice at high drives remains more robust against pinning forces as F_D is decreased. Second, for driving in the x-direction, there is a peak in H near K₂ = 0.2 corresponding to the smectic x-direction to y-direction realignment transition. Third, a crossover is observed at K₂ ≈ 1.0 [figure 9(a)] and K₂ ≈ 0.9 [figure 9(b)] for the higher and lower pinning density respectively. Below the crossover, the hysteresis is largest for driving in the x-direction, while above the crossover, the situation is reversed.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

A broader survey of pinning densities up to ρ_p = 0.4825/λ² and different directions of the drive is given in figure 10, which shows heat maps of the minimum value of P_d (i.e., for F_D = 1.5) as a function of K₂ versus ρ_p. For driving in the x-direction, the light blue line at K₂ ≈ 0.2 in figure 10(a) indicates the switching of the smectic from x-direction to y-direction alignment. The critical value of K₂ at which this transition occurs increases with increasing ρ_p, since a greater anisotropy is required to form chains aligned in the soft anisotropy direction that are also robust against being driven over a higher density of pinning sites. For K₂ > 1.0 and ρ_p > 0.1, the amount of disorder in the system increases dramatically and a nematic state replaces the smectic state. For driving in the y-direction, shown in figure 10(b), the smectic structures formed by the vortices are always aligned in the y-direction. Furthermore, the system generally remains in a smectic state even as the number of defects increases for K₂ > 1.0 and ρ_p > 0.19.

Zoom In Zoom Out Reset image size

We have examined the driven dynamics of vortices with twofold 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. Cycling the drive leads to a 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 as a function of anisotropy and pinning density, characterized by the defect density and the structure factor. Our results should be general to the broader class of driven systems with twofold anisotropy driven over random disorder, which includes electronic liquid crystals, colloidal particles, and magnetic skyrmions.

We are grateful to D Minogue, M W Olszewski and D Spulber for assistance with the molecular dynamics simulations and analysis. This research was supported in part by the Notre Dame Center for Research Computing. Work at the University of Notre Dame (EJR, MRE: MD simulations, data analysis) was supported by the US Department of Energy, Office of Basic Energy Sciences, under Award No. DE-SC0005051. Part of this work (CR, CJOR: code development) was carried out under support 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).

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