{'fraction_non_alphanumeric': 0.029809795706499574, 'fraction_numerical': 0.016573356568165807, 'mean_word_length': 4.608650447078395, 'pattern_counts': {'":': 0, '<': 1, '<?xml version=': 0, '>': 3, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 2, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Activated vortex lattice transition in a superconductor with combined sixfold and twelvefold anisotropic interactions

Numerical simulations are used to examine the transition dynamics between metastable and ground state vortex lattice phases in a system with combined sixfold and twelvefold contributions to the vortex-vortex interactions. The system is initially annealed using a a twelvefold anisotropy, yielding domain of two different orientations and separated by grain boundaries. The vortex-vortex interaction is then suddenly changed to a sixfold anisotropy, rendering the twelvefold state metastable. Applying a drive that mimics an oscillating magnetic field causes the metastable state to decay, indicated by the structure factor that evolves from twelve to six peaks. The results fits the behavior seen in recent small-angle neutron scattering studies of the vortex lattice in MgB2. At higher drive amplitudes, the decay exhibits a two step process in which the initial fast decrease is followed by a slower regime where avalanches or bursts are correlated with defect annihilation events. The results are compared to other types of metastable systems with quenched disorder that decay under a periodic external drive.


I. INTRODUCTION
Magnetic flux enters a type-II superconductor in the form of quantized vortices 1 . In the absence of strong pinning by impurities or thermally driven disordering, the mutual repulsion between vortices causes them to crystallize into an ordered vortex lattice (VL) 2 . In a perfectly isotropic superconductor the VL would be triangular and without a preferred orientation relative to the host lattice 3 ; however, real materials always posses some degree of anisotropy that causes an alignment of the VL relative to the crystalline axes, often giving rise to field and/or temperature driven symmetry transitions 4 . Such transitions lead to the formation of competing VL domains with different orientations. The VL can thus serve as a simple model system for domain nucleation and growth, with relevance for a wide range of physical systems including martensitic phase transitions 5 , domain switching in ferroelectrics 6 , and the skyrmion lattice in chiral magnets [7][8][9] .
For magnetic fields applied along the c axis of superconductors with a hexagonal crystal structure, the VL is empirically found to always have a triangular symmetry, but will often undergo continuous and sometimes even non-monotonic rotations as seen in MgB 2 10 and UPt 3 11-13 . In addition, the VLs in these materials exhibit significant metastability when cooled or heated across equilibrium phase transitions in a constant field 11,14 . To reach the equilibrium configuration the VL must be excited by an external perturbation, which causes vortex motion and allows the system to find the global minimum in the free energy 15 . In MgB 2 the transition from the metastable to the equilibrium state has been studied extensively using small-angle neutron scattering (SANS), revealing an activation barrier that increases as the metastable state is gradually sup-pressed 16,17 . This "work hardening" is speculated to arise from a proliferation of grain boundaries in the VL. Although spatially resolved SANS measurements provide support for such a scenario 18 , a real space study of the VL is desirable.
Spatial information about the VL can be obtained by molecular dynamics (MD) methods in which the vortices are treated as point particles with repulsive interactions. In addition to static VL configurations, MD simulations provide information about the dynamics of a large number of vortices over long times. This technique has been widely used to study how pinning affects the structural 19 and dynamical [20][21][22] properties of vortex matter or its thermal melting 23 . Due, however, to the added complexity introduced by a fully anisotropic interaction potential and the resulting nonradial forces, such studies were limited to systems with isotropic vortex-vortex interactions. Recently we implemented a phenomenological model that allows MD simulations of the VL in the presence of anisotropic interactions [24][25][26] . This was used to study systems with combined sixfold and twelvefold anisotropies relevant for MgB 2 25 . By varying the ratio of the six-and twelvefold contributions to the interaction potential, it was possible to reproduce the experimentally observed VL phases in this material 10,14 . Specifically, F and I phases are oriented along the two high symmetry directions within the basal plane (a and a * ), and connected by an intermediate, continuously rotating L phase 10,27 . The simulations also spontaneously produced VL dislocations as well as a range of vacancy and interstitial defects that were analyzed in terms of their energy cost 25 .
In this work we extend our MD simulations to study the gradual VL transition from a metastable L phase towards the equilibrium F phase. The metastable L phase is obtained by first annealing the system with a vortex-vortex interaction potential that has a twelvefold anisotropy and subsequently changing this to a sixfold anisotropy. Weak pinning is incuded to stabilize the initial twelvefold symmetric phase, and the temperature is kept low enough that thermal fluctuations alone are insufficient to overcome the activation barrier for a transition to the ground state. In the SANS experiments, the transition to the ground state is driven by magnetic small-amplitude field oscillations, leading to a varying vortex density 15 . Such a scenario is not straightforward to implement in the MD simulations, due to the periodic boundary conditions, and we instead employ a periodic variation of the superconducting penetration depth. This changes the effective vortex-vortex interaction strength and thus mimics a field oscillation. As successive penetration depth cycles are applied, the metastable L phase gradually relaxes towards the ground state F phase. Initially, the angle between the two VL domain orientations decreases as a power law with an exponent that increases with the amplitude of the penetration depth variation, similar to what was observed in the SANS experiments 17 . For the larger amplitudes there is a cross-over to a logarithmic behavior as the number of penetration depth cycles increases. Here, the VL domain rotation occur in bursts or avalanches that are associated with correlated topological annihilation events along the grain boundaries. We discus how these behaviors compare to twodimensional (2D) coarsening models [28][29][30][31] , and show that the log(t) decay could arise due to coarsening or glassy motion with quenched disorder 32 .
Our approach can be applied to other particlebased systems with anisotropic interactions, including skyrmion lattices [33][34][35] or colloidal particles with anisotropic interactions 36,37 . Furthermore, it will allow modeling of the dynamics associated with the generation and recombination of dislocations, grain boundaries and domain formation both in the VL [16][17][18] and in general 38,39 . We discuss the relevance of our results to the decay of metastable states in other systems with periodic perturbations, such as granular matter 40,41 , colloidal systems 32 , magnetic systems with return point memory 42,43 , and other types of periodically sheared systems that relax to a steady state under ac driving [44][45][46][47][48] .

II. METHODS
We follow the approach of Olszewski et al. and use a phenomenological vortex-vortex interaction potential that includes six-and twelvefold anisotropies 25 : Here r and θ are the distance and angle, respectively, between two vortices at positions r i and r j , namely r = |r i − r j | and θ = tan −1 (r y /r x ), with r = r i − r j , r x = r ·x, and r y = r ·ŷ. The radial part of the repulsive interaction potential is given by the amplitude A v times the zeroth order Bessel function K 0 (r) 49 . The parameters K 6 and K 12 are the magnitudes of the sixand twelvefold contributions to the potential, with an anisotropy ratio defined by κ = K 6 /K 12 . When |κ| > 4, the potential is dominated by the six-fold term, with minima along either the horizontal (κ < −4) or vertical (κ > 4) direction and at every 60 • . This is illustrated in Fig. 1(a). For |κ| ≤ 4, there are 12 minima in the interaction potential as shown in Fig. 1(b). The locations of the minima are given by Six of the minima rotate clockwise as κ is reduced while the other six rotate counterclockwise. Changing the sign of κ is equivalent to a 30 • rotation, which exchanges the locations of the minima and maxima.
It should be noted that the interaction potential in Eq. (1) is not based on an ab-initio derivation, but is instead motivated by an expression for the VL free energy obtained from Ginzburg-Landau (GL) theory 50 . As shown by Olszewski et al., this potential reproduces the experimentally observed macroscopic VL phases in MgB 2 10 and UPt 3 11-13 for a suitable choice of parameters K 6 and K 12 25 . It can thus be considered a minimal model for MD simulations of the VL in superconductors with a hexagonal crystal structure.
In the MD simulations the dynamics of each vortex is governed by an overdamped equation of motion: Here η is the damping constant which is set equal to unity. The force field from the surrounding vortices is given by F vv = −∇(U ) = (−∂U/∂x, −∂U/∂y), yielding components We consider a system of size L × L with L = 108λ, and apply periodic boundary conditions in both the x and y directions. Here λ is the London penetration depth. The sample contains N v = 5130 vortices interacting with N p = 4000 pinning sites placed randomly in non-overlapping positions and generating a pinning force F p i . Each pinning site is modeled as a parabolic trap of range r p = 0.25 that can exert a maximum pinning force of F p = 1.15 on a vortex. The vortex-pin interaction is isotropic, directed towards the center of the pinning site, and given by The values of r p and F p permit at most one vortex to occupy each pinning site, and were chosen in order to stabilize a metastable VL configuration yet still allow the system to be driven towards the equilibrium configuration.

III. RESULTS
The results are organized as follows. First, the protocol used to prepare the metastable VL configuration is described, and how the bulk behavior is characterized by the peak splitting in the structure factor. Second, it is demonstrated how the VL can be driven towards the equilibrium configuration by successive applications of an effective field oscillation, and the dynamics of this process is explored as a function of the oscillation amplitude. Finally, the evolution of VL domain boundaries and energetics is studied at the local scale as the system is driven towards the equilibrium state.
We fix the twelvefold anisotropy amplitude at K 12 = 0.0125 and vary the anisotropy ratio κ by changing K 6 , which was previously established to produce reliable results 25 . Varying κ changes the magnitude of the interaction potential experienced by each vortex, which is equivalent to introducing a change in the effective vortex density. To eliminate this variation, the effective magnetic field, defined by a two-dimensional integral of the interaction potential is kept at a constant value. Here, K 6 = 0 and A v = 2.0 is used as a reference, and the magnitude of the isotropic contribution to the interaction potential is adjusted such Intermediate metastable VL, Fig. 5 Fig. 4(a) FIG. 2. Protocol used to prepare the metastable VL state, plotted as a function of anisotropy ratio κ and temperature F T . From a high temperature molten state at κ = 0, the temperature is decreased to zero, the anisotropy ratio is then increased to κ = 5, and finally the temperature is increased to that A v = 2(2 + K 12 )/(2 + K 6 + K 12 ).

A. Preparing a metastable VL
The protocol used to prepare the VL in a metastable configuration is shown schematically in Fig. 2. The system begins in a molten state with temperature F T = 4.0 and κ = 0. It is then cooled gradually to F T = 0.0 by reducing the temperature by ∆F T = −0.05 every 2 × 10 4 simulation time steps, which is sufficiently slow to ensure that the system reaches an equilibrium state. The anisotropy ratio is then increased to κ = 5 in increments of ∆κ = 0.25 every 2 × 10 4 simulation time steps, while keeping F T = 0. Finally, the temperature is increased to F T = 0.35 using the same rate as for the initial cooling. The final temperature increase provides the system with the necessary thermal energy to overcome activation barriers and be driven towards the equilibrium configuration by the penetration depth oscillations.
The global VL configuration is characterized by computing the structure factor, S(k) ∝ | Nv i exp(−ik ·r i )| 2 . Figure 3(a) shows S(k) at κ = 0 following the initial annealing to F T = 0. Here the intensity is concentrated in 12 (Bragg) peaks lying on a circle, indicating the presence of VL domains that are oriented along one of two different directions corresponding to the minima in the interaction potential in Eq. (1). The radius agrees with the calculated value for a triangular lattice, shows the azimuthal dependence of the structure factor folded into one 60 • segment. This is fitted by a two-peak Voigt function which allows a determination of the peak splitting, ∆θ. A Voigt function is used because it provides overall better fits to the azimuthal intensity distribution than a Gaussian or a Lorentzian.
A metastable VL is established when κ is increased at zero temperature, since without thermal fluctuations the system is no longer able to evolve enough to remain in the changing equilibrium state. This is illustrated in the plots of ∆θ in Fig. 4(a). Although ∆θ decreases it always remains higher than the equilibrium value ∆θ eq = 60 • − 2θ min (κ), indicating that the system becomes trapped in a metastable state. The measured splitting angle has a kink near κ = 4 where ∆θ eq reaches zero. Once the final value of κ = 5 has been reached, the system is allowed to evolve for an additional 2 × 10 4 simulation time steps. This only gives rise to a modest additional reduction of the splitting angle, which stabilizes at a value of ∆θ ∼ 20 • . During the subsequent temperature increase at κ = 5, ∆θ remains largely unchanged as shown in Fig. 3(b), indicating that the thermal fluctua- tions alone are insufficient to permit the system to escape the metastable state.
Real space information about the VL is obtained from Voronoi constructions of the vortex positions, providing local information about domain formation, VL defects and grain boundaries. Figure 5(a) shows the Voronoi construction for an equilibrium VL at κ = 0 and F T = 0, corresponding to Fig. 3. Here the system breaks up into domains with two different orientations. Within each triangular domain, the vortices are sixfold coordinated, but due to the twelvefold anisotropy of the vortex-vortex interaction potential, there are two degenerate orientations. The grain boundaries between domains are decorated by edge dislocations consisting of fivefold and sevenfold coordinated vortices. Figure 5(b) shows the Voronoi construction after plot after the anisotropy ratio has been increased to κ = 5 at zero temperature. Although the number of edge dislocations has decreased, the system still show clear domain formation with easily discernible grain boundaries. Furthermore, the VL domain have rotated towards each other as indicated by a reduced peak splitting in the structure factor. In contrast, the ground state for this κ would be a single domain containing sixfold coordinated vortices. Heating the system to F T = 0.35 at κ = 5 cause some additional annihilation of defects, but the VL remains in a metastable configuration as seen from the Voronoi construction and structure factor in Fig. 5(c). This confirms that thermal effects at this temperature are not strong enough to permit vortices to hop over the activation barrier between the metastable equilibrium states.

B. Effective field oscillations
In the SANS experiments, the transition from the metastable to the equilibrium state is induced by the application of magnetic field oscillations, which change the number of vortices in the sample 16,17 . Changing the number of vortices is not possible in the MD simulations with periodic boundary conditions, and instead an effective change in the vortex density is modeled by by varying the penetration depth λ = λ 0 + ∆λ. This correspond to a change of the magnetic field where Φ 0 is the flux quantum and B 0 is the field for ∆λ = 0. The corresponding amplitude of the field variation, ∆B = B − B 0 , is given by After preparing the metastable state as described in the previous section, effective field oscillations are per-formed by changing λ in steps of 0.05 every 2,000 simulation time steps until the desired value of ∆λ is reached. The cycle is completed by changing the penetration depth back to the original value λ = λ 0 at the same rate. Values of ∆λ between −0.1λ 0 and −0.5λ 0 were studied, corresponding to ∆B/B 0 between 0.23 and 3.0. Figure 5(d) shows the Voronoi construction and structure factor after 1000 penetration depth oscillations were applied to the mestastable state in Fig. 5(c). The number of dislocations are further reduced and the VL domains have rotated closer to the equilibrium orientation as indicated by the smaller structure factor peak splitting. This shows that the penetration depth oscillations are able to mimic the behavior observed in the SANS experiments 16,17 .
One difference between the penetration depth oscillations, used in the MD simulations, and the experimental magnetic field oscillations is that the latter produces a gradient in the vortex density perpendicular to the sample boundary. This small gradient produces a current, which, in turn, causes vortex motion. In comparison, no current applied in the MD simulations since this uses periodic boundary conditions and there is no sample boundary. It is therefore not possible to directly compare the values of the effective field oscillations used in the MD simulations and the much smaller values used in the SANS experiments (∆B/B 0 = 0.001 − 0.003) 17 .

C. Metastable state relaxation dynamics
To study the relaxation dynamics of the metastable VL state MD simulations were performed in which the system was subjected to successive penetration depth oscillations. Simulations were performed for ∆B/B 0 = 0, 0.2346, 0.5625, 1.0408, 1.7778, and 3.0. Ten simulation runs were performed for each value of ∆B/B 0 . Figure 6(a) shows ∆θ determined from the structure factor peak splitting, averaged over the ten simulation runs, versus the cycle count n + 1, where n is the number of oscillations that have been applied. Using n+1 as the abscissa allows the inclusion of the initial metastable state (n = 0). For non-zero ∆B/B 0 an initial rapid drop of ∆θ is observed, with a magnitude that increases with the amplitude of the effective field oscillation. The data for n ≤ 25 is well fitted by the functional form ∆θ = a(n + 1) γ + b, as shown by the full lines in Fig. 6(a). A similar form was also used to fit the experimental SANS data 17 . Fitted values of γ as a function of ∆B/B 0 are shown in Fig. 6(b). For ∆B/B 0 = 0 (no oscillations), ∆θ decreases very slowly due only to thermal effects at F T = 0.35 and γ ∼ 0 although the uncertainty is large. For the lowest non-zero ∆B/B 0 , γ is still very small but the uncertainty is significantly reduced. As the amplitude of the field oscillations is increased, the magnitude of γ also grows until, for the highest value of ∆B/B 0 , γ = −0.5. The overall behavior of γ with increasing ∆B/B 0 is similar to that found in the SANS experiments 17 , where exponents ranging from γ = −0.051 to γ = −0.22 were obtained as shown in Fig. 6(b). It is likely that the largest values of ∆B/B 0 used in the MD simulations exceed the range explored in the SANS in the experiments, but as discussed a direct comparison is not possible. As evident from Fig. 6(a), the power law fits of ∆θ describe the data well for n 300; however, at large n they do not perform as well. Rather, the system exhibits a logarithmic behavior at large n for all values of ∆B/B 0 .
To gain a better understanding of the longer time behaviors, Fig. 7(a) shows ∆θ for ∆B/B 0 = 3.0 versus n+1 on a linear scale for a single simulation run. After the initial smooth decrease, a series of sharp downward jumps or avalanches appears for for n > 100. Figure 7(b) shows the corresponding fraction of sixfold coordinated vortices P 6 versus n + 1. From this it is evident that the downward jumps in ∆θ are accompanied by sudden increases of P 6 , meaning that the avalanches are due to the annihilation of groups of topological defects. This is different from the behavior for n + 1 ≤ 15 where P 6 undergoes a smooth continuous increase.
A more direct demonstration that the jumps in the splitting angle occurs by means of avalanches, in which multiple defect pairs are annihilated, is shown in Fig. 8. This shows a portion of the sample immediately before and after the jump highlighted in the insets of Fig. 7. In the Voronoi construction prior to the jump in Fig. 8(a), eight 5-7 edge dislocations are located near a VL grain boundary within the circled region. After the jump, in Fig. 8(c), about four of dislocations have annihilated. This is also reflected in the local VL energy shown in Figs. 8(b) and 8(d). After the jump, the overall energy of the system has dropped and become more uniform.
In summary, the dynamics can be broken into two parts with an fast initial power law dependence on n followed by a slower universal logarithmic dependence dominated by bursts or avalanches.

IV. DISCUSSION
The decay of a metastable state studied here and in experiments could have implications for many other nonthermal systems in metastable states. One example is granular matter confined in a vertical container. After the grains are initially poured in to the container, they have a fixed volume fraction; however, if the container is subjected to periodic vertical shaking or tapping, the grains become more compressed and the volume fraction V varies as V ∝ 1/(a + b log(n) where n is the number of taps applied 40 . In 2D systems, there have been several studies examining the coarsening dynamics of topological defects in the case where the equilibrium system would be free of defects. Over time, the density of defects decay according to a power law with an exponent in the range 1/4 to 1/3 [28][29][30][31] . The values of the exponents are consistent with some of the values we obtain at shorter times. Studies of systems that form 2D crystals in the absence of quenched disorder show that when quenched disorder is added, numerous defects appear in the lattice that form grain boundaries at zero temperature; however, if an ac driving is applied that shakes the particles back and forth over the substrate, the defects annihilate over time at a rate that becomes proportional to ln(t) at longer times 32,48,51 . If quenched disorder is present, there are generally some defects that remain permanently trapped.
One of the closest systems outside of superconducting vortices to which we can compare our results is the mechanical annealing of grain boundaries in 2D colloidal crystals, where experiments found that the density of grain boundaries decays as a power law as a function of time 51 . It may be possible that the coarsening dynamics has a faster time scale than the glassy dynamics which are typically associated with a ln(t) dependence, so that for the first 1 to 100 cycles the system behaves as if it is coarsening, but at longer times it reverts to ln(t) behavior when the dynamics becomes dominated by hopping or avalanches.
Although our results are not accurate enough to distinguish between the different types of coarsening, at early times they do show that the basic microscopic ingredients for the model we have presented captures the main features found in the experiments. Additionally, our results agree with the vortex-vortex interactions observed in the experiments, where a metastable twelvefold symmetric state is quenched into a sixfold symmetric state and the ac driving progressively reduces the splitting angle.
It could also be possible that in some cases an ordered system in the presence of quenched disorder would show increased defect formation under an applied drive because some particles become trapped while other particles shear past them. This also suggests that there could be an optimal ac amplitude or frequency for creating the maximum amount of dynamic annealing. In addition, introduction of quenched disorder can cause the system to remember the metastable states. If a cyclic driving current were applied, the vortices would move in a periodic manner, and it would be possible to examine whether they return to the same positions at the end of each cycle in a form of return point memory or random organization, which has been studied in numerous periodically driven systems [42][43][44][45][46][47] .

V. SUMMARY
The relaxation of a metastable vortex state was studied numerically for a system quenched from a twelvefold to a sixfold anisotropy of the vortex-vortex interaction. The system is initially in an equilibrium configuration for the twelvefold anisotropy, consisting of triangular VL domains with two different orientations and separated by grain boundaries decorated by edge dislocations. In com-parison, the ground state for the sixfold anisotropy is a single trianglar VL domain. The bulk VL configuration is characterized by the structure factor. Following the quench from the twelvefod to the sixfold anisotropy, the initial twelve equally spaced structure factor peaks forms six pairs with a decreasing splitting. However, the presence of grain boundaries traps the system is a metastable state with a non-zero splitting angle.
The system is driven towards the ground state by a periodic variation of the penetration depth, which mimics the magnetic field oscillations used in SANS experiments [16][17][18] . This causes a gradual decrease of the structure factor peak splitting with successive applications of the penetration depth cycles. Initially, the decrease follows a power law with an exponent that increases with the amplitude of the penetration depth oscillations, in agreement with the SANS experiments performed in a similar regime 17 . As the number of penetration oscillations grows, the decay of the peak splitting becomes logarithmic with a slope that is independent of the amplitude. This indicates a transition from individual defect annihilation at short times to coarsening dynamics at longer times, where the splitting angle decreases in bursts or avalanches that are correlated with the annihilation of groups of topological defects.