Structural Transitions in Vortex Systems with Anisotropic Interactions

We introduce a model of vortices in type-II superconductors with a four-fold anisotropy in the vortex-vortex interaction potential. Using numerical simulations we show that the vortex lattice undergoes structural transitions as the anisotropy is increased, with a triangular lattice at low anisotropy, a rhombic intermediate state, and a square lattice for high anisotropy. In some cases we observe a multi-$q$ state consisting of an Archimedean tiling that combines square and triangular local ordering. At very high anisotropy, domains of vortex chain states appear. We discuss how this model can be generalized to higher order anisotropy as well as its applicability to other particle-based systems with anisotropic particle-particle interactions.


I. INTRODUCTION
A type-II superconductor subjected to a magnetic field forms vortices that each carry one quantum of magnetic flux. Due to their mutual repulsion, these vortices arrange themselves in an ordered vortex lattice (VL) the vortex-vortex interactions dominate external influences such as pinning by impurities or thermal disordering. The VL is triangular in isotropic superconductors 1 ; however, anisotropy in the vortex-vortex interactions can cause changes in the VL symmetry. 2 One of the most studied systems exhibiting a transition to a square VL is borocarbide materials [3][4][5][6][7][8][9][10] . Transitions to square lattices have also been observed in high temperature superconductors, [11][12][13][14] heavy fermion materials, [15][16][17][18] and spin triplet superconductors such as Sr 2 RuO 4 . [19][20][21] Additionally, transitions from triangular to square VLs can arise in other superfluid systems including Bose-Einstein condensates 22 and dense nuclear matter in extreme conditions. 23 In superconducting systems, theoretical approaches used to study the VL structural transitions include modifications to the London model 24,25 , addition of four-fold symmetric terms to the Ginzburg-Landau free energy 26 , Eilenberger theory 27 , modified Ginzburg-Landau approaches 28 , and modifications to the vortex interactions produced by strain fields 29 . Notably absent from this list are molecular dynamics (MD) simulations, which treat the vortices as point particles with bulk Bessel function interactions or thin film Pearl interactions. Until now, MD methods have only been applied to isotropic pairwise potentials, which produce a triangular VL in clean systems [30][31][32][33][34][35][36][37][38][39][40] . One of the issues is that simply adding a stronger radial repulsive force along certain directions does not capture the triangular to square transitions in the vortex system. Instead, it is essential to consider the full anisotropic potential in which nonradial forces also arise. In addition to the equilibrium VL configurations, MD simulations give access to the statics and dynamics of a large number of vortices over long times.
Here we introduce a model for vortices in anisotropic superconductors where the vortex-vortex interaction potential is extended to include a fourfold anisotropy. Using MD simulations we show that as the anisotropy is increased, the VL symmetry changes from triangular to rhombic to square. Additionally, in some cases we find a multi−q lattice state consisting of an Archimedean tiling with a combination of square and triangular local ordering. For very large anisotropy, domains of vortex chains appear. Our model can be generalized to any anisotropy and provides a basis for modeling VL reorientation or other vortex symmetry transitions such as those observed in multigap superconductors. 41 The model can be applied not only to vortices in type-II superconductors, but also to other particle-based systems with anisotropic interactions where triangular to square transitions can arise. These include skyrmion lattices, where triangular to square transitions have recently been observed, 42,43 as well as colloidal particles with anisotropic interactions. 44,45

II. MODEL
In an isotropic bulk superconductor, the vortex-vortex interaction potential is typically isotropically repulsive and takes the form of a zeroth order Bessel function, U (r) = K 0 (r). 46 For anisotropic materials, we propose the following modified pair potential for the vortex interactions: where r = |r i − r j | is the distance between two vortices at positions r i and r j . The angle between the two vortices with respect to the positive x axis is given by θ = tan −1 (r y /r x ) where r = r i − r j , r x = r ·x, and r y = r ·ŷ. The rotation angle of the anisotropic axes is defined to be φ, and n a is the order of the anisotropic interaction. The prefactor A v represents the strength of the isotropic component of the vortex-vortex pair interaction force, and A a controls the amplitude of the anisotropic interaction. To fully understand the manner in which Eq. (1) produces nonradial interactions, we examine the nature of the anisotropic forces generated by U (r). For nonzero values of A a , the potential is stronger along certain interaction directions and weaker between these directions. The anisotropic order n a determines the number of strong interaction directions, which are evenly spaced in θ. For example, for n a = 6 and φ = 0, the interaction po-tential passes through local maxima at θ = 0 • , 60 • , 120 • , 180 • , 240 • , and 300 • . Setting φ = 15 • , the local maxima of the interaction potential are shifted to θ = 15 • , 75 • , 135 • , 195 • , 255 • , and 315 • . Between these values of θ, the interaction strength varies as a cosine. The force field from the potential is In the present work we focus on the n a = 4 interaction potential with fourfold anisotropy; however, by changing n a our model can be applied to twofold (n a = 2), sixfold (n a = 6), or other degrees of anisotropy. In Fig. 1(a,b) we plot the equipotential lines and force field for U (r) from Eq. 1 for the isotropic case with A a = 0, where the equipotential lines are circularly symmetric and the forces are strictly radial. At A a = 0.1 in Fig. 1(c,d), the fourfold symmetry is apparent and weak nonradial forces appear. In Fig. 1(e,f) at A a = 0.25, the fourfold symmetry is much more pronounced and the nonradial forces are clearly visible.
To investigate the VL ground states that emerge as the anisotropy of the potential increases, we perform MD simulations of N = 441 vortices in a two-dimensional system of size L × L with L = 36λ and periodic boundary conditions in the x and y directions. Distances are measured in units of the London penetration depth λ. The dynamics of vortex i is governed by an overdamped equation of motion: Here η is the damping constant which we set equal to unity. Thermal forces are modeled by Langevin kicks where k B is the Boltzmann constant. We perform simulated annealing by starting in a high temperature molten state and gradually cooling the system to T = 0. We use an initial temperature of F T = 4.0 and decrement the temperature by ∆F T = −0.01 every 40,000 simulation time steps, which is long enough to ensure that the system reaches an equilibrium state.

