Activated vortex lattice transition in a superconductor with combined sixfold and twelvefold anisotropic interactions

,


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 always to 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,12 .In addition, the VL in these materials exhibit significant metastability when cooled or heated across equilibrium phase transitions in a constant field 11,13 .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 4 .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 suppressed 14 .This "work hardening" is speculated to arise from a proliferation of grain boundaries in the VL.Although spatially resolved SANS measurements provide some support for such a scenario 15 , 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 16 and dynamical [17][18][19] properties of vortex matter.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 [20][21][22] .This was used to study systems with combined sixfold and twelvefold anisotropies relevant for MgB 2 21 .By varying the ratio of the six-and twelvefold contributions to the interaction potential, it was possible to reproduce the three different triangular VL phases observed experimentally in this material 10,13 .The VL phases, denoted F, L and I, differ only in their orientation relative to the MgB 2 crystalline axes.Specifically, F and I phases are oriented along the two high symmetry directions within the hexagonal basal plane (a along the nearest neighbor direction and a * along the next nearest neighbor direction), and are connected by an intermediate, continuously rotating L phase 10,23 .The orientation of the VL phases relative to a and a * is shown in Ref. 23, and the directions of the crystalline axes relative to the interaction potential are indicated in Fig. 1(a).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 21 .
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 included 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.We do not observe metastable effects in simulations that does not include pinning, suggesting that at least some weak quenched disorder plays a role in the SANS experiments [13][14][15] .Here, the quenched disorder will decrease thermal hopping effects so that only the mechanical driving from changing the effective field, as discussed below, will play a role.
In the SANS experiments, the transition to the ground state is driven by small-amplitude magnetic field oscillations, leading to a temporary variation in the vortex density akin to a breathing mode 14 .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 rotates towards the ground state F phase.Initially, the angle between the two degenerate L phase domain orientations decreases as a power law with an exponent that increases with the amplitude of the penetration depth cycles, similar to what was observed in the SANS experiments 14 .For the larger amplitudes there is a cross-over to a logarithmic behavior as the number of penetration depth cycles increases.Here, the VL relaxation occurs in bursts or avalanches that are associated with correlated dislocation annihilation events along the grain boundaries.We discuss how these behaviors compare to two-dimensional (2D) coarsening models [24][25][26][27] , and show that the log(t) decay could arise due to coarsening or glassy motion with quenched disorder 28 .
Our approach can be applied to other particlebased systems with anisotropic interactions, including skyrmion lattices [29][30][31] or colloidal particles with anisotropic interactions 32,33 .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 14,15 and in general 34,35 .We discuss the relevance of our results to the decay of metastable states in other systems with periodic perturbations, such as granular matter 36,37 , colloidal systems 28 , magnetic systems with return point memory 38,39 , and other types of periodically sheared systems that relax to a steady state under ac driving [40][41][42] .

II. METHODS
We follow the approach of Olszewski et al. and use a phenomenological vortex-vortex interaction potential that includes six-and twelvefold anisotropies 21 : 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) 43 .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, for a fixed r, 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), and consequently there are two degenerate triangular VL FIG. 1. Vortex-vortex interaction potential given by Eq. (1).Heat maps illustrate the spatial variation of U (r, θ) for κ = 5 (a) and κ = 0 (b).In both cases K12 = 0.0125.Arrows in (a) indicate the directions of the MgB2 in-plane crystalline axes relative to the interaction potential.(c) Angular dependence of U (r, θ) for a fixed radius and different values of κ.For each curve, the corresponding ground state VL phase is indicated on the right.The FL-and LI-transitions occur at κ = ±4.
configurations.The locations of the minima change as κ is varied and are given by Six of the minima rotate clockwise as κ is reduced while the other six rotate counterclockwise, as shown in Fig. 1(c).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 in Eq. ( 1) is not based on an ab-initio derivation, and, within Ginzburg-Landau (GL) theory, it is strictly speaking not possible to formulate an interaction potential for a single vortex 44 .Instead Eq. ( 1) is motivated by an expression for the entire VL free energy obtained from GL theory 45 .As shown by Olszewski et al., this potential reproduces the experimentally observed macroscopic VL phases in MgB 2 10 and UPt 3 11,12 for a suitable choice of parameters K 6 and K 12 21 .It can thus be considered a minimal and conceptually simple phenomenological 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 For both F vv i and F p i the unit of length is λ.The thermal forces F T i are modeled by Langevin kicks with the properties ⟨F T ⟩ = 0 and where k B is the Boltzmann constant.The Langevin kicks are implemented by applying random forces in the x and y directions drawn from a Gaussian distribution with zero mean and a standard deviation of F T .In our dimensionless units F T is the strength of the thermal forces and referred to as the temperature of the system.The values of r p and F p permit at most one vortex to occupy each pinning site, and were chosen by trial and error in order to stabilize a metastable VL configuration yet still allow the system to be driven towards the equilibrium configuration.
In the two-dimensional MD approximation, the vortices are assumed to be stiff lines spaced far enough apart that their cores do not overlap.This is appropriate as MgB 2 is a three dimensional (i.e., not layered) materials with weak pinning and therefore little or no vortex bending is expected.It should also be noted that MD simulations have been used extensively to model the dy-namics of stiff vortices in systems with pinning [46][47][48] .