III. VORTEX STRUCTURES
We use a Voronoi construction to obtain the local coordination number z i of each vortex, and compute the fraction P n of vortices with coordination number n using P n = N −1 N i=1 δ(z i − n) for n = 4, 5, 6, and 7. Figure 2(a) shows P 4 , P 5 , P 6 , and P 7 versus F T obtained during the annealing process in an isotropic system with A a = 0 and A v = 2.0. Initially the system is in a high temperature molten state, and as F T is reduced the vortices order into a triangular lattice with P 6 = 1.0. In Fig. 2(b) at A a = 0.039, the vortices freeze into a state with rhombic ordering, while in Fig. 2(c) at A a = 0.099, the vortices form a square lattice with P 4 = 1.0. Together, these results show how the equilibrium (T = 0) value of P 6 is suppressed and the value of P 4 grows as the anisotropy is increased. Compared to the drop in P 6 , the rise of P 4 is more gradual and happens at a lower F T .
We can further characterize the final F T = 0 state us- ing the structure factor S(k) Fig. 3(a) we illustrate the final real space vortex positions in a sample with A a = 0 and A v = 2.0, and in Fig. 3(b) we plot the corresponding |S(k)| as a heightfield. The lattice angle θ l is defined to be the angle between adjacent first-order peaks in |S(k)|, as illustrated in Fig. 3(b), and for the A a = 0 triangular lattice, θ l ≈ 60 • . The triangular lattice is slightly distorted by our perfectly square simulation box, which gives us the minor deviation from the ideal value θ l = 60 • . In Fig. 3(c,d), at A a = 0.039 the vortices form a rhombic lattice with θ l ≈ 76 • . At A a = 0.099 in Fig. 3(e,f), a square lattice appears with θ l = 90 • . In Fig. 4 we plot the lattice angle θ l versus A a for samples with A v = 1.0, 1.5, 2.0, and 2.5. In all cases, at low A a 0.3 the vortex ordering is triangular, at intermediate A a a rhombic lat- tice structure appears, and at large A a 0.5 a square VL emerges. In Fig. 5 we plot a structural phase diagram indicating where the triangular, rhombic, and square VLs appear as a function of A a versus A v . The evolution of the VL symmetry is nearly independent of the value of A v , showing that the triangular-to-square transition is a robust feature of our MD simulations.
For some combinations of A a and A v what we term a multi-q state appears in which the vortices exhibit simultaneous square and triangular ordering. Figure 6 shows both the real space Voronoi polygons and the corresponding |S(k)| for multi-q states in samples with A a = 0.04. For both A v = 2.0 in Fig. 6(a,b) and A v = 2.8 in Fig. 6(c,d) we find the same combination of square and triangular ordering in the Voronoi polygons, while |S(k)| exhibits multiple sets of peaks. As illustrated in Fig. 6, the multi-q ordering can be oriented along either the x or the y direction. The multi-q VL structure closely resembles an Archimedean tiling in which space is filled with a combination of square and triangular tiles. 47 Archimedean ordering of this type has also been observed for colloidal assemblies driven over quasiperiodic substrates, where it arises due to the competition between the ordering imposed by the substrate and the triangular ordering that minimizes the colloid-colloid interaction energy. 48,49 In our system the competition responsible for producing this structure is between the square and triangular orderings favored by the anisotropy. We note that the multi-q state only appears occasionally in regions of A a and A v that are dominated by the rhombic state, suggesting that it could be metastable. Additionally, the two different orientations of the multi-q state that we find indicate that in diffraction experiments on macroscopic samples, domains of different orientations will most likely coexist. This would smear out the peaks in |S(k)|, making it difficult to deconvolute the signal from an individual domain orientation. As a result, local probe imaging techniques may be the best method for observing multi-q states.
For anisotropies A a larger than those discussed above, we find that the square lattice gradually transforms into domains of vortex chains. This process is illustrated in

IV. SUMMARY
We have introduced a model for vortices with anisotropic pairwise interactions, focusing on the case of four-fold asymmetry. Using MD simulations we show that this model captures a transition from a triangular lattice at low anisotropy to a square VL at high anisotropy, with an intermediate rhombic phase. We also find that in some cases a multi-q state with Archimedean ordering appears in which the vortices have both square and triangular local ordering. For the highest anisotropy values, domains of vortex chains form. Our model could be applied to study the dynamics near the VL transitions, where nonequilibrium phenomena can arise. 53,54 Additionally, it can be generalized for higher order anisotropy in order to capture other types of symmetry and reorientation transitions in VLs. [55][56][57][58][59] . It is also possible to use the model with different isotropic pairwise interactions to investigate hexagonal to square transitions in other particle-based systems, such as skyrmions or colloids with anisotropic interactions.