III. RESULTS
The results are organized as follows.First, the protocol used to prepare the metastable VL configuration is described, as is 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 the effective field oscillation, and the dynamics of this process is explored as a function of the oscillation amplitude.Finally, the evolution of VL grain 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 21 .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 (5) 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 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 F T = 0.35.Labels indicate the locations and/or trajectories in parameter space illustrated in Figs. 3, 4 and 5.

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 with κ = 0 and a temperature F T = 4.0, which is sufficiently high to put it in a molten state.It is then cooled gradually to F T = 0.0 by reducing the temperature by ∆F T = −0.05every 2 × 10 4 simulation time steps, which is sufficiently slow to ensure that the system reaches an equilibrium state 21 .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 effective field oscillations described in Sect.III B.
The global VL configuration is characterized by computing the structure factor, S(k 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,  imuthal intensity distribution than simple Gaussian or Lorentzian.The peak splitting, ∆θ, is defined as the difference between the two fitted peak positions.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 plot 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. 4(b), indicating that the thermal fluctuations 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.Within vs the anisotropy ratio as κ is increased at zero temperature.After reaching κ = 5 the simulation was continued for an additional 2 × 10 4 simulation time steps (shaded area) with κ fixed to allow the system to reach a static configuration.The red line indicates ∆θeq, expected for an equilibrium VL.(b) ∆θ vs temperature under warming with fixed κ = 5.After reaching F T = 0.35, the simulation was continued for an additional 2 × 10 4 simulation time steps (shaded area) with F T fixed to allow the system to reach a static configuration.
each Voronoi cell, all points are closer to a given vortex than to any other.Furthermore, the number of edges yields the coordination of the vortex within the cell.The Voronoi constructions provide local information about domain formation, VL defects and grain boundaries.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 dislocations consisting of fivefold and sevenfold coordinated vortices.
Figure 5(b) shows the Voronoi construction plotted after the anisotropy ratio has been increased to κ = 5 at zero temperature.Although the number of dislocations has decreased, the system still shows clear domain formation with easily discernible grain boundaries.Some grain boundaries have shifted relative to those in Fig. 5(a) and indicated by the black lines.Furthermore, the VL domains have rotated towards each other as indicated by a reduced peak splitting in the structure factor.In other words, the relaxation towards the equilibrium configuration occurs primarily through a continuous rotation of the metastable VL domains rather than the nucleation and growth of separate ground state domains.The ground state for this κ would be a single domain containing sixfold coordinated vortices.Heating the system to F T = 0.35 at κ = 5 causes some additional annihilation of dislocations and grain boundary motion, 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 14 .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 varying the penetration depth λ = λ 0 + ∆λ.This corresponds to a change in both the system size L = 108λ and the range of the vortex-vortex interaction in Eq. ( 4a) and (4b) as r is measured in units of λ.As such, it is equivalent 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 Changing the range of the vortex-vortex interactions by varying λ modifies the force balance between the unevenly spaced vortices near the grain boundaries.This creates creates mechanical instabilities that allow for vortex jumps or rearrangements to occur.It is important to note that whereas varying λ in an experiment is associated with a change in temperature, all the effective field oscillations are performed with a constant F T .An alternative to the effective field oscillations would have been to drive the system would be with an applied current.Such an approach was previously used to study the dynamical annealing of disordered vortex states 49 .
Each effective field oscillation is performed 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.It was verified that this rate is slow enough that the results do not change if the frequency is further reduced.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 reduced and the VL domains have rotated closer to the equilibrium orientation as indicated by the smaller structure factor peak splitting.Furthermore, more grain boundary motion is evident and also some fracturing of the VL into smaller domains.This shows that the penetration depth oscillations are able to move the grain boundaries and thus mimic the behavior observed in the SANS experiments 14 .Notably, while all simulations were performed with a different configuration of the pin locations the system evolve toward the same equilibrium state.
One difference between the penetration depth oscillations, used in the MD simulations, and the experimental magnetic field oscillations is that the latter produces a small gradient in the vortex density perpendicular to the sample boundary.This gradient gives rise to a current, which, in turn, causes vortex motion.In comparison, no current is applied in the MD simulations since the simulations use 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) 14 .

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 penetration depth 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 14 .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 14 , 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 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 dislocations.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 occur by means of avalanches, in which multiple dislocation 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 dislocations are located near a VL grain boundary within the circled region.After the jump, in Fig. 8(c), about four of the 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 a 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 36 .In 2D systems, there have been several studies examining the coarsening dynamics of topological defects in the case where the equilibrium system would be defect free.Over time, the density of defects decays according to a power law with an exponent in the range 1/4 to 1/3 [24][25][26][27] .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 28,40,50 .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 50 .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 capture 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 38,39,41,42 .

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 dislocations.In comparison, 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 twelvefold to the sixfold anisotropy, the initial twelve equally spaced structure factor peaks form six pairs with a decreasing splitting.However, the presence of grain boundaries and weak pinning traps the system in 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 14 .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 14 .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 annihilation of individual dislocation pairs 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 larger groups of dislocations.
metastable VL, Fig. 5(c) κ κ Equilibrium VL, Figs. 3, 5(a) F T change, Fig. 4(b) Intermediate metastable VL, Fig. 5(b) κ change, Fig. 4(a) Figure 3(b) shows the azimuthal dependence of the structure factor folded into one 60 • segment, fitted by a two-peak Voigt function.Voigt fits (a convolution of a Gaussian function and a Lorentzian) are used throughout this work because they were found to provide better fits to the az-

FIG. 3 .
FIG. 3. Bulk characterization of the VL configuration.(a) False color plot of the structure factor amplitude |S(k)| as a function of ky vs kx following an initial annealing to F T = 0 with κ = 0. (b) Azimuthal dependence of the structure factor amplitude |S(θ)| folded into the 60 • segment indicated by a white outline in panel (a).The line is a fit to a two-peak Voigt function.

FIG. 4 .
FIG.4.Evolution of the structure factor peak splitting ∆θ during the preparation of the metastable VL.(a) ∆θ (circles) vs the anisotropy ratio as κ is increased at zero temperature.After reaching κ = 5 the simulation was continued for an additional 2 × 10 4 simulation time steps (shaded area) with κ fixed to allow the system to reach a static configuration.The red line indicates ∆θeq, expected for an equilibrium VL.(b) ∆θ vs temperature under warming with fixed κ = 5.After reaching F T = 0.35, the simulation was continued for an additional 2 × 10 4 simulation time steps (shaded area) with F T fixed to allow the system to reach a static configuration.
Figure 5(a) shows the Voronoi construction for an equilibrium VL at κ = 0 and F T = 0, corresponding to Fig. 3.

FIG. 5 .
FIG. 5.Voronoi constructions of the vortex positions for different system configurations.Red, white, and blue Voronoi cells correspond to five-, six-, and seven-fold coordinated vortices, respectively.The grey and white shading indicate the VL orientation.Inserts show the structure factor amplitude in the 60 • segment of reciprocal space indicated in Fig.3(a).(a) Equilibrium VL for κ = 0 and F T = 0, following an annealing from the initial molten state.Domains are separated by well-defined grain boundaries decorated by 5-7 dislocations.(b) Metastable VL following an increase of the anisotropy ratio to κ = 5 at fixed F T = 0, corresponding to the rightmost data point in Fig. 4(a).Grain boundaries have fewer dislocations and some have shifted slightly.(c) Metastable VL following a subsequent heating to F T = 0.35 with κ = 5, corresponding to the rightmost data point in Fig. 4(b).Only minor changes are observed with respect to the zero temperature state in (b).(d) After the application of 1000 penetration depth oscillations at F T = 0.35 and κ = 5.The number of dislocations is further reduced and more grain boundaries have shifted.Black lines in panels (b), (c) and (d) indicate the location of VL grain boundaries in (a).

HFIG. 6 .
FIG. 6. Metastable state relaxation dynamics.(a) Structure factor peak splitting ∆θ versus the penetration depth cycle count n + 1.The data are an average over ten simulation runs.The solid lines are fits to the form ∆θ = a(n + 1) γ + b for n + 1 ≤ 25 .For large n, all data show a logarithmic behavior as indicated by the dashed line.(b) The exponent γ from the fits in (a) versus ∆B/B0.Also shown for comparison is the exponent versus ac field amplitude Hac/H0, obtained from SANS experiments on MgB2 14 .The line is a guide to the eye.

FIG. 7 .
FIG. 7. A single simulation run showing sudden, stepwise relaxations for n ≳ 200.(a) Splitting angle ∆θ versus n + 1 for ∆B/B0 = 3.0.This shows a series of discontinuous jumps down at large n.(b) Corresponding fraction of sixfold coordinated particles P6.The sudden decreases in the splitting angle are accompanied by increases in P6, corresponding to the annihilation of dislocations.The insets highlight the jump near n + 1 = 450.

FIG. 8 .
FIG. 8. Portion of the sample (a,b) just before and (c,d) after the jump in ∆θ and P6 near n + 1 = 450 in Fig. 7.The circled region indicates the location of a VL grain boundary.(a,c) Voronoi constructions of the vortex positions.Colors are the same as in Fig. 5. (b,d) Corresponding heatmaps showing the local energy of each vortex based on its interactions with the surrounding VL.Yellow and red is higher than the average and blue is lower.The lines give an approximate indication of the local domain orientation.